Magnetically propagating Hund’s exciton in van der Waals antiferromagnet NiPS3

Magnetic van der Waals (vdW) materials have opened new frontiers for realizing novel many-body phenomena. Recently NiPS3 has received intense interest since it hosts an excitonic quasiparticle whose properties appear to be intimately linked to the magnetic state of the lattice. Despite extensive studies, the electronic character, mobility, and magnetic interactions of the exciton remain unresolved. Here we address these issues by measuring NiPS3 with ultra-high energy resolution resonant inelastic x-ray scattering (RIXS). We find that Hund’s exchange interactions are primarily responsible for the energy of formation of the exciton. Measuring the dispersion of the Hund’s exciton reveals that it propagates in a way that is analogous to a double-magnon. We trace this unique behavior to fundamental similarities between the NiPS3 exciton hopping and spin exchange processes, underlining the unique magnetic characteristics of this novel quasiparticle.


INTRODUCTION
Two-dimensional (2D) vdW materials provide an ideal platform for combining strong electronic correlations, lowdimensional magnetism, and weak dielectric screening to realize novel electronic quasiparticles and functionality [1][2][3].Recent years have seen the identification of excitons in a host of closely related vdW compounds such as NiPS 3 , CrSBr, NiI 2 , and MnPS 3 [4][5][6][7].Within this family, the NiPS 3 exciton exhibits several fascinating properties, including strong interactions between the exciton lifetime and magnetic order [4], thickness-dependent properties in the few-layer-limit [8][9][10], coupling between magnetism and exciton polarization [8,11,12], and unconventional exciton-driven metallic behavior [13].These observations suggest that excitons in magnetic vdW materials such as NiPS 3 might have fundamentally different character from other types of exciton, such as the Frenkel, Wannier, and Hubbard varieties.Frenkel and Wannier excitons form via Coulomb interactions between electrons and holes in different Bloch states and propagate according to the detailed form of the band-structure and electron-hole attraction [14].Hubbard excitons, on the other hand, form from strongly correlated many-body states and their propagation is expected to involve the scattering of spin waves [15,16].
In previous works, the NiPS 3 exciton has been described as a "Zhang-Rice" mode [4].This terminology derives from studies of cuprate superconductors, and refers to a specific form of hybridized wavefunctions that have one hole on the transition metal (Ni) site and one hole on the ligand (S) site and describes the "Zhang-Rice exciton" as a transition from a high-spin triplet to a low-spin singlet [17].However, the way the exciton changes with applied magnetic field has been argued to be incompatible with this picture [18].A Zhang-Rice scenario also does not address why the exciton has such a narrow linewidth [8,11,19].The unsettled and possibly unconventional exciton electronic character suggests that the exciton may also propagate in an exotic manner different to regular excitons, but, to date, this has never been measured.Here, we use ultra-high energy resolution RIXS to directly detect the NiPS 3 exciton momentum dispersion and discover it propagates magnetically in a similar way to the doublemagnon excitation.Through detailed analysis of the exciton wavefunction, we further reveal the different interactions involved in its formation and establish that its primary character is that of a Hund's exciton, distinct from the Zhang-Rice and other scenarios.

Electronic character of the exciton
We start by measuring the incident energy dependence of the Ni L 3 -edge RIXS spectrum of NiPS 3 to identify the different spectral features present (see Fig. 1a and Meth-ods).The most intense peaks centered around 1.0, 1.1, and 1.7 eV energy loss are dd excitations in which electrons transition between different Ni 3d orbitals.A remarkably sharp (almost resolution limited) peak is apparent at an energy loss matching the known energy of the exciton at 1.47 eV.This excitation resonates strongly at 853.4 eV and is well separated from other dd and the higher energy charge-transfer excitations.We identify this feature as the NiPS 3 exciton, consistent with previous reports [4].
To facilitate our understanding of the exciton and its interplay with magnetism, we constructed an effective NiS 6 cluster model representing NiPS 3 (See Supplementary Note 1A).Our model includes Coulomb repulsion, Hund's coupling, crystal field, and Ni-S and S-S hopping.As explained in the Methods section, the rich spectrum, including the detailed splitting of the two dd-excitations at 1.0 and 1.1 eV allows us to obtain a well-constrained model Hamiltonian for NiPS 3 (see Supplementary Note 1B and Supplementary Fig. 2).
To better understand the nature of the exciton, we plot several expectation values describing the NiPS 3 wavefunction in Fig. 1d-f, which reveal that NiPS 3 has dominant Hund's rather than Zhang-Rice character.The plotted expectation values include the hole occupations of the Ni 3d and ligand orbitals and the weights of the d 8 , d 9 L, and d 10 L 2 configurations (L stands for a ligand hole) that make up each state.We also calculate the expectation values for the total spin operator squared Ŝ2 , which, for two holes, has a maximum value of S(S + 1) = 2.The ground state is close to this pure high spin state; while the spin is strongly reduced in the exciton.We therefore confirm that the exciton is dominantly a triplet-singlet excitation.In the Zhang-Rice scenario, the leading character of the ground state would be d 9 L. The dominant component is, in fact, d 8 revealing that the state has dominant Hund's character.
The ground state and the exciton wavefunctions obtained from our model are illustrated in Fig. 1g,h and are described in detail in Supplementary Note 2. We see that substantial charge redistribution occurs during exciton formation, which is partly crystal field and partly a Ni-S charge transfer in nature (see Fig. 1d,e).As explained later in the discussion section, the fact that the ground states have dominant Ni character (rather than dominant Zhang-Rice character) plays a leading role in the energy for exciton formation.

Exciton dispersion
Having clarified the character of the exciton, we study its propagation by tuning the incident x-ray energy to the exciton resonance at 853.4 eV and mapping out the in-plane dispersion with high energy resolution (Fig. 2a  and b).The two high-symmetry in-plane reciprocal space directions both exhibit a small upward dispersion away from the Brillouin zone center with similar bandwidths of ∼ 15 meV.This non-zero dispersion suggests that the exciton excites low-energy quasiparticles as it propagates through the lattice.We consequently, mapped out the low energy excitations in Fig. 2c,d.The strongest feature is the magnon, which was found to be consistent with the prediction based on prior inelastic neutron scattering measurements [20,21] (the white line, see Methods).Intriguingly, we also found another broad low-energy dispersive feature at an energy scale roughly twice as large as the magnon peak (the orange dots in Fig. 2c and d).As we justify in detail later, the observed low-energy feature, in fact, corresponds to "double-magnon" excitations, which is a process in which two spins are flipped on each site making up the excitation to create a pair of magnons with the same spin.By fitting the exciton and doublemagnon energy, we see that these two excitations show similar dispersion despite their drastically different energy scales (Fig. 2e and f), indicating that the exciton propagates in a way that is similar to the double-magnon.Since the dispersive effects are subtle, we confirmed the calibration of the spectra by verifying that the magnon energy exactly reproduces the dispersion obtained in prior inelastic neutron scattering experiments [20,21].We similarly confirmed that the Brillouin zone center exciton energy is consistent with values from optical measurements [4,8,11].

Temperature dependence
To substantiate the identity of the low energy excitations we observe in NiPS 3 , we need to consider the RIXS cross-section.Since RIXS is a photon-in photon-out scattering process with each photon carrying one unit of angular momentum, it can couple to processes involving either zero, one, or two spin flips [22][23][24][25][26].While the spin-flip processes are necessarily magnetic, the zero spin flip process can either correspond to a phonon or to a so called "bi-magnon" in which two magnons with opposite spins are created on neighboring sites to make an excitation with a net spin of zero.The most straightforward way to distinguish magnetic and non-magnetic excitations is to measure the dispersion above the Néel order temperature of T N = 159 K, as presented in Fig. 3. Above T N the exciton remains visible, but it becomes weaker and more diffuse compared to the data at 40 K (see Fig. 3a and the linecuts in Supplementary Fig. 10).Consequently, no dispersion is detectable.The double-magnon peak that was observed at 40 K is replaced by a diffuse, over-damped tail of intensity, contrary to what would be expected for a phonon and corroborating its magnetic origin (see Fig. 3b and the linecuts in Supplementary Fig. 9).The residual intensity arises from short-range spin fluctuations, which are expected to persist well above T N in quasi-2D magnets as long as the thermal energy scale is well below the energy scale of the magnetic interactions [27].We also note that optical phonons in the 70-100 meV energy range in NiPS 3 are known to be minimally-dispersive [28], which again suggests a magnetic origin.

Identifying double-magnons through their resonant profile
Having established a magnetic origin for the dispersion, we can use the energy-dependent resonant profile of RIXS to distinguish different magnetic processes and substantiate our assignment of the low-energy feature as the doublemagnon.Figure 4a plots the energy dependence of the measured double-magnon feature alongside the magnon and the exciton.We compare this with calculations of the RIXS cross-section based on our model, which includes processes involving zero, one, and two spin flips denoted by ∆m S = 0, 1, and 2, respectively, as well as coupling to the exciton.∆m S = 0 reflects the cross section for either elastic scattering or a bi-magnon, ∆m S = 1 corresponds to a magnon, and ∆m S = 2 corresponds to the doublemagnon.Notably, since the double-magnon involves an exchange of two units of spin angular momentum, RIXS is especially suitable for detecting this process.The main resonance around 853 eV shows partial overlap between the double-magnon and bi-magnon resonance, so it does not clearly distinguish between them.However, the presence of the excitation at the satellite resonance at 857.7 eV is only compatible with a ∆m S = 2 double-magnon process.We therefore suggest that the dominant character of the excitation around 80 meV is that of double-magnons substantiating our spectral assignment, although we cannot exclude a small sub-leading contribution.Our assignment is also supported by RIXS studies of NiO where the double-magnon was also found to be much more intense than the bi-magnon [22,24,26,29].The satellite resonance is generated by the exchange part of the corevalence Coulomb interaction on the Ni site.This same interaction creates the double-magnon because the angular momentum state of the core-hole needs to vary in the intermediate state in order to allow two subsequent spin-flip processes to occur in the photon absorption and emission processes.The core-valence exchange interaction facilitates this change in the core-hole state by mixing of the core-hole and valence eigenstates.This conclusion is also borne out in studies of NiO where the same process occurs [22,24,26,29].

DISCUSSION
In this work we use high resolution RIXS to assess the formation and propagation of the excitonic state of NiPS 3 .By combining RIXS and exact diagonalization (ED) calculations, we reveal that the primary mechanism behind the exciton formation is the Hund's interaction.As illustrated in Fig. 1d-h the exciton forms from a ground state with dominant d 8 character and involves significant charge transfer and crystal field changes.As such, the state we identify is quite different from prior descriptions of a Zhang-Rice exciton [4].As we discuss later, the distinction between these models is crucial as it corresponds to a different majority component of the wavefunction, different interactions playing the leading role in the exciton energy, and the possibility of realizing a model with physically reasonable parameters.These issues will be central to efforts to manipulate the exciton energy and cross section.We found that the difference in the prior identification of the exciton arises from using an under-constrained model.If one considers just the exciton energy and assumes that Hund's coupling can take any value, there are a range of different Hund's interactions and charge transfer energy parameters that predict a 1.47 eV exciton.If one adds a further constraint that the 1.0, 1.1, and 1.7 eV ddexcitations must be reproduced within an accuracy similar to their width, properly constrained solutions can be identified (see Supplementary Note 1F).Importantly, the solution found here also yields a physically reasonable value for Hund's coupling J H = 1.24 eV corresponding to 87% of the atomic value.This is relevant because the pure triplet and singlet Zhang-Rice components of the wavefunctions are energetically split by rather weak Ni-S exchange pro- cesses, so it is difficult to justify the 1.47 eV energy scale of the exciton within a model with dominant Zhang-Rice character.In Ref. [4] an unphysically large Hund's coupling corresponding to 120% of the atomic value was required.
Based on the wavefunction extraction performed in this study, we can determine which electronic interactions play the leading role in the exciton energy.To do this, we factorized the wavefunctions into the singlet and triplet components of the d 8 , d 9 L, and d 10 L 2 , configurations as described in Supplementary Note 2. We can then compute the contributions of Hund's, charge transfer, and crystal field to each of these components.We find that the primary contribution to the exciton energy comes from singlet-triplet splitting of the d 8 component of the wavefunction, which means that it is best thought of as a Hund's exciton since its energy of formation is mostly driven by Hund's exchange.Our exciton character derived here, also retains partial magnetic character com-ing from a sizable contribution of states beyond the three d-electron, Zhang-Rice, and ligand, singlet states.This means that the exciton is expected to vary with applied magnetic field compatible with recent observations [18].It also implies that future efforts to realize similar symmetry excitons with different energies should target means to modify the on-site Ni Hund's exchange coupling and not the Ni-S exchange processes that would be the leading contribution to the energy of a Zhang-Rice exciton.
Our detection of exciton dispersion in NiPS 3 proves that the exciton is an intrinsic propagating quasiparticle and excludes prior suggestions that the exciton might be a localized phenomenon associated with defects [30].The most common form of exciton propagation in weakly correlated transition metal chalcogenides involves excitations that are composed of bound pairs of specific Bloch states [31,32].The NiPS 3 exciton is quite different from the more conventional excitons since it is bound by the local Hund's interactions described previously, rather than long-range Coulomb attraction.Recent calculations that work for many more conventional excitons, indeed fail to capture our measured exciton dispersion [33].We also note that the two holes in the Ni e g manifold represent what is in some sense the simplest way to realize a tripletsinglet excitation.Consequently, NiPS 3 has relatively few excitations compared to other materials and the exciton is energetically well separated from other transitions, such that it interacts exclusively by magnons and not other types of excitations.The exciton character is also dominated by processes that rearrange spins on the Ni site, rather than moving charge between more extended states, which would tend to reduce the coupling of the exciton to phonons.These factors may contribute to the long lifetime and narrow linewidth of the exciton.
A key result of this study is that the NiPS 3 exciton propagates like the double-magnon, even though the average energies of the exciton and double-magnon differ by more than one order of magnitude.This remarkable similarity can be understood by analyzing the exchange processes involved in the motion of the quasiparticles.We start by considering the spin-superexchange processes involved in exciton motion, finding that the exciton can swap its position with a spin (see Supplementary Note 4B for details).Consequently, when the exciton moves through the antiferromagnetic lattice, it generates a string of misaligned spins.Given that the exciton appears to propagate freely, we should consider processes that heal the misaligned spins in the wake of the exciton, which leads to the image in Fig. 5a that illustrates the spin flips involved in exciton motion.Similar considerations can be applied to the motion of a double-magnon as shown in Fig. 5b.Importantly, both exciton and double-magnon motion involve four spin exchanges.If we consider the sequence of different overlap integrals involved on the Ni and ligand states, the amplitude of the exciton hopping and double-magnon exchange processes are expected to be quite similar (see Supplementary Notes 4A and 4C).These considerations help us rationalize the similarities in the propagation of the two quasiparticles.A simple empirical tight-binding model fit to the exciton dispersion (see Supplementary Note 5) reveals that the third nearest neighbor interaction is the leading term in determining the exciton dispersion, consistent with the third nearest neighbor spin exchange being the dominant term in the spin Hamiltonian [20,21].We note in passing that a similar picture for a well-known free single magnon propagation in an antiferromagnet would require generating two spin exchanges [34].
The coupling between the exciton and magnetism might be relevant to why the 1.47 eV exciton feature is visible in optics experiments.Since the exciton involves transitions between d-orbitals, it is expected to be nominally optically forbidden in a centrosymmetric crystal due to dipole selection rules.However, these rules can be lifted by perturbations that break the Ni-site symmetry, which includes exciton-spin or exciton-lattice interactions.The fact that there is a strong change in the optical cross section through T N [35], whereas the RIXS is only modestly broadened, also supports the interpretation that excitonspin interactions play a key role in the optical cross section.
Overall, our measurements reveal a Hund's excitonic quasiparticle in NiPS 3 that propagates in a similar manner to a two-magnon excitation.Coming years will likely see further instrumental developments that allow RIXS, and exciton microscopy, measurement of NiPS 3 to be extended to the ultrafast pump-probe regimes [36,37].We believe that this has outstanding potential for understanding new means of using magnetic Hund's excitons to realize new forms of controllable transport of magnetic information.

Sample information
NiPS 3 bulk single crystal samples were procured from 2D Semiconductors, which synthesized the crystals by the chemical vapor transport method.The full unit cell of NiPS 3 has a monoclinic symmetry (Space Group C2/m, #12) with lattice parameters a = 5.8 Å, b = 10.1 Å, c = 6.6 Å, and β = 107.0• .We adopted this monoclinic-unit-cell convention, and index reciprocal space using scattering vector Q = (H, K, L) in reciprocal lattice units (r.l.u.).Therefore, the reciprocal lattice vector c * is perpendicular to the ab plane.
Ni 2+ ions in NiPS 3 lie on a honeycomb lattice in the ab plane and form with ABC-type stacking of the layers.Such stacking breaks the three-fold rotational symmetry of the monolayer structure which can be detected by measuring structural Bragg peaks such as (0, 2, 4).NiPS 3 is prone to characteristic twinning involving three equivalent domains rotated by 120 • in the ab plane [20,38].Laboratory single-crystal x-ray diffraction measurements confirmed the presence of three twin domains in our samples.Therefore, measured quantities should be a weighted average over the three twin domains.We included these domain-averaging effects in the ED and magnon energy calculations.The apparent similarity in the measured dispersion between the two distinct directions (i.e., antiferromagnetic across the zig-zag chain and ferromagnetic along the chain) can be ascribed to the domain averaging effect due to the presence of structural twinning in the measured sample.
NiPS 3 orders magnetically below a Néel temperature of T N = 159 K [38].The magnetic unit cell is the same as the structural unit cell.It has a collinear magnetic structure consisting of antiferromagnetically coupled zigzag chains.

RIXS measurements
Ultra-high-energy-resolution RIXS measurements were performed at the SIX 2-ID beamline of the National Synchrotron Light Source II [39].The surface normal of the sample was c * axis (i.e., L direction).The in-plane orientation was determined by Laue diffraction.Then the pre-aligned sample was cleaved with Scotch tape in air to expose a fresh surface and immediately transferred into the RIXS sample chamber.Ni L 3 -edge RIXS measurements were taken with linear horizontal (π) polarization with a scattering plane of either (H0L) or (0KL).The main resonance energy (around 853 eV) is common for Nicontaining compounds [40][41][42] but different from the previous report [4].This difference comes from the absolute energy calibration of the beamline, but does not affect the RIXS measurements and interpretations, which depend only on the relative changes.The spectrometer was operated with an ultrahigh energy resolution of 31 meV fullwidth at half-maximum (FWHM).The temperature of the sample was kept at T = 40 K except for the temperaturedependent measurements.Since the interlayer coupling in NiPS 3 is relatively weak, the dispersion measurement was taken at a scattering angle of 2Θ = 150 • while varying the incident angle of the x-rays.A self-absorption correction [43] was applied to the RIXS spectra, which, however, does not affect peak positions.

Fitting of the RIXS spectra
In order to quantify the dispersion of the excitons and double-magnons, we fitted these peaks in the measured RIXS spectra to extract their peak positions.Although, in principle both the exciton and double-magnon can contain detailed substructure, we found that both features in our spectra can be accurately fit with simple peak shapes.
In the double-magnon region, we used a Gaussian function for the elastic peak, a damped harmonic oscillator (DHO) model for both the magnon and double-magnon peaks, and a constant background.The width of the elastic peak was fixed to the energy resolution, which was determined by a reference measurement on a multilayer heterostructure sample designed to produce strong elastic scattering.The DHO equation for RIXS intensity S(Q, ω) as a function of Q and energy ω is ) where f Q is the undamped energy, χ Q is the oscillator strength, z Q is the damping factor, k B is the Boltzmann constant, and T is temperature.This DHO was then convoluted with a resolution function (a Gaussian function with peak width fixed to the energy resolution) to describe the magnons and double-magnons.We used the fitted value of f Q to represent the magnon or double-magnon peak energy and used z Q to characterize the peak width.
In the exciton region, we resolved both the main exciton and an additional exciton sideband separated by ∼ 40 meV as shown in Supplementary Fig. 3 (e.g., the spectrum at K = −1.22r.l.u.), consistent with a previous report [4].Therefore, we fitted the data with two Voigt peaks with a third order polynomial background.The width of the Gaussian component was fixed to the energy resolution, and we constrained the widths of these two peaks to be the same.For the scans where the minor peak is not obvious, we further fixed the spacing between these peaks to the average value of 40 meV obtained from other scans.

Exact diagonalization RIXS calculations
Our NiPS 3 data were interpreted using standard ED methods for computing the RIXS intensity [23].The Kramers-Heisenberg formula for the cross section was used.This is derived by treating the interaction between the photon and the material within second order perturbation theory (as is required for scattering via an intermediate state resonance).We use the polarization-dependent dipole approximation for the photon absorption and emission interactions and simulate the presence of a core hole in the intermediate state with a core hole potential.In strongly correlated insulators like NiPS 3 , accurate treatment of the electron-electron interactions are particularly important and the brief presence of a core hole means that local processes dominate the scattering.These factors means that cluster approximations are particularly appropriate and are widely used for this reason [23].We therefore perform calculations for a NiS 6 cluster, which can be projected onto Anderson impurity model (AIM) with essentially no loss of accuracy.As explained in detail in Supplementary Note 1, we were able to extract a well-constrained effective model for NiPS 3 from the data.This was used to generate Fig. 1b based on the parameters specified in Tab.I.The detailed definitions of these parameters can be found in Supplementary Note 1.
Since the ED method employed involves directly computing wavefunctions, the model can be used to extract the wavefunctions plotted in Fig. 1g,h as outlined in Supplementary Note 2. All the calculations were done using the open-source software EDRIXS [44].

Magnetic cross-section calculations
With the validated model in hand, it can be used to compute the resonance profile of the magnetic crosssection.A small Zeeman interaction was applied to the total spin angular momentum of the system, serving as the effective molecular magnetic field in the magnetically ordered state.The initially degenerate ground triplet consequently splits into three levels separated based on the spin state.After diagonalizing the Hamiltonian matrix, we used the Kramers-Heisenberg formula in the dipole approximation to calculate the incident-energy dependent RIXS cross section for transitions between the different elements of the triplet.The experimental geometry was explicitly included in the calculations.Three-fold twinning was also accounted for, which in fact has no effect on the final results due to the preserved cubic symmetry.We also note that these results are independent of the magnitude of the applied effective molecular magnetic field since this interaction is much smaller than the splitting between the ground state triplet and the lowest energy dd-excitation.
which includes the on-site energies Êd and Êp for the Ni 3d and S 3p orbitals, respectively.Here, Ûdd and Ûpp are on-site Coulomb interactions and Ûdp are the intersite interactions, and Ûq is the interaction between the core and valence holes.Regarding the intersite hybridization, Tpp indicates p-p hopping among S sites and Tdp stands for p-d hopping between Ni and S sites.The charge-transfer energy ∆ is included implicitly as discussed below.ζ represents spin-orbit coupling, which we include for the core and valence states.
The crystal field associated with the trigonal distortion in NiPS 3 is on the order of only 1 meV [1,2], We therefore assumed cubic symmetry for simplicity.The t 2g − e g -orbital energy splitting is further set by the Ni crystal field, and is specified by 10D q .To describe Ni-S hopping, we use the Slater-Koster parameters of V ppσ , V ppπ , V pdπ , and V pdσ , which denote hopping between Ni d and S p states with either π or σ orbital symmetry.We denote the energy of the S 3p orbitals measured relative to the Ni e g orbitals by ϵ p and neglect the splitting between the 3p σ and 3p π orbitals in the cluster.
The Coulomb interactions in the model are specified using Slater integrals.Following standard methods, we include F 0 dd , F 2 dd , and F 4 dd to describe interactions on Ni 3d orbitals.Similarly, we used F 0 pp and F 2 pp for S 3p orbitals.F 0 dp is the intersite Coulomb interaction between Ni 3d and S 3p orbitals.In the RIXS intermediate states, we used F 0 dp , F 2 dp , G 1  dp , and G 3 dp to describe the Coulomb interactions between Ni 3d and Ni 2p holes.Lastly, we also included spin-orbit coupling terms of Ni 3d orbitals for the initial and intermediate states (ζ i and ζ n ), as well as the much larger core-hole coupling (ζ c ).An inverse core-hole lifetime Γ c = 0.6 eV half-width at half-maximum (HWHM) was used in order to fit the observed width of the resonance and the final state energy loss spectra are broadened using a Gaussian function with a full-width at half-maximum (FWHM) of 0.05 eV, in order to match the observed width of the final states.

B. Determining the cluster model parameters
Although there are several parameters in the original Hamiltonian, they can be constrained by physical considerations and by exploiting the richness of the NiPS 3 RIXS spectra.We start by outlining the different constraints on the parameters: 1. Coulomb interactions: The Slater integrals for Ni 3d Coulomb interactions can be recast into onsite Coulomb repulsion U and Hund's coupling J H . Since we work in the hole language, the on-site Coulomb repulsion for the ligands is not crucial since double hole occupation on ligands is unlikely.Indeed, we get nearly zero probability for the d 10 L 2 (L stands for a ligand hole) configuration in our calculations, as shown in Fig. 1e in the main text.Therefore, F 0 pp and F 2 pp do not influence our conclusions and are fixed to standard values appearing in the literature [3].Due to the distance between the Ni and S atoms the strength of the intersite Coulomb interaction F 0 dp is expected to be much smaller than the Ni on-site Coulomb interactions and thus plays a negligible role.We therefore fix this value to 1 eV, consistent with values often found in the literature [4].The precise value of this interaction is not crucial because its effect on the charge-transfer energy can also be absorbed into the S 3p orbital energy ϵ p .
Supplementary Table II 2. Hopping integrals: Making use of the Slater-Koster scheme, all the hopping integrals can be derived from two parameters, V pdσ and V ppσ .Here, we fixed V pdπ = −V pdσ /2 and V ppπ = −V ppσ /4, which is common for transition metal compounds [5].

Charge-transfer energy:
In a multi-orbital system, this energy is defined as Once the Ni 3d on-site Coulomb interactions and 10D q are chosen, the charge-transfer energy is only affected by the energy difference between the S 2p states and the Ni 3d e g states, i.e. ϵ p (recall we set the latter to zero).
4. Core-hole potential: The core-hole interactions are known to be only weakly screened in the solid state.We initially set the values to 80% of their atomic values and further refined them based on the x-ray absorption spectrum (XAS) (Supplementary Fig. 1) after we determined other above parameters from the RIXS spectra.

Spin-orbit coupling:
The spin-orbit coupling terms for the Ni 3d orbitals are weak and have negligible effects of the spectra.We consequently simply fixed them to their atomic values.The core-hole spin-orbit coupling parameter is adjusted slightly to match the L 2 peak position.
Aside from the core-hole potential, we have 10D q , U , J H , V pdσ , V ppσ , and ϵ p as the only free parameters, which have distinct effects on the RIXS spectra.The energy of the excitation around 1.0 eV is mostly sensitive to J H , while the position of the feature around 1.1 eV is primarily determined by 10D q .The exciton energy is correlated with a combination of J H , 10D q , and U .Once J H and 10D q are determined by the excitations around 1.0 and 1.1 eV, respectively, we can tune the onsite Coulomb repulsion U to match the exciton energy.ϵ p is constrained by the resonant energy as well as the intensity of the exciton, while the hopping integrals can be determined mainly by the resonant energy of satellite peaks.Using the above strategy, we successfully determine these parameters with estimated error bars of ∼ 1 eV for all the Coulomb interactions and ∼ 0.2 eV for hopping integrals and 10D q .The final parameters are provided in Tab.I.

C. Anderson impurity model
The NiS 6 cluster model can be effectively projected to an Anderson impurity model (AIM) using ligand field theory [6,7].The S 3p orbitals form 10 ligand orbitals that have the same symmetry as the Ni 3d orbitals through p-p hopping.The ligand orbitals further hybridize with the Ni 3d orbitals, which can be parameterized based on p-d hopping integrals.Correspondingly, the energies of the ligand orbitals and their hopping integrals with the Ni 3d orbitals can be evaluated as shown in Tab.II.The missing piece is the onsite Coulomb interactions for the ligand orbitals, and we refined these interactions using the same approach as for the cluster model.The full list of the parameters are shown in Tab.III.We will show later that the AIM calculation results are, as expected, essentially identical to the NiS 6 cluster results.

D. Single site atomic model
To test whether a simpler model might be able to capture the relevant physics, we also considered an atomic model composed of a single Ni atom with 10 effective Ni 3d orbitals.In this case, these spin orbitals represent effective hybridized Ni-S orbitals, so they are consequently quite different from the Ni orbitals in the above models.As expected, with appropriate values of crystal field and Hund's coupling, this model has a triplet ground state and a singlet excited Supplementary Table III.Full list of parameters used for the AIM in the ED calculations.Here, ζL is the spin-orbit coupling strength for the ligand orbitals, U dL is intersite Coulomb interactions, and F 0 LL , F state.We refined the parameters using a methodology similar to the one we used for our cluster model.However, we found that the exciton energy could not be accurately reproduced with this model, which we will discuss in more depth when we compare the models.The obtained values of these parameters are listed in Tab.IV.

E. Discussion of an alternative model
Reference [8] (which we will refer to as Ref. "A") used a NiS 6 cluster to calculate the RIXS spectra as we did here but adopting different parameters.Most of the parameters are listed in that paper but two of them are not, namely the energy different between Ni e g and S 3p orbitals and the core-hole potential parameter F 0 dp , which can be indirectly inferred from ∆ = ϵ(d 9 L) − ϵ(d 8 ) = 0.95 eV and ∆ ′ = ϵ(d 10 cL) − ϵ(d 9 c) = −0.55eV where c stands for a core hole.It can be easily concluded that ϵ p − ϵ d (e g ) = 6.37 eV but it is tricky to obtain F 0 dp since it depends on how core-hole potential is included in ∆ ′ .We choose it to be 7.79 eV.Similarly, we inferred that a similar (but not identical) inverse core-hole lifetime of 0.4 eV HWHM was used and their spectra were broadened using a Lorentzian function with FWHM of 0.1 eV.Table V lists all the parameters.We have confirmed that our approach reproduces Ref. A's results [8] when we adopt their parameters.
Supplementary Table V. Full list of parameters used for the NiS6 cluster model in Ref.A [8].The italic values are parameters that are not directly listed in the paper but inferred from other parameters.Units are eV.

On-site orbital energies and spin-orbit coupling
On  therefore define Similar states exist with two ligand holes as ( Using this basis, the ground state triplet can be described in terms of the d 8 triplet state mixed with a Zhang-Rice triplet, with only a small contribution from other states (which we denote | . . .⟩).The wavefunctions are The wavefunction of the exciton can be expressed as which demonstrates that the exciton is primarily composed of the singlet components of the d 8 , d 9 L, and d 10 L 2 states.We also see a noticeable increase in the ligand character compared to the ground state.A further component arises from mixing with a large number of other states.This wavefunction analysis can further be used to understand the interactions underlying the exciton energy.While the Zhang-Rice singlet and triplet are only split by quite weak exchange interactions, the exciton state has an appreciable fraction of doubly occupied holes, which involve much stronger Hund's exchange.To more accurately quantify the contribution from Hund's interaction J H , we compute the expectation value of the operator ĴH for the exciton wavefunction, i.e., ⟨E| ĴH |E⟩.The calculated expectation value is 1.60 eV, which is indeed the leading factor to set the exciton energy scale and compensated slightly by the charge-transfer process to give the exciton energy of 1.47 eV.consequently use the magnon peak position to perform the fine alignment of the spectra.The expected magnon peak energy is calculated using spin wave theory [12] based on the weighted sum of the magnon branches predicted in the Hamiltonian obtained in Ref. [13].In this process, structural domain averaging introduces an error of ∼ 1 meV and the errors of the neutron results themselves are ∼ 2 meV [13,14].For the spectra at |H(K)| ≤ 0.1, because the magnon and elastic peaks are too close to each other, to make the fit converge, we have to manually fix the energy zero to the value that gives the lowest χ 2 .In this case, we assign an estimated error bar of 5 meV for the energy zero because the fits would severely deviate from the measured lineshapes if we shifted the energy zero by ±5 meV or more.The error bars for the fitted double-magnon and exciton peak positions shown in Fig. 2 include not only their own fitting error (0.5-1 meV for the exciton and ∼ 3 meV for the double-magnon) but also the fitting error of the energy zero (∼ 1 meV) as well as the uncertainty of the calculated magnon energies (model differences and the twinning effect as explained above).

B. Best fits to the low temperature momentum-dependent RIXS spectra
We present the best fits to the linecuts of the RIXS spectra measured at various in-plane momentum transfer at T = 40 K in Supplementary Figs.3-5.Particularly for the low-energy region, we display the three components used in the fits to clarify the identify of different spectral features.

C. Cross-checks of the exciton and double-magnon dispersions
To verify the fitted exciton dispersion, we performed three tests.The first one is to check the consistency with the reported zone-center exciton energy from optical measurements.The exciton energy from previous photoluminescence and optical absorption spectra studies is 1.475 eV with uncertainty below 1 meV [8], which is in line with our fitting results.Even after taking into account the uncertainty from the energy calibration of our spectrometer (∼ 2.5 meV at this energy), this zone-center exciton energy taken from optical measurements still has smaller error bar than our fitting results in this region and can help to see an energy difference between the zone-center exciton and excitons at higher Q in Fig. 2. Since the fitted exciton energy heavily depends on the fitting quality of the low-energy region, we next inspect the fitted elastic and magnon peak intensities.As shown in Supplementary Fig. 6, the fitted magnon intensity matches calculated values quite well.The elastic peak becomes stronger near the zone center as expected for specular reflection.The last check we performed is to test the null hypothesis in which we presume that the exciton is, in fact, non-dispersive, and that the apparent dispersion comes from calibration errors.We find that this indeed leads to an unphysical result in which the elastic line intensity drops at the specular (0, 0) position, when it would be expected to increase or stay the same (Supplementary Fig. 6).
Similarly, we also did the null hypothesis test for the double-magnons.Even though the double-magnon peak is broader than the exciton peak and therefore has larger fitting error for the peak positions, we can still exclude the null hypothesis after carefully inspecting the fitted curves.This is most evident in the spectra near the Brillouin zone center, where the fits assuming non-dispersive double-magnons clearly deviate from the best fits and fail to describe the experimental data (Supplementary Fig. 7).

D. Magnon and double-magnon fits with Voigt functions
Although the DHO model are generally used to fit low-energy magnetic excitations in RIXS spectra, to avoid uncertainties inherent in the model we used, we also investigate the Voigt model to fit the magnons and double-magnons.The fitting result turns out to be quite similar to the original DHO fits (Supplementary Fig. 8).Since the magnon peak is quite narrow, both DHO and Voigt models give nearly identical fitting results.Therefore, the change to the energy zero calibration is minimal, so as to the fitted exciton peak positions.For the double-magnon peak, Voigt function gives similar fitting quality but consistently lower peak positions (by ∼ 2 meV) as a consequence of increased peak width.However, the existence of the double-magnon dispersion and its resemblance to the exciton dispersion are still valid in the Voigt model fits.
Despite that the low-temperature data set is unable to distinguish the two models, the DHO model is more appropriate to fit the high-temperature data set, where the anti-stoke peaks on the negative energy loss side become more apparent (Supplementary Fig. 9).For consistency, we therefore adopt the DHO model throughout the manuscript to fit the magnon and double-magnon peaks.

E. Fits to the high-temperature spectra
Fitting the high-temperature spectra at 190 K is more difficult since the magnon peak is softened and cannot be used for energy zero determination.Thanks to the enhanced elastic peak, we are able to fit the spectra at several large Q positions (e.g., H = 0.69 r.l.u.) where the elastic, magnon, and double-magnon peaks are well separated.The magnon peak is not only softened, but also broadened compared to the low temperature data.We then fix its damping factor z Q to an average value of 16 meV.Supplementary Fig. 9 is the best fits we obtained for the low-energy region, although the error bars for the energy zero determination could be as large as ∼ 5 meV.Then we fit the exciton peak as shown in Supplementary Fig. 10, which has a clear softening and broadening.However, the existence of a dispersive mode is ambiguous here due to the large error bars on the fitted peak positions.

F. Fits to the incident energy dependent spectra
To quantify the resonant behaviors of excitons and double-magnons, we also fit the incident energy dependence of the RIXS spectra to extract their spectral weights (integrated intensities).As expected for a Raman-like process in RIXS the peak energies were seen to be independent of incident energy, so this was used as an additional fitting constraint.In Supplementary Fig. 11, we can see that the magnon peak intensity has only one maximum near the main resonance peak around 853 eV.On the contrary, both the double-magnon peak (Supplementary Fig. 12) and the exciton peak (Supplementary Fig. 13) have two maxima, one near the main resonance peak and the other around 857.7 eV.Such distinct resonance behaviors have been summarized in Fig. 4a.

Supplementary Note 4. Double-magnon and exciton propagation
This section provides further analysis of the propagation of the exciton and double-magnon excitations illustrated schematically in Fig. 5 of the main text.For clarity, we will refer to double-magnon propagation via exchange and exciton propagation via hopping.

A. Double-magnon propagation
NiPS 3 exhibits antiferromagnetic order on a honeycomb lattice with both easy-plane and easy-axis anisotropy in addition to first (J 1 ), second (J 2 ), and third (J 3 ) nearest-neighbor isotropic exchange interactions, with the latter being the largest interaction [13,14].Despite these complexities, we can obtain a significant amount of insight by considering a simplified picture, which nonetheless captures the essentials of the interactions at play.Propagation within NiPS 3 's zig-zag antiferromagnetic ground state can be conceptualized as propagation along either a ferro-or antiferro-magnetic chain direction.We simplify matters further by considering only nearest-neighbor exchange processes in the discussion that follows.(The same procedure applies to longer range exchange processes, as outlined below, so the effects of these terms can also be deduced easily).
For a ferromagnetic chain, the double-magnon propagation occurs via the exchange of |1, 1⟩ and |1, −1⟩, where the states are labeled as |S, m S ⟩ following Fig. 5 of the main text.The amplitude for this process can be calculated using second-order perturbation theory in the spin flip processes, with the intermediate state being nearest-neighbor single magnon states with an energy cost proportional to the sum of the single-ion anisotropy ∝ D and the spin exchange ∝ J 1 .Altogether, we obtain that such a double-magnon can propagate along the ferromagnetic direction at the scale of the order of ∝ J 2 1 /(2D + nJ 1 ) ∼ J 1 , where n is the number of broken magnetic bonds in the intermediate state of perturbation theory.Although such a calculation is in principle only valid in the Ising-like limit, Ref. [15] has shown that it can be extended also to the isotropic case.
The other propagation direction in the NiPS 3 plane is the antiferromagnetic direction (see Fig. 5 of the main text).For this case, we need to invoke a fourth order process in the spin flip terms with three intermediate states, with the highest energy being equivalent to the cost of having two 'extra' double-magnons in a line.This multiplicity is due to the fact that the double-magnon propagation is a composition of two processes along the antiferromagnetic direction, namely 1) the exchange of the |1, 1⟩ and |1, −1⟩ states (illustrated in row 2 and 3 of Fig. 5b of the main text), and 2) propagation of a double-magnon via an intermediate state with two 'extra' double-magnons on the nearest-neighbor sites.A similar case for the magnon propagation in the antiferromagnet was discussed in Ref. [16].Due to spin exchanges entering both the denominator and numerator of the perturbative formulae, such a process again leads to an effective propagation at the scale of the spin exchange J 1 .
These estimations for the exchange process have implicitly assumed that the double-magnon is a bound state, which is formally justified only in the limit of large Ising anisotropy.For NiPS 3 , the double-magnon probably has high decay rates into two nearest-neighbor single-magnon states.Fortunately, the latter should remain bound (due to attractive interactions between magnons on nearest neighbor sites), so our analysis should remain a good order-of-magnitude estimate.
We note that the leading J 3 exchange process connects only antiferromagnetically aligned Ni spins, which would suppress any difference between propagation along different directions in the lattice.This is in addition to the fact that NiPS 3 is structurally and magnetically twinned, meaning that the ferro-and antiferro-magnetic directions in the lattice are not empirically distinguishable (see the methods section).Although technically very challenging, the development of ultrahigh energy resolution RIXS under strain could be implemented to study specific magnetic monodomains to more directly test whether the double-magnon is or is not a bound state.

B. Exciton hopping
Before discussing the amplitude of the exciton hopping, we first consider possible spin exchanges for the |1, 0⟩ i state with |1, ±1⟩ j between sites i and j, following the discussion in Ref. [17].In this case, the relevant terms in the Hamiltonian are where J i,j is the exchange coupling between sites i and j and S ± i are the raising and lowering spin operators on site i.We now turn to the exciton hopping process.Since NiPS 3 is a magnetically ordered insulator, an exciton at site i, denoted by |0, 0⟩ i , can only hop via an interchange process with the S = 1 states at a neighboring site j (denoted here as |1, ±1⟩ j ).The relevant terms in the Hamiltonian for exciton hopping are similar to Eq. ( 8), i.e.: with the prefactors in Eqs. ( 8) and ( 9) being identical.The underlying reason for this key observation is that the dominant processes involved are (super)exchange processes that take place on the ligand orbitals rather than the nickel orbitals [17].As a result, individual s = 1/2 spin flips cannot happen on the nickel atoms, leading to vanishing amplitudes for terms that would differentiate between the |0, 0⟩ singlet hopping and the |1, 0⟩ triplet exchange.

C. Exciton propagation
We start by discussing exciton propagation along the ferromagnetic direction in the zigzag antiferromagnetic ground state.In fact, once the exciton hopping is known [see Eq. ( 9) above], this propagation is the easiest to understand: it merely amounts to the free hopping of the exciton in the ferromagnetic background with a hopping amplitude equal to the spin exchange, similar to the double-magnon case.
A more complex situation is encountered for the antiferromagnetic direction.Here exciton propagation is also due to the hopping process described by Eq. ( 9) but obtaining a coherent propagation is slightly intricate-as shown in Fig. 5a of the main text.This whole process can essentially be divided into two steps.First, the exciton interchanges twice with the spin background, just as in the ferromagnetic background (second and third row of Fig. 5a of the main text).This leads to the creation of an intermediate state with two 'extra' double-magnons situated next to each other.Second, the double-magnons can be annihilated by two spin-exchange processes, in a similar manner as described in Supplementary Note 4 A above.In total, it requires four spin exchanges for the exciton to freely move to the next-nearest-neighbor site, just as the double-magnon does.Such exciton propagation is also at the energy scale proportional to the spin exchange.
Altogether, we observe strong similarities between the way the exciton and the double-magnon move through the spin S = 1 zigzag antiferromagnetic honeycomb lattice-in particular, both motions require the same number of spin exchanges and the energy scales are in both cases the same and proportional to the spin exchanges.As discussed in the following section, the dominant exciton hopping is through the third nearest neighbors connected by antiferromagnetically aligned spins, therefore, the dispersions along the in-plane H and K directions would be expected to be overall rather similar, just like the double-magnon case.Supplementary Note 5. Tight-binding model fit to the exciton dispersion Tight-binding model approaches have been widely used to model exciton dispersion in molecular solids [18].Although the way the Hund's exciton we studied here moves is different to the mechanism for Frenkel excitons (spin exchange processes as discussed above instead of dipole-dipole interactions), tight-binding models can still be useful at a phenomenological level to provide a simple, but informative empirical approach to extract the lengthscale of the effective interactions governing the exciton dispersion.
Using the monoclinic unit-cell notation, we formulated a simplified effective tight-binding model on a two-dimensional honeycomb lattice with isotropic effective "hopping" terms t n .We obtain three terms in the exciton dispersion (E t1 , E t2 , and E t3 ) associated with first, second, and third nearest neighbor effective interactions (t 1 , t 2 and t 3 ), where By co-fitting the measured low-temperature exciton dispersion along both H and K directions, we obtain the results as shown in Supplementary Fig. 14.For E t1 and E t3 , we selected the sign of the solution based on the empirically observed upward-dispersion trend near the Brillouin zone center.We tested the individual contribution of the three different functional forms and found that the fits with third nearest neighbor effective hopping alone can capture the observed periodicity of the dispersion, indicating that third nearest neighbor interactions play the leading role in the exciton dispersion.As stated before, this phenomenological model is an effective parameterization of a process that fundamentally arises from magnetic exchange and not real hopping.Our observation of leading third nearest neighbor interactions is consistent with the fact the third nearest neighbor spin exchange is dominant in the spin Hamiltonian of NiPS 3 [13,14].Energy loss (meV) Supplementary Figure 9. RIXS spectra measured in the (H0L) plane at T = 190 K with an energy window chosen to isolate the magnon and double-magnon.Each panel displays the spectra at a specific in-plane momentum transfer H measured in the (H0L) plane at T = 190 K, above its magnetic ordering temperature.We used linear horizontal π polarization of the incident x-rays at the resonant energy of 853.4 eV for excitons.These data are the same as the intensity maps shown in Fig. 3b and are provided to show the linecuts directly.The solid red lines are best fits to the data with three components, i.e., the elastic line (gray), the magnon peak (blue), and the double-magnon peak (orange).These are provided to clarify the identity of different spectral features.The blue dashed line shows fits to the 40 K data for comparison.The vertical dashed line in each panel labels the energy zero.Error bars represent one standard deviation.

Fig. 1 .
Fig. 1.Electronic character of the NiPS3 exciton.a, RIXS intensity map as a function of incident photon energy through the Ni L3 resonance.The exciton is visible at an energy loss of 1.47 eV and reaches a maximum intensity at an incident energy of 853.4 eV.These data were taken at 40 K with π-polarized x-rays incident on the sample at θ = 22.6 • and scattered to 2Θ = 150 • .b, RIXS calculations for NiPS3 that capture the energy and resonant profile of the dd-transitions and exciton in the material.c, Calculated unbroadened RIXS intensity (vertical lines) and broadened RIXS spectra (solid curve) at the main resonant incident energy of the exciton peak (i.e., Ei = 853.4eV).d-f, Description of the ground and excited states in NiPS3.d shows the hole occupations of Ni 3d (denoted by d) and ligand (denoted by L) orbitals.e displays probabilities of having d 8 , d 9 L, and d 10 L 2 configurations.f gives the expectation value of the total spin operator squared ⟨ Ŝ2 ⟩.The orange (green) vertical lines in d-f indicate the energy for the double-magnons (excitons).g,h, Wavefunction illustrations extracted from b for g the exciton and h the ground state.The size of each orbital (3d for the central Ni site and 3p for the six neighboring S sites) is proportional to its hole occupation.The color represents the expectation value of the spin operator along the z axis ⟨ Ŝz⟩, again calculated separately for the Ni and S states.Therefore, the change in spin state and the partial transfer of holes involved in the exciton transition is encoded in the change in color and size of orbitals, respectively.We represent the ground state by only the down-spin configuration, omitting the up-spin and spin-zero elements of the triplet.

Fig. 2 .
Fig.2.Low temperature exciton dispersion and comparison with double-magnons.a,b, RIXS intensity maps as a function of the H and K in-plane momentum transfer, respectively, with an energy window chosen to isolate the exciton dispersion.The overlaid green squares mark the peak positions of the exciton.c and d show the low energy dispersion at equivalent momenta with the observed inelastic feature, including magnons (white lines) and double-magnons (orange circles).Panels e and f, show that both the exciton and double-magnon have similar dispersion with an energy offset of ∼ 1.4 eV.All the measurements were taken at T = 40 K using π-polarized incident x-rays at an incident energy of 853.4 eV corresponding to the exciton resonance.The asterisks in panels e and f denote the reported exciton energy from optical measurements[4, 8, 11]   with error bars from our instrument energy calibration (one standard deviation).All other error bars are 1-σ confidence intervals evaluated from the fitting as explained in the Methods.Detailed linecuts showing the fitting are provided in Supplementary Figs.3-5.

Fig. 3 .
Fig. 3. High temperature exciton dispersion and comparison with double-magnons above the Néel transition.The RIXS intensity maps with energy window chosen to isolate the excitons (a) and magnons/double-magnons (b) as a function of in-plane momentum transfer H measured in the (H0L) scattering plane.All the measurements were taken at T = 190 K with linear horizontal π polarization of the incident x-rays at the resonant energy of 853.4 eV for excitons.The same data are provided as linecuts in Supplementary Figs. 9 and 10.

2 Fig. 4 .
Fig. 4. Resonance behavior of magnons, doublemagnons and excitons.a, The measured RIXS spectral weights of the magnon, double-magnon and exciton extracted by fitting experimental RIXS spectrum for each incident energy.Error bars represent one standard deviation.Data were taken at 40 K at a scattering angle of 2Θ = 150 • and an incident angle of θ = 22.6 • .b, The calculated RIXS spectral weights of the exciton and low-energy zero-, one-, and twospin-flip transitions (∆ms = 0, 1, 2) as a function of incident energy.The curves in both panels are scaled for clarity.

Fig. 5 .
Fig. 5. Illustration of exciton and double-magnon propagation based on perturbation theory.The top row represents the antiferromagnetic background and subsequent rows show the time evolution of the state.a, After the singlet |0, 0⟩ exciton forms (second row from the top) it exchanges spin with neighboring sites such that it moves while flipping spins and breaking magnetic bonds; free propagation to the next nearest neighbor site (bottom row) is possible after four spin exchanges, involving up to four magnons created in the intermediate state (middle rows).b, After the double-magnon excitation is created on the same site (second row from the top), it can freely move to the next nearest neighbor (bottom row) by four spin exchanges and exciting four magnons in the intermediate state (middle rows).The similarities between the propagation in a and b rationalize the experimentally observed similar dispersion relation of the exciton and double-magnon.These processes are mediated by the different spin-exchange interactions, with the third nearest neighbor exchange process playing the leading role.

Supplementary Figure 4 . 10 fFigure 7 .Supplementary Figure 8 .
RIXS spectra measured in the (H0L) plane at T = 40 K with an energy window chosen to isolate the magnon and double-magnon.Each panel displays the spectra at a specific in-plane momentum transfer H measured in the (H0L) plane.The solid lines are best fits to the data with three components, i.e., the elastic line (gray), the magnon peak (blue), and the double-magnon peak (orange).The detailed description of the fitting can be found in the Methods Section and Supplementary Note 3. The vertical dashed line in each panel labels the energy zero.Error bars represent one standard deviation.Supplementary Representative fitting comparisons between the best fits (blue lines), which give a dispersive double-magnon mode, and the fits assuming a null hypothesis of a non-dispersive double-magnon (red lines).a-c, RIXS spectra measured in the (H0L) plane at T = 40 K at three representative in-plane momentum transfers H. d-f, RIXS spectra measured in the (0KL) plane at T = 40 K at three representative in-plane momentum transfers K.We can clearly see that red lines assuming non-dispersive double-magnons are worse compared to the best fits (blue lines) which give dispersive double-magnons.The vertical dashed line in each panel labels the energy zero.Error bars represent one standard deviation.Fit comparison between the DHO and Voigt models used for the magnon and doublemagnon peaks.a-b, fitted peak positions of the excitons (green squares) and double-magnons (orange circles) using the DHO model for the magnon and double-magnon peaks.The data is the same as the dispersion plots shown in Fig.2e and f. c-d, fitted peak positions of the excitons (green squares) and double-magnons (orange circles) using the Voigt function for the magnon and double-magnon peaks.The fitting results based on the two different functional forms for the magnons and double-magnons are quite similar, validating the robustness of the exciton and double-magnon dispersions.

Supplementary Figure 14 .
Tight-binding model fits to the exciton dispersion.Each column displays tight-binding model fits with only first (a-b), second (c-d), and third (e-f ) nearest neighbor interactions, respectively.The green squares are the measured exciton dispersion as a function of the H and K in-plane momentum transfer, respectively.The black lines are fitted curves.The best fitted values are t1 = 1.4(5) meV, t2 = −0.5(1)meV, and t3 = 1.7(3) meV, respectively.Error bars represent one standard deviation.

TABLE I .
Full list of parameters used in the AIM calculations.Units are eV.

Table I .
Full list of parameters used for the NiS6 cluster model in the ED calculations.Units are eV.
. Projection of the NiS6 cluster model to the AIM.Here, ϵL refers to the ligand orbital energies determined by the S 3p orbital energy ϵp and Tpp = |Vppσ − Vppπ|.The hopping integrals V between Ni 3d orbitals and ligand orbitals can be evaluated from p-d hoppings.Units are eV.

.
2LL , and F4LL are the Slater integrals for ligand orbital onsite Coulomb interactions.Units are eV.Full list of parameters used for the single site atomic model in the ED calculations.All parameters are in units of eV.