# Supplementary Information for: Magnetically Propagating Hund’s
exciton in van der Waals antiferromagnet NiPS3

W. He<EMAIL_ADDRESS>Y. Shen  K. Wohlfeld  J. Sears  J. Li  J. Pelliciari  M.
Walicki  S. Johnston  E. Baldini  V. Bisogni  M. Mitrano  M. P. M. Dean
<EMAIL_ADDRESS>
###### Contents

  1. Supplementary Note 1 Exact diagonalization calculations
    1. Supplementary Note 1.1 Cluster model
    2. Supplementary Note 1.2 Determining the cluster model parameters
    3. Supplementary Note 1.3 Anderson impurity model
    4. Supplementary Note 1.4 Single site atomic model
    5. Supplementary Note 1.5 Discussion of an alternative model
    6. Supplementary Note 1.6 Comparison of all the models
  2. Supplementary Note 2 Wavefuction analysis
  3. Supplementary Note 3 Detailed fitting procedures for the RIXS spectra
    1. Supplementary Note 3.1 Energy zero determination
    2. Supplementary Note 3.2 Best fits to the low temperature momentum-dependent RIXS spectra
    3. Supplementary Note 3.3 Cross-checks of the exciton and double-magnon dispersions
    4. Supplementary Note 3.4 Magnon and double-magnon fits with Voigt functions
    5. Supplementary Note 3.5 Fits to the high-temperature spectra
    6. Supplementary Note 3.6 Fits to the incident energy dependent spectra
  4. Supplementary Note 4 Double-magnon and exciton propagation
    1. Supplementary Note 4.1 Double-magnon propagation
    2. Supplementary Note 4.2 Exciton hopping
    3. Supplementary Note 4.3 Exciton propagation
  5. Supplementary Note 5 Tight-binding model fit to the exciton dispersion

##  Supplementary Note 1 Exact diagonalization calculations

Here we outline the detailed approach and values used in our resonant
inelastic x-ray scattering (RIXS) calculations. Throughout this section, we
represent the model in hole (rather than electron) language.

###  Supplementary Note 1.1 Cluster model

Our cluster representation of NiPS3 uses six S atoms surrounding a central Ni
atom. The orbital basis includes the full set of ten Ni $3d$ spin-orbitals and
six S $3p$ spin-orbitals for each of the six neighboring S ions. The Ni
$2p_{3/2}$ core states are included to simulate the resonant pathway involved
in RIXS. The Hamiltonian for this model is

$\mathcal{H}=\hat{E}_{d}+\hat{E}_{p}+\hat{U}_{dd}+\hat{U}_{pp}+\hat{U}_{dp}+\hat{U}_{q}+\hat{T}_{pp}+\hat{T}_{dp}+\hat{\zeta},$
(1)

which includes the on-site energies $\hat{E}_{d}$ and $\hat{E}_{p}$ for the Ni
$3d$ and S $3p$ orbitals, respectively. Here, $\hat{U}_{dd}$ and
$\hat{U}_{pp}$ are on-site Coulomb interactions and $\hat{U}_{dp}$ are the
intersite interactions, and $\hat{U}_{q}$ is the interaction between the core
and valence holes. Regarding the intersite hybridization, $\hat{T}_{pp}$
indicates $p$–$p$ hopping among S sites and $\hat{T}_{dp}$ stands for $p$–$d$
hopping between Ni and S sites. The charge-transfer energy $\Delta$ is
included implicitly as discussed below. $\hat{\zeta}$ represents spin-orbit
coupling, which we include for the core and valence states.

The crystal field associated with the trigonal distortion in NiPS3 is on the
order of only 1 meV [1, 2], We therefore assumed cubic symmetry for
simplicity. The $t_{\mathrm{2g}}-e_{\mathrm{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\sigma}$,
$V_{pp\pi}$, $V_{pd\pi}$, and $V_{pd\sigma}$, which denote hopping between Ni
$d$ and S $p$ states with either $\pi$ or $\sigma$ orbital symmetry. We denote
the energy of the S $3p$ orbitals measured relative to the Ni $e_{g}$ orbitals
by $\epsilon_{p}$ and neglect the splitting between the $3p_{\sigma}$ and
$3p_{\pi}$ 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 ($\zeta_{\mathrm{i}}$ and
$\zeta_{\mathrm{n}}$), as well as the much larger core-hole coupling
($\zeta_{\mathrm{c}}$). An inverse core-hole lifetime $\Gamma_{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.

###  Supplementary Note 1.2 Determining the cluster model parameters

Supplementary Table 1: Full list of parameters used for the NiS6 cluster model in the ED calculations. Units are eV. On-site orbital energies and spin-orbit coupling |  | On-site Coulomb interactions
---|---|---
$\epsilon_{d}(e_{\mathrm{g}})$ | $\epsilon_{d}(t_{\mathrm{2g}})$ | $\epsilon_{p}$ | $\zeta_{\mathrm{i}}$ | $\zeta_{\mathrm{n}}$ | $\zeta_{\mathrm{c}}$ |  | $F^{0}_{dd}$ | $F^{2}_{dd}$ | $F^{4}_{dd}$ | $F^{0}_{pp}$ | $F^{2}_{pp}$
0 | 0.42 | 7.5 | 0.083 | 0.102 | 11.4 |  | 7.88 | 10.68 | 6.68 | 3.2 | 5
Hopping integrals |  | Intersite Coulomb interactions |  | Core-hole potential
$V_{pd\sigma}$ | $V_{pp\sigma}$ |  | $F^{0}_{dp}$ | $F^{2}_{dp}$ | $G^{1}_{dp}$ | $G^{3}_{dp}$ |  | $F^{0}_{dp}$ | $F^{2}_{dp}$ | $G^{1}_{dp}$ | $G^{3}_{dp}$
0.95 | 0.8 |  | 1 | 0 | 0 | 0 |  | 7.45 | 6.56 | 4.92 | 2.80

Although there are several parameters in the original Hamiltonian, they can be
constrained by physical considerations and by exploiting the richness of the
NiPS3 RIXS spectra. We start by outlining the different constraints on the
parameters:

  1. 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}\underline{L}^{2}$
($\underline{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 $\epsilon_{p}$.

  2. 2.

Hopping integrals: Making use of the Slater-Koster scheme, all the hopping
integrals can be derived from two parameters, $V_{pd\sigma}$ and
$V_{pp\sigma}$. Here, we fixed $V_{pd\pi}=-V_{pd\sigma}/2$ and
$V_{pp\pi}=-V_{pp\sigma}/4$, which is common for transition metal compounds
[5].

  3. 3.

Charge-transfer energy: In a multi-orbital system, this energy is defined as
$\Delta=\epsilon(d^{9}\underline{L})-\epsilon(d^{8})=[\epsilon_{d}(e_{\mathrm{g}})+\epsilon_{p}+F^{0}_{dp}]-[\epsilon_{d}(e_{\mathrm{g}})\times
2+F^{0}_{dd}-8F^{2}_{dd}/49-F^{4}_{dd}/49]$. 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. $\epsilon_{p}$ (recall we set the latter to zero).

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

  5. 5.

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\sigma}$, $V_{pp\sigma}$, and $\epsilon_{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. $\epsilon_{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 $\sim 1$ eV for all the Coulomb interactions and $\sim 0.2$ eV for hopping
integrals and $10D_{q}$. The final parameters are provided in Tab. 1.

###  Supplementary Note 1.3 Anderson impurity model

The NiS6 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. 2. 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. 3. We will show later
that the AIM calculation results are, as expected, essentially identical to
the NiS6 cluster results.

Supplementary Table 2: Projection of the NiS6 cluster model to the AIM. Here, $\epsilon_{L}$ refers to the ligand orbital energies determined by the S $3p$ orbital energy $\epsilon_{p}$ and $T_{pp}=|V_{pp\sigma}-V_{pp\pi}|$. The hopping integrals $V$ between Ni $3d$ orbitals and ligand orbitals can be evaluated from $p$–$d$ hoppings. Units are eV. Symmetry | $\epsilon_{d}$ | $\epsilon_{L}$ | $V$
---|---|---|---
$e_{\mathrm{g}}$ | $\epsilon_{d}(e_{\mathrm{g}})$ | $\epsilon_{p}-T_{pp}$ | $\sqrt{3}|V_{pd\sigma}|$
$t_{\mathrm{2g}}$ | $\epsilon_{d}(t_{\mathrm{2g}})$ | $\epsilon_{p}+T_{pp}$ | $2|V_{pd\pi}|$
Supplementary Table 3: Full list of parameters used for the AIM in the ED calculations. Here, $\zeta_{\mathrm{L}}$ is the spin-orbit coupling strength for the ligand orbitals, $U_{dL}$ is intersite Coulomb interactions, and $F^{0}_{LL}$, $F^{2}_{LL}$, and $F^{4}_{LL}$ are the Slater integrals for ligand orbital onsite Coulomb interactions. Units are eV. On-site orbital energies |  | Spin-orbit coupling |  |  | Hopping integral
---|---|---|---|---|---
$\epsilon_{d}(e_{\mathrm{g}})$ | $\epsilon_{d}(t_{\mathrm{2g}})$ | $\epsilon_{L}(e_{\mathrm{g}})$ | $\epsilon_{L}(t_{\mathrm{2g}})$ |  | $\zeta_{\mathrm{i}}$ | $\zeta_{\mathrm{n}}$ | $\zeta_{\mathrm{c}}$ | $\zeta_{\mathrm{L}}$ |  |  | $V_{pd\sigma}$
0 | 0.42 | 6.5 | 8.5 |  | 0.083 | 0.102 | 11.4 | 0 |  |  | 0.95
Coulomb interactions |  | Core-hole potential
$F^{0}_{dd}$ | $F^{2}_{dd}$ | $F^{4}_{dd}$ | $F^{0}_{LL}$ | $F^{2}_{LL}$ | $F^{4}_{LL}$ | $U_{dL}$ |  | $F^{0}_{dp}$ | $F^{2}_{dp}$ | $G^{1}_{dp}$ | $G^{3}_{dp}$
7.88 | 10.68 | 6.68 | 0.46 | 2.59 | 1.62 | 1 |  | 7.45 | 6.56 | 4.92 | 2.80

###  Supplementary Note 1.4 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 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. 4.

Supplementary Table 4: Full list of parameters used for the single site atomic model in the ED calculations. All parameters are in units of eV. $\epsilon_{d}(e_{\mathrm{g}})$ | $\epsilon_{d}(t_{\mathrm{2g}})$ | $\zeta_{\mathrm{i}}$ | $\zeta_{\mathrm{n}}$ | $\zeta_{\mathrm{c}}$ | $F^{0}_{dd}$ | $F^{2}_{dd}$ | $F^{4}_{dd}$ | $F^{0}_{dp}$ | $F^{2}_{dp}$ | $G^{1}_{dp}$ | $G^{3}_{dp}$
---|---|---|---|---|---|---|---|---|---|---|---
0 | 1.07 | 0.083 | 0.102 | 11.2 | -0.70 | 5.26 | 3.29 | 0.26 | 6.18 | 2.89 | 1.65

###  Supplementary Note 1.5 Discussion of an alternative model

Reference [8] (which we will refer to as Ref. “A”) used a NiS6 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_{\mathrm{g}}$ and S $3p$ orbitals
and the core-hole potential parameter $F^{0}_{dp}$, which can be indirectly
inferred from $\Delta=\epsilon(d^{9}\underline{L})-\epsilon(d^{8})=0.95$ eV
and
$\Delta^{\prime}=\epsilon(d^{10}\underline{c}\underline{L})-\epsilon(d^{9}\underline{c})=-0.55$
eV where $\underline{c}$ stands for a core hole. It can be easily concluded
that $\epsilon_{p}-\epsilon_{d}(e_{\mathrm{g}})=6.37$ eV but it is tricky to
obtain $F^{0}_{dp}$ since it depends on how core-hole potential is included in
$\Delta^{\prime}$. 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 5 lists all the parameters. We have confirmed that our
approach reproduces Ref. A’s results [8] when we adopt their parameters.

Supplementary Table 5: 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-site Coulomb interactions
---|---|---
$\epsilon_{d}(e_{\mathrm{g}})$ | $\epsilon_{d}(t_{\mathrm{2g}})$ | $\epsilon_{p}$ | $\zeta_{\mathrm{i}}$ | $\zeta_{\mathrm{n}}$ | $\zeta_{\mathrm{c}}$ |  | $F^{0}_{dd}$ | $F^{2}_{dd}$ | $F^{4}_{dd}$ | $F^{0}_{pp}$ | $F^{2}_{pp}$
0 | 1 | $\mathit{6.37}$ | 0.08 | 0.08 | 11.5 |  | 8 | 14.68 | 9.12 | 1 | 0
Hopping integrals |  | Intersite Coulomb interactions |  | Core-hole potential
$V_{pd\sigma}$ | $V_{pd\pi}$ |  | $F^{0}_{dp}$ | $F^{2}_{dp}$ | $G^{1}_{dp}$ | $G^{3}_{dp}$ |  | $F^{0}_{dp}$ | $F^{2}_{dp}$ | $G^{1}_{dp}$ | $G^{3}_{dp}$
0.74 | -0.4 |  | 0 | 0 | 0 | 0 |  | $\mathit{7.79}$ | 6.18 | 4.63 | 2.63

###  Supplementary Note 1.6 Comparison of all the models

Supplementary Fig. 2 summarizes the calculation results for all the three
models, as well as the one from Ref. A [8]. The exciton at around $1.47$ eV
can only be captured in the models with ligand orbitals. Compared with the
combination of parameters used in Ref. A [8], our model (NiS6 cluster or
equivalently AIM) better captures the $dd$ excitations at 1.0 and 1.1 eV. We
also see that our model has the secondary satellite resonance at 857.7 eV,
consistent with the data, whereas the model of Ref. A [8] predicts this
feature around 857 eV. Moreover, our model is compatible with the previous
optical work, which reported absorption peaks at $1.1$, $1.7$, $2.2$, $3.5$,
and $4.6$ eV [9]. On the contrary, the calculated spectra using the model in
Ref. A [8] cannot reproduce the energies of the last three excitations well.

If one considers the exciton in isolation, without requiring a precise
reproduction of the shape of the $dd$-excitations, our model and Ref. A [8]
both obtain 1.47 eV consistent with the experimental energy of this mode. This
observation suggests that there is some ambiguity in determining the model
parameters from the exciton energy alone. By analyzing the different models,
we found that this ambiguity arises from an anticorrelation between Hund’s
interactions and charge transfer energy, which produces a ranges of parameters
that obtain the observed exciton energy. We note, however, that the Hund’s
energy is known to be weakly screened in solids and generally falls in a
reasonably narrow range of values around 70-90% of the values calculated for
isolated atoms [10]. In fact, the onsite Coulomb interaction parameters for Ni
$3d$ orbitals in Ref. A [8], $F^{2}_{dd}$ and $F^{4}_{dd}$, are 120% of atomic
values, which is unphysical. Our parameters are 87% of atomic values, which
falls within the normal range for correlated insulators.

Although the ground state and exciton in all these models share similar
overall symmetry, their electronic composition varies appreciably. For
example, our charge-transfer energy ($2.50$ eV) is larger than that reported
in Ref. A ($0.95$ eV) [8]. Correspondingly, the exciton in our model, which is
a singlet, has $1.26$ versus $0.74$ holes occupying the Ni and S sites,
respectively. Conversely, it has $1.08$ versus $0.92$ holes occupying the Ni
and S sites, respectively, in the model of Ref. A [8]. As noted in the main
text, the preferential Ni occupation of the exciton is vital for understanding
the interactions underlying the exciton formation.

##  Supplementary Note 2 Wavefuction analysis

As discussed in the main text, we performed a detailed wavefunction analysis
for the AIM model to investigate the nature of the exciton, but results are
similar for both the AIM and cluster models.

Figure 1d–f plots several expectation values for the NiPS3 wavefunction, which
reveals a significant portion of doubly occupied $d^{8}$ configuration and
deviation from the pure $d^{9}\underline{L}$ states. A full understanding of
these effects can be obtained by projecting the real wavefunctions onto those
devised by Zhang and Rice to describe the $d^{9}\underline{L}$ manifold [11].
The Zhang-Rice singlet can be written as

$\displaystyle|\mathrm{ZRS}\rangle$
$\displaystyle=\frac{1}{2}(|d_{x^{2}-y^{2},\downarrow}L_{x^{2}-y^{2},\uparrow}|-|d_{x^{2}-y^{2},\uparrow}L_{x^{2}-y^{2},\downarrow}|+|d_{z^{2},\downarrow}L_{z^{2},\uparrow}|-|d_{z^{2},\uparrow}L_{z^{2},\downarrow}|).$
(2)

The Zhang-Rice triplet states are

$\begin{split}|\mathrm{ZRT},S_{z}=1\rangle&=\frac{1}{\sqrt{2}}(|d_{x^{2}-y^{2},\uparrow}L_{z^{2},\uparrow}|-|d_{z^{2},\uparrow}L_{x^{2}-y^{2},\uparrow}|)\\\
|\mathrm{ZRT},S_{z}=0\rangle&=\frac{1}{2}(|d_{x^{2}-y^{2},\uparrow}L_{z^{2},\downarrow}|+|d_{x^{2}-y^{2},\downarrow}L_{z^{2},\uparrow}|-|d_{z^{2},\uparrow}L_{x^{2}-y^{2},\downarrow}|-|d_{z^{2},\downarrow}L_{x^{2}-y^{2},\uparrow}|)\\\
|\mathrm{ZRT},S_{z}=-1\rangle&=\frac{1}{\sqrt{2}}(|d_{x^{2}-y^{2},\downarrow}L_{z^{2},\downarrow}|-|d_{z^{2},\downarrow}L_{x^{2}-y^{2},\downarrow}|).\end{split}$
(3)

Since $d^{8}$ states have a significant amount of weight in the wavefunction,
these states must also be considered. We therefore define

$\begin{split}|^{1}A_{1},d\rangle&=\frac{1}{\sqrt{2}}(|d_{x^{2}-y^{2},\uparrow}d_{x^{2}-y^{2},\downarrow}|+|d_{z^{2},\uparrow}d_{z^{2},\downarrow}|)\\\
|^{3}A_{2},d,S_{z}=1\rangle&=|d_{x^{2}-y^{2},\uparrow}d_{z^{2},\uparrow}|\\\
|^{3}A_{2},d,S_{z}=0\rangle&=\frac{1}{\sqrt{2}}(|d_{x^{2}-y^{2},\downarrow}d_{z^{2},\uparrow}|+|d_{x^{2}-y^{2},\uparrow}d_{z^{2},\downarrow}|)\\\
|^{3}A_{2},d,S_{z}=-1\rangle&=|d_{x^{2}-y^{2},\downarrow}d_{z^{2},\downarrow}|.\end{split}$
(4)

Similar states exist with two ligand holes as

$\begin{split}|^{1}A_{1},L\rangle&=\frac{1}{\sqrt{2}}(|L_{x^{2}-y^{2},\uparrow}L_{x^{2}-y^{2},\downarrow}|+|L_{z^{2},\uparrow}L_{z^{2},\downarrow}|)\\\
|^{3}A_{2},L,S_{z}=1\rangle&=|L_{x^{2}-y^{2},\uparrow}L_{z^{2},\uparrow}|\\\
|^{3}A_{2},L,S_{z}=0\rangle&=\frac{1}{\sqrt{2}}(|L_{x^{2}-y^{2},\downarrow}L_{z^{2},\uparrow}|+|L_{x^{2}-y^{2},\uparrow}L_{z^{2},\downarrow}|)\\\
|^{3}A_{2},L,S_{z}=-1\rangle&=|L_{x^{2}-y^{2},\downarrow}L_{z^{2},\downarrow}|.\end{split}$
(5)

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 $|\dots\rangle$). The
wavefunctions are

$\begin{split}|\mathrm{GS},S_{z}=-1\rangle&=0.757\
|^{3}A_{2},d,S_{z}=-1\rangle+0.163\ |^{3}A_{2},L,S_{z}=-1\rangle+0.625\
|\mathrm{ZRT},S_{z}=-1\rangle+0.099\ |\dots\rangle\\\
|\mathrm{GS},S_{z}=0\rangle&=0.757\ |^{3}A_{2},d,S_{z}=0\rangle+0.163\
|^{3}A_{2},L,S_{z}=0\rangle+0.625\ |\mathrm{ZRT},S_{z}=0\rangle+0.099\
|\dots\rangle\\\ |\mathrm{GS},S_{z}=1\rangle&=0.757\
|^{3}A_{2},d,S_{z}=1\rangle+0.163\ |^{3}A_{2},L,S_{z}=1\rangle+0.625\
|\mathrm{ZRT},S_{z}=1\rangle+0.099\ |\dots\rangle.\end{split}$ (6)

The wavefunction of the exciton can be expressed as

$|\mathrm{E}\rangle=0.430\ |^{1}A_{1},d\rangle+0.209\
|^{1}A_{1},L\rangle+0.766\ |\mathrm{ZRS}\rangle+0.430\ |\dots\rangle,$ (7)

which demonstrates that the exciton is primarily composed of the singlet
components of the $d^{8}$, $d^{9}\underline{L}$, and $d^{10}\underline{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
$\hat{J}_{H}$ for the exciton wavefunction, i.e., $\langle
E|\hat{J}_{H}|E\rangle$. 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.

##  Supplementary Note 3 Detailed fitting procedures for the RIXS spectra

As noted in Methods Section, we use Voigt functions to fit the excitons, and
use a damped harmonic oscillator (DHO) convoluted with a Gaussian resolution
function to fit the magnon as well as the double-magnon peaks. In this
section, we provide further details on our fitting procedures.

###  Supplementary Note 3.1 Energy zero determination

To ensure a reliable measurement of the subtle exciton dispersion we observe
here, we require an accurate calibration of the energy zero, which is often
done by fitting the elastic line. However, in the current case of NiPS3 we
only observe quite weak elastic scattering and it is only possible to
precisely fit the elastic line in a subset of the spectra. We 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 $\sim 1$ meV and the errors of the neutron results themselves are $\sim 2$
meV [13, 14]. For the spectra at $|H(K)|\leq 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 $\chi^{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 $\pm 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 $\sim 3$ meV for
the double-magnon) but also the fitting error of the energy zero ($\sim 1$
meV) as well as the uncertainty of the calculated magnon energies (model
differences and the twinning effect as explained above).

###  Supplementary Note 3.2 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.

###  Supplementary Note 3.3 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 ($\sim
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).

###  Supplementary Note 3.4 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 $\sim 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.

###  Supplementary Note 3.5 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 $\sim 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.

###  Supplementary Note 3.6 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.

###  Supplementary Note 4.1 Double-magnon propagation

NiPS3 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 [14, 13]. 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 NiPS3’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\rangle$ and $|1,-1\rangle$, where the states are labeled as
$|S,m_{S}\rangle$ 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 $\propto D$ and the spin exchange $\propto J_{1}$. Altogether, we
obtain that such a double-magnon can propagate along the ferromagnetic
direction at the scale of the order of $\propto J_{1}^{2}/(2D+nJ_{1})\sim
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 NiPS3 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\rangle$ and
$|1,-1\rangle$ 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 NiPS3, 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 NiPS3 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.

###  Supplementary Note 4.2 Exciton hopping

Before discussing the amplitude of the exciton hopping, we first consider
possible spin exchanges for the $|1,0\rangle_{\bf i}$ state with $|1,\pm
1\rangle_{\bf j}$ between sites $\bf{i}$ and $\bf{j}$, following the
discussion in Ref. [17]. In this case, the relevant terms in the Hamiltonian
are

$\begin{split}h^{\text{spin-flip}}&=\frac{J_{{\bf i},{\bf
j}}}{2}\left(S^{+}_{\bf i}S^{-}_{\bf j}+\text{h.c.}\right)\\\ &=J_{{\bf
i},{\bf j}}\big{(}|1,1\rangle_{\bf j}|1,0\rangle_{\bf i}~{}_{\bf j}\langle
1,0|_{\bf i}\langle 1,1|+|1,-1\rangle_{\bf i}|1,0\rangle_{\bf j}~{}_{\bf
i}\langle 1,0|_{\bf j}\langle 1,-1|+\text{h.c.}\big{)},\end{split}$ (8)

where $J_{{\bf i},{\bf j}}$ is the exchange coupling between sites $\bf{i}$
and $\bf{j}$ and $S^{\pm}_{\bf i}$ are the raising and lowering spin operators
on site ${\bf i}$.

We now turn to the exciton hopping process. Since NiPS3 is a magnetically
ordered insulator, an exciton at site ${\bf i}$, denoted by $|0,0\rangle_{\bf
i}$, can only hop via an interchange process with the $S=1$ states at a
neighboring site ${\bf j}$ (denoted here as $|1,\pm 1\rangle_{\bf j}$). The
relevant terms in the Hamiltonian for exciton hopping are similar to Eq. (8),
i.e.:

$\displaystyle h^{\text{hop}}=J_{{\bf i},{\bf j}}\big{(}|1,1\rangle_{\bf
i}|0,0\rangle_{\bf j}~{}_{\bf i}\langle 0,0|_{\bf j}\langle
1,1|+|1,-1\rangle_{\bf i}|0,0\rangle_{\bf j}~{}_{\bf j}\langle 1,-1|_{\bf
i}\langle 0,0|+\text{h.c.}\big{)},$ (9)

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\rangle$ singlet hopping and the
$|1,0\rangle$ triplet exchange.

###  Supplementary Note 4.3 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.1 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

$\begin{split}E_{t1}=&t_{1}\sqrt[\pm]{4\cos^{2}(\pi H)+4\cos(\pi H)\cos(\pi
K)+1}\\\ E_{t2}=&2t_{2}\\{\cos(2\pi H)+\cos[\pi(H+K)]+\cos[\pi(H-K)]\\}\\\
E_{t3}=&t_{3}\sqrt[\pm]{4\cos^{2}(2\pi H)+4\cos(2\pi H)\cos(2\pi
K)+1}.\end{split}$ (10)

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 NiPS3 [13, 14].

Supplementary Figure 1: XAS spectra. The spectra was taken at 40 K using
$\pi$-polarized x-rays in the total fluorescence yield mode. The sample
geometry was the same as the measurement for the RIXS energy map (i.e.,
$\theta=22.6^{\circ}$, $2\Theta=150^{\circ}$, $(0KL)$ scattering plane) except
that the photo diode detector was placed slightly above the scattering plane.
Supplementary Figure 2: Comparison of RIXS energy maps calculated using
different models. a–e, The measured and calculated RIXS energy maps. f–j,
Zoom-in of the measured and calculated RIXS maps in the region outlined by the
dashed box in panels a–e. The calculations assume a constant lifetime for the
excitations, and omit the continuum of states needed to accurately simulate
the fluorescent processes that occur at high energies, so the transitions
above 2 eV appear far narrower in the calculations compared to the theory.
Supplementary Figure 3: RIXS spectra measured at $\bm{T=40}$ K with an energy
window chosen to isolate the exciton. a, The exciton peak at various in-plane
momentum transfer $H$ measured in the ($H0L$) plane. b, The exciton peak at
various in-plane momentum transfer $K$ measured in the ($0KL$) plane. All the
measurements were taken with linear horizontal $\pi$ polarization of the
incident x-rays at the resonant energy of $853.4$ eV for excitons. The data is
the same as the intensity maps shown in Fig. 2a and b. The solid lines are
fits to the data with details in Methods Section. Data are shifted vertically
for clarity. Error bars represent one standard deviation. Supplementary Figure
4: RIXS spectra measured in the ($\bm{H0L}$) plane at $\bm{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 Figure 5: RIXS spectra measured in the ($\bm{0KL}$)
plane at $\bm{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 $K$ measured in the ($0KL$) 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 Figure 6: Fitting result
comparison between the best fit (blue circle), which gives a dispersive
exciton mode, and the fit assuming a null hypothesis of a non-dispersive
exciton (red square) coming from a miscalibration of the spectra. a–b,
Integrated Intensity of the magnon peak measured in ($H0L$) plane a and
($0KL$) plane b, respectively. The black solid lines are the calculated magnon
intensities as described in the Methods Section, with the same scaling factor
for both panels. c–d, Integrated Intensity of the elastic peak measured in
($H0L$) plane c and ($0KL$) plane d, respectively. We see that the best fit is
in good accord with the calculated magnon intensity, whereas the null
hypothesis can be excluded as it implies an unphysical decrease in the elastic
line intensity at $(0,0)$, which is the point at which the elastic line
intensity would usually be maximized. Supplementary Figure 7: 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. Supplementary Figure 8: Fit comparison between the DHO and Voigt
models used for the magnon and double-magnon 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 9: RIXS spectra measured in the ($\bm{H0L}$)
plane at $\bm{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 $\pi$ 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. Supplementary Figure 10: RIXS
spectra measured at $\bm{T=190}$ K with an energy window chosen to isolate the
exciton. The measurements were taken at various in-plane momentum transfer $H$
in the ($H0L$) plane at $T=190$ K, above its magnetic ordering temperature. We
used linear horizontal $\pi$ 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. 3a and are provided to show the linecuts
directly. The solid lines are fits to the data. The dashed lines show the
fitted curves to the corresponding 40 K data for comparison. Data are shifted
vertically for clarity. Error bars represent one standard deviation.
Supplementary Figure 11: Incident energy dependence of the RIXS spectra with
an energy window chosen to highlight the intensity change of the magnon peak.
Each panel displays the spectra with a specific incident photon energy through
the Ni $L_{3}$ resonance measured in the ($0KL$) plane at $T=40$ K with linear
horizontal $\pi$ polarization of the incident x-rays. These data are the same
as the intensity maps shown in Fig. 1a 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). The vertical dashed line in each panel labels the energy zero.
Error bars represent one standard deviation. Supplementary Figure 12: Incident
energy dependence of the RIXS spectra with an energy window chosen to
highlight the intensity change of the double-magnon peak. Each panel displays
the spectra with a specific incident photon energy through the Ni $L_{3}$
resonance measured in the ($0KL$) plane at $T=40$ K with linear horizontal
$\pi$ polarization of the incident x-rays. These data are the same as that
shown in Supplementary Fig. 11 but zoomed in on the double-magnon peak
(orange). Error bars represent one standard deviation. Supplementary Figure
13: Incident energy dependence of the RIXS spectra with an energy window
chosen to highlight the intensity change of the exciton peak. The measurements
were taken at various incident photon energy through the Ni $L_{3}$ resonance
measured in the ($0KL$) plane at $T=40$ K with linear horizontal $\pi$
polarization of the incident x-rays. These data are the same as the intensity
maps shown in Fig. 1a and are provided to show the linecuts directly. The
solid lines are the fits to the data with details in Methods Section. Data are
shifted vertically for clarity. Error bars represent one standard deviation.
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 $t_{1}=1.4(5)$ meV, $t_{2}=-0.5(1)$ meV, and
$t_{3}=1.7(3)$ meV, respectively. Error bars represent one standard deviation.

## References

  * [1] Afanasiev, D. _et al._ Controlling the anisotropy of a van der Waals antiferromagnet with light.  _Science Advances_ 7, eabf3096 (2021).
  * [2] Chandrasekharan, N. & Vasudevan, S.  Magnetism and exchange in the layered antiferromagnet NiPS3.  _Journal of Physics: Condensed Matter_ 6, 4569 (1994).
  * [3] Wang, Y., Fabbris, G., Dean, M. & Kotliar, G.  EDRIXS: An open source toolkit for simulating spectra of resonant inelastic x-ray scattering.  _Computer Physics Communications_ 243, 151–165 (2019).
  * [4] Neudert, R. _et al._ Four-band extended Hubbard Hamiltonian for the one-dimensional cuprate ${\mathrm{Sr}}_{2}{\mathrm{CuO}}_{3}$: Distribution of oxygen holes and its relation to strong intersite Coulomb interaction.  _Phys. Rev. B_ 62, 10752–10765 (2000).
  * [5] Mizuno, Y. _et al._ Electronic states and magnetic properties of edge-sharing Cu-O chains.  _Phys. Rev. B_ 57, 5326–5335 (1998).
  * [6] Eskes, H., Tjeng, L. H. & Sawatzky, G. A.  Cluster-model calculation of the electronic structure of CuO: A model material for the high-${T}_{c}$ superconductors.  _Physical Review B_ 41, 288–299 (1990).
  * [7] Haverkort, M. W., Zwierzycki, M. & Andersen, O. K.  Multiplet ligand-field theory using Wannier orbitals.  _Physical Review B_ 85, 165113 (2012).
  * [8] Kang, S. _et al._ Coherent many-body exciton in van der Waals antiferromagnet NiPS3.  _Nature_ 583, 785–789 (2020).
  * [9] Kim, S. Y. _et al._ Charge-Spin Correlation in van der Waals Antiferromagnet NiPS3.  _Physical Review Letters_ 120 (2018).
  * [10] Khomskii, D.  _Transition metal compounds_ (Cambridge University Press, 2014).
  * [11] Zhang, F. C. & Rice, T. M.  Effective Hamiltonian for the superconducting Cu oxides.  _Physical Review B_ 37, 3759–3761 (1988).
  * [12] Toth, S. & Lake, B.  Linear spin wave theory for single-Q incommensurate magnetic structures.  _Journal of Physics Condensed Matter_ 27 (2015).
  * [13] Scheie, A. _et al._ Spin wave hamiltonian and anomalous scattering in ${\mathrm{NiPS}}_{3}$.  _Phys. Rev. B_ 108, 104402 (2023).
  * [14] Wildes, A. R. _et al._ Magnetic dynamics of NiPS3.  _Physical Review B_ 106, 174422 (2022).
  * [15] Zheng, W., Oitmaa, J. & Hamer, C. J.  Series studies of the spin-$\frac{1}{2}$ Heisenberg antiferromagnet at $T=0$: Magnon dispersion and structure factors.  _Phys. Rev. B_ 71, 184440 (2005).
  * [16] Verresen, R., Pollmann, F. & Moessner, R.  Quantum dynamics of the square-lattice Heisenberg model.  _Phys. Rev. B_ 98, 155102 (2018).
  * [17] Takubo, K. _et al._ Unusual superexchange pathways in an NiS2 triangular lattice with negative charge-transfer energy.  _Physical Review Letters_ 99, 037203 (2007).
  * [18] Cudazzo, P., Sottile, F., Rubio, A. & Gatti, M.  Exciton dispersion in molecular solids.  _Journal of Physics: Condensed Matter_ 27, 113204 (2015).


N^{2}t}$. It follows from the triangle inequality and the above bound that

$\sup_{S,S^{\prime}}\mathbb{E}_{\xi_{[T]}}\left[\left|D(\bar{w}_{T},W^{*})-D(\bar{w}^{\prime}_{T},W^{*})\right|\right]\leq\sup_{S,S^{\prime}}\mathbb{E}_{\xi_{[T]}}\left[\|\bar{w}_{T}-\bar{w}^{\prime}_{T}\|\right]\leq\frac{20\sqrt{2LM}(1+\log(T))}{\lambda\rho
nT}.$

Since $\xi_{[T]}$ are independent on $S$, as a direct consequence of applying
McDiarmid’s inequality with $c_{i}\equiv
c=\frac{20\sqrt{2LM}(1+\log(T))}{\lambda\rho nT}$ to
$h(S):=D(\bar{w}_{T},W^{*})$, we can show that with probability at least
$1-\delta$ over the randomness of $S$,

$\mathbb{E}_{\xi_{[T]}}\left[D(\bar{w}_{T},W^{*})-\mathbb{E}_{S}\left[\mathbb{E}_{\xi_{[T]}}\left[D(\bar{w}_{T},W^{*})\right]\right]\right]\leq
c\sqrt{\frac{nT\log(1/\delta)}{2}}=\frac{20\sqrt{LM\log(1/\delta)}(1+\log(T))}{\lambda\rho\sqrt{nT}}.$

We next derive a bound for $\mathbb{E}_{S}\left[D(\bar{w}_{T},W^{*})\right]$.
In view of Jensen’s inequality and the quadratic growth property of $F$ we
have

$\displaystyle\mathbb{E}_{S}\left[\mathbb{E}_{\xi_{[T]}}\left[D(\bar{w}_{T},W^{*})\right]\right]=$
$\displaystyle\mathbb{E}_{\xi_{[T]}}\left[\mathbb{E}_{S}\left[D(\bar{w}_{T},W^{*})\right]\right]$
$\displaystyle\leq$
$\displaystyle\mathbb{E}_{\xi_{[T]}}\left[\sqrt{\mathbb{E}_{S}\left[D^{2}(\bar{w}_{T},W^{*})\right]}\right]$
$\displaystyle\leq$
$\displaystyle\mathbb{E}_{\xi_{[T]}}\left[\sqrt{\frac{2}{\lambda}\mathbb{E}_{S}\left[R(\bar{w}_{T})-R^{*}\right]}\right]$
$\displaystyle\lesssim$
$\displaystyle\mathbb{E}_{\xi_{[T]}}\left[\sqrt{\frac{\rho\left[R(w_{0})-R^{*}\right]}{\lambda
T^{2}}+\frac{L}{\lambda^{2}\rho nT}R^{*}+\frac{\sqrt{\epsilon}}{\lambda
T^{2}}\left(\frac{L}{\lambda\rho}+G\sqrt{\frac{1}{\lambda\rho}}\right)}\right]$
$\displaystyle=$
$\displaystyle\sqrt{\frac{\rho\left[R(w_{0})-R^{*}\right]}{\lambda
T^{2}}+\frac{L}{\lambda^{2}\rho nT}R^{*}+\frac{\sqrt{\epsilon}}{\lambda
T^{2}}\left(\frac{L}{\lambda\rho}+G\sqrt{\frac{1}{\lambda\rho}}\right)},$

where in the last inequality we have invoked Theorem 4. Therefore, based on
the previous two inequalities we obtain that with probability at least
$1-\delta$ over $S$,

$\displaystyle\mathbb{E}_{\xi_{[T]}}\left[D(\bar{w}_{T},W^{*})\right]$
$\displaystyle\lesssim$
$\displaystyle\frac{\sqrt{LM\log(1/\delta)}\log(T)}{\lambda\rho\sqrt{nT}}+\sqrt{\frac{\rho\left[R(w_{0})-R^{*}\right]}{\lambda
T^{2}}+\frac{L}{\lambda^{2}\rho nT}R^{*}+\frac{\sqrt{\epsilon}}{\lambda
T^{2}}\left(\frac{L}{\lambda\rho}+G\sqrt{\frac{1}{\lambda\rho}}\right)},$

which gives the desired bound. ∎

###  C.3 Proof of Theorem 7

Here we prove the following restated Theorem 7. See 7

We need the following lemma essentially from Bousquet et al. (2020, Corollary
8) that gives a near-tight generalization bound for a learning algorithm that
is uniformly stable with respect to loss function.

######  Lemma 12 (Bousquet et al. (2020)).

Suppose that a learning algorithm $\mathcal{A}_{w}$, parameterized by $w$,
satisfies
$|\ell(\mathcal{A}_{w_{S}}(x),y)-\ell(\mathcal{A}_{w_{S^{\prime}}}(x),y)|\leq\varrho$
for any $(x,y)\in\mathcal{X}\times\mathcal{Y}$ and $S\doteq S^{\prime}$.
Assume the loss function satisfies $0\leq\ell(y^{\prime},y)\leq M$ for all
$y,y^{\prime}$. Then for any $\delta\in(0,1)$, with probability at least
$1-\delta$ over an i.i.d. data set $S$ of size $N$,

$\left|R(\mathcal{A}_{w_{S}})-R_{S}(\mathcal{A}_{w_{S}})\right|\lesssim\varrho\log(N)\log\left(\frac{1}{\delta}\right)+M\sqrt{\frac{\log\left(1/\delta\right)}{N}}.$

With this lemma in place, we can prove the main result in Theorem 7

###### Proof of Theorem 7.

Let $S=\\{S_{t}\\}_{t\in[T]}$ and $S^{\prime}=\\{S^{\prime}_{t}\\}_{t\in[T]}$
be two sets of data minibatches satisfying $S\doteq S^{\prime}$. Note that
$\gamma_{t}\equiv\gamma=\sqrt{\frac{T}{n}}$. Then according to Proposition 1
the average output $\bar{w}_{T}$ and $\bar{w}^{\prime}_{T}$ respectively
generated by Algorithm 1 over $S$ and $S^{\prime}$ satisfy

$\sup_{S,S^{\prime}}\|\bar{w}_{T}-\bar{w}^{\prime}_{T}\|\leq\frac{4\sqrt{2LM}}{n\gamma}+\sum_{t=1}^{T}2\sqrt{\frac{2\epsilon_{t}}{\gamma}}\leq\frac{5\sqrt{2LM}}{n\gamma}=\frac{5\sqrt{2LM}}{\sqrt{nT}}.$

where in the last but one inequality we have used the condition
$\epsilon_{t}\leq\frac{LM}{4nT^{2}\sqrt{N}}$. It follows that

$|\ell(\bar{w}_{T};z)-\ell(\bar{w}^{\prime}_{T};z)|\leq\sqrt{2ML}\|\bar{w}_{T}-\bar{w}^{\prime}_{T}\|\leq\frac{10LM}{\sqrt{nT}},$

where we have used $\ell(\cdot;\cdot)$ is $\sqrt{2LM}$-Lipschitz with respect
to its first argument (which is implied by Lemma 2). In view of Assumption 3
we have

$|r(\bar{w}_{T})-r(\bar{w}^{\prime}_{T})|\leq
G\|\bar{w}_{T}-\bar{w}^{\prime}_{T}\|\leq\frac{5G\sqrt{2LM}}{\sqrt{nT}}.$

This preceding two inequalities indicate that M-SPP is
$\frac{10LM+5G\sqrt{2LM}}{\sqrt{nT}}$-uniformly stable with respect to the
composite loss function $\ell+r$. By invoking Lemma 12 to M-SPP we obtain that

$\left|R(w_{S})-R_{S}(w_{S})\right|\lesssim\frac{(LM+G\sqrt{LM})\log(nT)}{\sqrt{nT}}\log\left(\frac{1}{\delta}\right)+M\sqrt{\frac{\log\left(1/\delta\right)}{nT}}.$

The proof is concluded. ∎

## References

  * Agarwal et al. (2012) Alekh Agarwal, Peter L. Bartlett, Pradeep Ravikumar, and Martin J. Wainwright.  Information-theoretic lower bounds on the oracle complexity of stochastic convex optimization.  _IEEE Trans. Inf. Theory_ , 58(5):3235–3249, 2012.
  * Allen-Zhu (2017) Zeyuan Allen-Zhu.  Katyusha: The first direct acceleration of stochastic gradient methods.  _J. Mach. Learn. Res._ , 18:221:1–221:51, 2017.
  * Anitescu (2000) Mihai Anitescu.  Degenerate nonlinear programming with a quadratic growth condition.  _SIAM J. Optim._ , 10(4):1116–1135, 2000.
  * Asi and Duchi (2019a) Hilal Asi and John C Duchi.  The importance of better models in stochastic optimization.  _Proceedings of the National Academy of Sciences_ , 116(46):22924–22930, 2019a.
  * Asi and Duchi (2019b) Hilal Asi and John C. Duchi.  Stochastic (approximate) proximal point methods: Convergence, optimality, and adaptivity.  _SIAM J. Optim._ , 29(3):2257–2290, 2019b.
  * Asi et al. (2020) Hilal Asi, Karan N. Chadha, Gary Cheng, and John C. Duchi.  Minibatch stochastic approximate proximal point methods.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , virtual, 2020.
  * Bach and Moulines (2011) Francis R. Bach and Eric Moulines.  Non-asymptotic analysis of stochastic approximation algorithms for machine learning.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 451–459, Granada, Spain, 2011.
  * Bartlett et al. (2005) Peter L Bartlett, Olivier Bousquet, and Shahar Mendelson.  Local rademacher complexities.  _The Annals of Statistics_ , 33(4):1497–1537, 2005.
  * Bassily et al. (2020) Raef Bassily, Vitaly Feldman, Cristóbal Guzmán, and Kunal Talwar.  Stability of stochastic gradient descent on nonsmooth convex losses.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , virtual, 2020.
  * Bertsekas (2011) Dimitri P. Bertsekas.  Incremental proximal methods for large scale convex optimization.  _Math. Program._ , 129(2):163–195, 2011.
  * Bottou and Bousquet (2007) Léon Bottou and Olivier Bousquet.  The tradeoffs of large scale learning.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 161–168, Vancouver, Canada, 2007.
  * Bottou et al. (2018) Léon Bottou, Frank E. Curtis, and Jorge Nocedal.  Optimization methods for large-scale machine learning.  _SIAM Rev._ , 60(2):223–311, 2018.
  * Bousquet and Elisseeff (2002) Olivier Bousquet and André Elisseeff.  Stability and generalization.  _J. Mach. Learn. Res._ , 2:499–526, 2002.
  * Bousquet et al. (2020) Olivier Bousquet, Yegor Klochkov, and Nikita Zhivotovskiy.  Sharper bounds for uniformly stable algorithms.  In _Proceedings of the Conference on Learning Theory (COLT)_ , pages 610–626, Virtual Event [Graz, Austria], 2020.
  * Chadha et al. (2022) Karan N. Chadha, Gary Cheng, and John C. Duchi.  Accelerated, optimal and parallel: Some results on model-based stochastic optimization.  In _Proceedings of the International Conference on Machine Learning (ICML)_ , pages 2811–2827, Baltimore, MD, 2022.
  * Collobert et al. (2001) Ronan Collobert, Samy Bengio, and Yoshua Bengio.  A parallel mixture of svms for very large scale problems.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 633–640, Vancouver, Canada], 2001.
  * Crammer et al. (2006) Koby Crammer, Ofer Dekel, Joseph Keshet, Shai Shalev-Shwartz, and Yoram Singer.  Online passive-aggressive algorithms.  _J. Mach. Learn. Res._ , 7:551–585, 2006.
  * Davis and Drusvyatskiy (2019) Damek Davis and Dmitriy Drusvyatskiy.  Stochastic model-based minimization of weakly convex functions.  _SIAM J. Optim._ , 29(1):207–239, 2019.
  * Defazio et al. (2014) Aaron Defazio, Francis R. Bach, and Simon Lacoste-Julien.  SAGA: A fast incremental gradient method with support for non-strongly convex composite objectives.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 1646–1654, Montreal, Canada, 2014.
  * Deng and Gao (2021) Qi Deng and Wenzhi Gao.  Minibatch and momentum model-based methods for stochastic weakly convex optimization.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 23115–23127, virtual, 2021.
  * Dieuleveut et al. (2017) Aymeric Dieuleveut, Nicolas Flammarion, and Francis R. Bach.  Harder, better, faster, stronger convergence rates for least-squares regression.  _J. Mach. Learn. Res._ , 18:101:1–101:51, 2017.
  * Drusvyatskiy and Lewis (2018) Dmitriy Drusvyatskiy and Adrian S. Lewis.  Error bounds, quadratic growth, and linear convergence of proximal methods.  _Math. Oper. Res._ , 43(3):919–948, 2018.
  * Duchi and Ruan (2018) John C. Duchi and Feng Ruan.  Stochastic methods for composite and weakly convex optimization problems.  _SIAM J. Optim._ , 28(4):3229–3259, 2018.
  * Duchi et al. (2010) John C. Duchi, Shai Shalev-Shwartz, Yoram Singer, and Ambuj Tewari.  Composite objective mirror descent.  In _Proceedings of the 23rd Conference on Learning Theory (COLT)_ , pages 14–26, Haifa, Israel, 2010.
  * Feldman and Vondrák (2018) Vitaly Feldman and Jan Vondrák.  Generalization bounds for uniformly stable algorithms.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 9770–9780, Montréal, Canada, 2018.
  * Feldman and Vondrák (2019) Vitaly Feldman and Jan Vondrák.  High probability generalization bounds for uniformly stable algorithms with nearly optimal rate.  In _Proceedings of the Conference on Learning Theory (COLT)_ , pages 1270–1279, Phoenix, AZ, 2019.
  * Frostig et al. (2015) Roy Frostig, Rong Ge, Sham M. Kakade, and Aaron Sidford.  Competing with the empirical risk minimizer in a single pass.  In _Proceedings of The 28th Conference on Learning Theory (COLT)_ , pages 728–763, Paris, France, 2015.
  * Ghadimi and Lan (2012) Saeed Ghadimi and Guanghui Lan.  Optimal stochastic approximation algorithms for strongly convex stochastic composite optimization I: A generic algorithmic framework.  _SIAM J. Optim._ , 22(4):1469–1492, 2012.
  * Guyon et al. (2004) Isabelle Guyon, Steve R. Gunn, Asa Ben-Hur, and Gideon Dror.  Result analysis of the NIPS 2003 feature selection challenge.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 545–552, Vancouver, Canada], 2004.
  * Hardt et al. (2016) Moritz Hardt, Ben Recht, and Yoram Singer.  Train faster, generalize better: Stability of stochastic gradient descent.  In _Proceedings of the 33nd International Conference on Machine Learning (ICML)_ , pages 1225–1234, New York City, NY, 2016.
  * Hu et al. (2009) Chonghai Hu, James T. Kwok, and Weike Pan.  Accelerated gradient methods for stochastic optimization and online learning.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 781–789, Vancouver, Canada, 2009.
  * Jaggi et al. (2014) Martin Jaggi, Virginia Smith, Martin Takác, Jonathan Terhorst, Sanjay Krishnan, Thomas Hofmann, and Michael I. Jordan.  Communication-efficient distributed dual coordinate ascent.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 3068–3076, Montreal, Canada, 2014.
  * Johnson and Zhang (2013) Rie Johnson and Tong Zhang.  Accelerating stochastic gradient descent using predictive variance reduction.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 315–323, Lake Tahoe, NV, 2013.
  * Jothimurugesan et al. (2018) Ellango Jothimurugesan, Ashraf Tahmasbi, Phillip B. Gibbons, and Srikanta Tirthapura.  Variance-reduced stochastic gradient descent on streaming data.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 9928–9937, Montréal, Canada, 2018.
  * Karimi et al. (2016) Hamed Karimi, Julie Nutini, and Mark Schmidt.  Linear convergence of gradient and proximal-gradient methods under the polyak-łojasiewicz condition.  In _Proceedings of the European Conference on Machine Learning and Knowledge Discovery in Databases (ECML/PKDD), Part I_ , pages 795–811, Riva del Garda, Italy, 2016.
  * Kasai (2017) Hiroyuki Kasai.  SGDLibrary: A MATLAB library for stochastic optimization algorithms.  _J. Mach. Learn. Res._ , 18:215:1–215:5, 2017.
  * Klochkov and Zhivotovskiy (2021) Yegor Klochkov and Nikita Zhivotovskiy.  Stability and deviation optimal risk bounds with convergence rate $o(1/n)$.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 5065–5076, virtual, 2021.
  * Koltchinskii (2006) Vladimir Koltchinskii.  Local rademacher complexities and oracle inequalities in risk minimization.  _The Annals of Statistics_ , 34(6):2593–2656, 2006.
  * Kulis and Bartlett (2010) Brian Kulis and Peter L. Bartlett.  Implicit online learning.  In _Proceedings of the 27th International Conference on Machine Learning (ICML)_ , pages 575–582, Haifa, Israel, 2010.
  * Kulunchakov and Mairal (2019) Andrei Kulunchakov and Julien Mairal.  A generic acceleration framework for stochastic composite optimization.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 12556–12567, Vancouver, Canada, 2019.
  * Lan (2012) Guanghui Lan.  An optimal method for stochastic composite optimization.  _Math. Program._ , 133(1-2):365–397, 2012.
  * Lee et al. (2017) Jason D. Lee, Qihang Lin, Tengyu Ma, and Tianbao Yang.  Distributed stochastic variance reduced gradient methods by sampling extra data with replacement.  _J. Mach. Learn. Res._ , 18:122:1–122:43, 2017.
  * Lehmann and Casella (2006) Erich L Lehmann and George Casella.  _Theory of point estimation_.  Springer Science & Business Media, 2006.
  * Lei and Ying (2020) Yunwen Lei and Yiming Ying.  Fine-grained analysis of stability and generalization for stochastic gradient descent.  In _Proceedings of the 37th International Conference on Machine Learning (ICML)_ , pages 5809–5819, Virtual Event, 2020.
  * Li et al. (2014) Mu Li, Tong Zhang, Yuqiang Chen, and Alexander J. Smola.  Efficient mini-batch training for stochastic optimization.  In _Proceedings of the 20th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD)_ , pages 661–670, New York, NY, 2014.
  * McDiarmid (1989) C. McDiarmid.  Surveys in combinatorics, 1989: On the method of bounded differences.  1989\.
  * Mei et al. (2018) Song Mei, Yu Bai, and Andrea Montanari.  The landscape of empirical risk for nonconvex losses.  _The Annals of Statistics_ , 46(6A):2747–2774, 2018.
  * Mukherjee et al. (2006) Sayan Mukherjee, Partha Niyogi, Tomaso A. Poggio, and Ryan M. Rifkin.  Learning theory: stability is sufficient for generalization and necessary and sufficient for consistency of empirical risk minimization.  _Adv. Comput. Math._ , 25(1-3):161–193, 2006\.
  * Negahban et al. (2012) Sahand N Negahban, Pradeep Ravikumar, Martin J Wainwright, Bin Yu, et al.  A unified framework for high-dimensional analysis of $m$-estimators with decomposable regularizers.  _Statistical Science_ , 27(4):538–557, 2012.
  * Nemirovski et al. (2009) Arkadi Nemirovski, Anatoli B. Juditsky, Guanghui Lan, and Alexander Shapiro.  Robust stochastic approximation approach to stochastic programming.  _SIAM J. Optim._ , 19(4):1574–1609, 2009.
  * Nemirovskii and Nesterov (1985) Arkaddii S Nemirovskii and Yu E Nesterov.  Optimal methods of smooth convex minimization.  _USSR Computational Mathematics and Mathematical Physics_ , 25(2):21–30, 1985.
  * Patrascu and Necoara (2017) Andrei Patrascu and Ion Necoara.  Nonasymptotic convergence of stochastic proximal point methods for constrained convex optimization.  _J. Mach. Learn. Res._ , 18:198:1–198:42, 2017.
  * Rakhlin et al. (2012) Alexander Rakhlin, Ohad Shamir, and Karthik Sridharan.  Making gradient descent optimal for strongly convex stochastic optimization.  In _Proceedings of the 29th International Conference on Machine Learning (ICML)_ , Edinburgh, Scotland, UK, 2012.
  * Ravikumar et al. (2009) Pradeep Ravikumar, John Lafferty, Han Liu, and Larry Wasserman.  Sparse additive models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 71(5):1009–1030, 2009.
  * Renegar and Grimmer (2022) James Renegar and Benjamin Grimmer.  A simple nearly optimal restart scheme for speeding up first-order methods.  _Found. Comput. Math._ , 22(1):211–256, 2022\.
  * Robbins and Monro (1951) Herbert Robbins and Sutton Monro.  A stochastic approximation method.  _The Annals of Mathematical Statistics_ , pages 400–407, 1951.
  * Robbins and Siegmund (1971) Herbert Robbins and David Siegmund.  A convergence theorem for non negative almost supermartingales and some applications.  In _Optimizing methods in statistics_ , pages 233–257. 1971.
  * Schmidt et al. (2011) Mark Schmidt, Nicolas Le Roux, and Francis R. Bach.  Convergence rates of inexact proximal-gradient methods for convex optimization.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 1458–1466, Granada, Spain, 2011.
  * Shalev-Shwartz et al. (2010) Shai Shalev-Shwartz, Ohad Shamir, Nathan Srebro, and Karthik Sridharan.  Learnability, stability and uniform convergence.  _J. Mach. Learn. Res._ , 11:2635–2670, 2010.
  * Shamir et al. (2014) Ohad Shamir, Nathan Srebro, and Tong Zhang.  Communication-efficient distributed optimization using an approximate newton-type method.  In _Proceedings of the 31th International Conference on Machine Learning (ICML)_ , pages 1000–1008, Beijing, China, 2014.
  * Srebro et al. (2010) Nathan Srebro, Karthik Sridharan, and Ambuj Tewari.  Smoothness, low noise and fast rates.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 2199–2207, Vancouver, Canada, 2010.
  * Toulis and Airoldi (2017) Panos Toulis and Edoardo M Airoldi.  Asymptotic and finite-sample properties of estimators based on stochastic gradients.  _The Annals of Statistics_ , 45(4):1694–1727, 2017.
  * Toulis et al. (2016) Panos Toulis, Dustin Tran, and Edoardo M. Airoldi.  Towards stability and optimality in stochastic gradient descent.  In _Proceedings of the 19th International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pages 1290–1298, Cadiz, Spain, 2016.
  * Tsybakov (2008) Alexandre B Tsybakov.  _Introduction to nonparametric estimation_.  Springer Science & Business Media, 2008.
  * Van de Geer (2008) Sara A Van de Geer.  High-dimensional generalized linear models and the lasso.  _The Annals of Statistics_ , 36(2):614–645, 2008\.
  * Vapnik (1999) Vladimir Vapnik.  An overview of statistical learning theory.  _IEEE Trans. Neural Networks_ , 10(5):988–999, 1999.
  * Wainwright (2009) Martin J. Wainwright.  Sharp thresholds for high-dimensional and noisy sparsity recovery using l1-constrained quadratic programming (lasso).  _IEEE Trans. Inf. Theory_ , 55(5):2183–2202, 2009.
  * Wang et al. (2017a) Jialei Wang, Mladen Kolar, Nathan Srebro, and Tong Zhang.  Efficient distributed learning with sparsity.  In _Proceedings of the 34th International Conference on Machine Learning (ICML)_ , pages 3636–3645, Sydney, Australia, 2017a.
  * Wang et al. (2017b) Jialei Wang, Weiran Wang, and Nathan Srebro.  Memory and communication efficient distributed stochastic optimization with minibatch prox.  In _Proceedings of the 30th Conference on Learning Theory (COLT)_ , pages 1882–1919, Amsterdam, The Netherlands, 2017b.
  * Woodworth and Srebro (2021) Blake E. Woodworth and Nathan Srebro.  An even more optimal stochastic optimization algorithm: Minibatching and interpolation learning.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 7333–7345, virtual, 2021.
  * Xiao and Zhang (2014) Lin Xiao and Tong Zhang.  A proximal stochastic gradient method with progressive variance reduction.  _SIAM J. Optim._ , 24(4):2057–2075, 2014.
  * You et al. (2020) Yang You, Jing Li, Sashank J. Reddi, Jonathan Hseu, Sanjiv Kumar, Srinadh Bhojanapalli, Xiaodan Song, James Demmel, Kurt Keutzer, and Cho-Jui Hsieh.  Large batch optimization for deep learning: Training BERT in 76 minutes.  In _Proceedings of the 8th International Conference on Learning Representations (ICLR)_ , Addis Ababa, Ethiopia, 2020.
  * Yuan and Li (2020) Xiao-Tong Yuan and Ping Li.  On convergence of distributed approximate newton methods: Globalization, sharper bounds and beyond.  _J. Mach. Learn. Res._ , 21:206:1–206:51, 2020.
  * Yuan and Li (2022) Xiao-Tong Yuan and Ping Li.  Stability and risk bounds of iterative hard thresholding.  _IEEE Trans. Inf. Theory_ , 68(10):6663–6681, 2022.
  * Zhang et al. (2017) Lijun Zhang, Tianbao Yang, and Rong Jin.  Empirical risk minimization for stochastic convex optimization: $o(1/n)$-and $o(1/n^{2})$-type of risk bounds.  In _Conference on Learning Theory_ , pages 1954–1979, 2017.
  * Zhang (2003) Tong Zhang.  Leave-one-out bounds for kernel methods.  _Neural Comput._ , 15(6):1397–1437, 2003.
  * Zhang (2004) Tong Zhang.  Solving large scale linear prediction problems using stochastic gradient descent algorithms.  In _Proceedings of the Twenty-first International Conference on Machine Learning (ICML)_ , Banff, Alberta, Canada, 2004.
  * Zhang and Lin (2015) Yuchen Zhang and Xiao Lin.  DiSCO: Distributed optimization for self-concordant empirical loss.  In _Proceedings of the 32nd International Conference on Machine Learning (ICML)_ , pages 362–370, Lille, France, 2015.
  * Zhou et al. (2022) Kaiwen Zhou, Lai Tian, Anthony Man-Cho So, and James Cheng.  Practical schemes for finding near-stationary points of convex finite-sums.  In _Proceedings of the International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pages 3684–3708, Virtual Event, 2022\.
  * Zhou et al. (2019) Pan Zhou, Xiaotong Yuan, Huan Xu, Shuicheng Yan, and Jiashi Feng.  Efficient meta learning via minibatch proximal update.  In _Advances in Neural Information Processing Systems (2019)_ , pages 1532–1542, Vancouver, Canada, 2019.