# Stable Bosonic Topological Edge Modes in the Presence of Many-Body
Interactions

Niclas Heinsdorf  Max-Planck-Institut für Festkörperforschung,
Heisenbergstrasse 1, D-70569 Stuttgart, Germany  Department of Physics and
Astronomy & Stewart Blusson Quantum Matter Institute, University of British
Columbia, Vancouver BC, Canada V6T 1Z4  Darshan G. Joshi  Tata Institute of
Fundamental Research, Hyderabad 500046, India  Hosho Katsura  Department of
Physics, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo,
Tokyo 113-0033, Japan  Institute for Physics of Intelligence, The University
of Tokyo, 7-3-1 Hongo, Tokyo 113-0033, Japan  Trans-scale Quantum Science
Institute, The University of Tokyo, 7-3-1, Hongo, Tokyo 113-0033, Japan
Andreas P. Schnyder  Max-Planck-Institut für Festkörperforschung,
Heisenbergstrasse 1, D-70569 Stuttgart, Germany

###### Abstract

Many magnetic materials are predicted to exhibit bosonic topological edge
modes in their excitation spectra, because of the nontrivial topology of their
magnon, triplon or other quasi-particle band structures. However, there is a
discrepancy between theory prediction and experimental observation, which
suggests some underlying mechanism that intrinsically suppresses the expected
experimental signatures, like the thermal Hall current. Many-body interactions
that are not accounted for in the non-interacting quasi-particle picture are
most often identified as the reason for the absence of the topological edge
modes. Here we report stable bosonic edge modes at the boundaries of a ladder
quantum paramagnet with gapped triplon excitations in the presence of the full
many-body interaction. For the first time, we use tensor network methods to
resolve topological edge modes in the time-dependent spin-spin correlations
and the dynamical structure factor, which is directly accessible
experimentally. We further show that these edge modes have anomalously long
time coherence, discuss the topological phase diagram of the model,
demonstrate the fractionalization of its low-lying excitations, and propose
potential material candidates.

The enormous success of the description of electronic topological phases of
matter quickly inspired research that lead to the generalization of the theory
to bosonic quasi-particles like photons [1, 2, 3, 4, 5, 6], phonons [7, 8, 9,
10], magnons [11, 12, 13, 14, 15, 16, 17, 18, 19] or triplons [20, 21, 22].
Analogous to the Quantum or Spin Hall effect for electrons, the nontrivial
topology of magnetic excitations contributes to unusal responses like the
Thermal Hall or Spin Nernst effect in the form of bosonic edge currents.

Figure 1: (a) Schematic of the spin model on a ladder geometry. On each rung
of the ladder, two spin-1/2 sites are (strongly) coupled through
antiferromagnetic Heisenberg interaction. Along the side rails, the spins
interact (weakly) antiferromagnetically. In this limit the spins form singlet
pairs along the rungs of the ladder. (b) The same model but with anti-
symmetric (DM) and symmetric (pseudo-dipolar) exchange in the $y$-direction
given by $D_{y}$ and $\Gamma_{y}$. These terms introduce a winding and open a
topological gap in the excitation spectrum.

These topologically protected edge spin waves are a key ingredient for many
future spintronic applications, which are a promising alternative to
conventional computing, but with a minimal ecological footprint[23]. There is
a multitude of blueprints that make use of their properties to build next-
generation devices like spin-wave diodes, beam splitters, interferometers and
others [24, 25, 26, 23].

Moreover, systems that are useful for storing and processing quantum
information are highly sought-after. Typically, the quantum information of a
many-body state is destroyed rapidly. Overcoming fast decoherence of quantum
states is one of the main challenges of quantum computing, and has driven
extensive theoretical interest in finding ways to increase their coherence
times using many-body localization [27, 28, 29, 30, 31], prethermalization
[32, 33, 34], quantum many-body scarring [35, 36, 37] or the states’
topological properties[38, 39, 40, 41]. The evolution of topological states
over time - specifically in the presence of defects or nonlinearities[42, 43,
44] \- as well as imaging and detection techniques[45, 46, 47] are thriving
fields of research and crucial for bringing real-world applications on their
way.

Even though a significant amount of material candidates has been proposed to
host nontrivial excitations [48, 49, 11, 50, 51, 14, 25, 52, 53, 54, 20],
their detection remains elusive. While some experiments indeed show signatures
of the predicted bosonic edge modes[55, 56, 57, 58, 59, 60], others have not
been able to produce evidence of topological quasi-particles [61, 62], and for
most of the proposed candidates a clear indication of their presence is yet to
be discovered.

Many reasons that could make the effects of the topological edge modes fail to
materialize have been theorized. Among other things, they include thermal
dampening, magnon-phonon coupling and domain formation, but what is most often
identified as the source of the suppression of thermal Hall effect are many-
body effects. In contrast to their fermionic counterparts, even small many-
body interactions [61, 63, 64] or exchange anisotropies that are typically
present in realistic models [65] seem to substantially affect bosonic
topological transport properties.

In this work we report bosonic topological edge modes that are stable in the
presence of the full many-body interaction using Density Matrix
Renormalization Group (DMRG) and time-evolution methods.

Model and Harmonic Approximation.—We consider a quantum $S=1/2$ Heisenberg
model on a ladder with spin-orbit interaction and an external magnetic field
as shown schematically in Fig. 1. This model, first considered in Ref. [21] is
given by the following Hamiltonian:

$\displaystyle\hat{H}$
$\displaystyle=\hat{H}_{\mathrm{rung}}+\hat{H}_{\mathrm{rail}}+\hat{H}_{\mathrm{SOC}}+\hat{H}_{Z},$
(1) $\displaystyle\hat{H}_{\mathrm{rung}}$
$\displaystyle=J\sum_{i}\hat{\bm{S}}_{li}\cdot\hat{\bm{S}}_{ri},$ (2)
$\displaystyle\hat{H}_{\mathrm{rail}}$
$\displaystyle=K\sum_{i}\left(\hat{\bm{S}}_{li}\cdot\hat{\bm{S}}_{li+1}+\hat{\bm{S}}_{ri}\cdot\hat{\bm{S}}_{ri+1}\right),$
(3) $\displaystyle\hat{H}_{\mathrm{SOC}}$
$\displaystyle=D_{y}\sum_{\alpha=l,r}\sum_{i}\left(\hat{S}^{z}_{\alpha
i}\hat{S}^{x}_{\alpha i+1}-\hat{S}^{x}_{\alpha i}\hat{S}^{z}_{\alpha
i+1}\right),$
$\displaystyle+\Gamma_{y}\sum_{\alpha=l,r}\sum_{i}\left(\hat{S}^{z}_{\alpha
i}\hat{S}^{x}_{\alpha i+1}+\hat{S}^{x}_{\alpha i}\hat{S}^{z}_{\alpha
i+1}\right),$ (4) $\displaystyle\hat{H}_{Z}$
$\displaystyle=h_{y}\sum_{\alpha=l,r}\sum_{i}\left(\hat{S}^{y}_{\alpha
i}+\hat{S}^{y}_{\alpha i}\right),$ (5)

where $\hat{\bm{S}}_{\alpha i}=\left(\hat{S}_{\alpha i}^{x},\hat{S}_{\alpha
i}^{y},\hat{S}_{\alpha i}^{z}\right)^{\intercal}$ are the usual spin operators
with $\alpha=l,r$ corresponding to the left and right rail of the ladder and
$i$ running over all rungs at positions $r_{i}$.

In the limit of strong anti-ferromagnetic rung coupling, this model is the
prototypical example of a gapped quantum paramagnet [66, 67]. In that limit
the ground state of the model is close to a product state of spin singlets

$\displaystyle|\psi_{0}\rangle\sim\bigotimes_{i=1}^{L}(|\uparrow\downarrow\rangle-|\downarrow\uparrow\rangle)/\sqrt{2},$
(6)

with low-lying triplet excitations, the corresponding quasiparticles being
triplons. In the presence of SU(2) symmetry ($D_{y}=\Gamma_{y}=h_{y}=0$), the
triplons are degenerate, but can be split by applying an external magnetic
field due to their different magnetic quantum number. Alternatively, the
degeneracy can be lifted through spin-orbit interactions. Its anti-symmetric
part - similar to the magnetic field term - splits the triplon modes while
retaining U(1) symmetry. The pseudo-dipolar, or $\Gamma$-term, which is the
symmetric part of the spin-orbit interaction breaks U(1)-symmetry.

Figure 2: Real part of the spin-spin correlation function at time $t=50/J$ in
the (a) trivial and (b) the topological case with $K=0.01J$ and field
strengths $h_{y}/D_{y}=0.5$ and $h_{y}/D_{y}=1.5$ respectively on a ladder
with $L=12$ rungs. The topological auto-correlation at this time step has
peaks at the boundaries, which are absent in the trivial case. (c) Time
dependence of the boundary auto-correlation $\mathrm{Re}\ C^{zz}_{11}$ for
$K=0.01J$, $D_{y}=\Gamma_{y}=0.1J$ and no magnetic field on a ladder with
$L=8$ rungs for different boundary conditions. For PBC (no in-gap mode), it
decays rapidly whereas for OBC (in-gap mode) it does not decay. (d) The decay
envelopes of (c) for different values of magnetic field. For finite field
strengths the correlator decays also for OBC, but more slowly than for PBC.

In the paramagnetic phase of the model, an effective low-energy spectrum can
be calculated by first rewriting the model in terms of triplons through bond-
operator theory[68], and then discarding all resulting terms that are not
bilinear, which has been shown to be a controlled approximation [69, 70].
Within this harmonic approximation, it was shown in Ref. [21] that the so-
obtained excitation spectrum can be characterized by a topological winding
number that classifies the topology of the magnet’s bulk and enforces - if
nontrivial - triplon modes at the ends of the ladder with fractional particle
number.

In contrast to the case of topological magnons of ordered magnets [59, 71, 72,
73, 74, 75] \- where there is a gapless Goldstone mode - the triplon
excitations in the quantum paramagnet are gapped. As a consequence of this
bulk gap the triplons, including the topological edge modes, do not have any
low-energy states to decay into. Therefore their lifetime is not reduced
significantly even in the presence of interactions. By applying an external
magnetic field or through spin-anisotropy exchanges, magnon systems can become
gapless too. However, a bulk gap induced by these effects is typically much
smaller than the rung coupling $J$, which is the relevant spin exchange energy
scale that determines the size of the gap in the ladder system.

The spin model presented herein is suitable to describe the magnetism of
potentially many materials and belongs to the well studied class of two-leg
ladder compounds. Examples of spin ladder materials are BiCu2PO6[76],
NaV2O5[77, 78] or multiple cuprate ladders[79, 80, 81, 82]. The recipe for
finding topological edge modes in these systems requires mainly two
ingredients: (i) strong antiferromagnetical rung coupling (stronger than along
the rails), and (ii) spin-orbit interaction. Because the Heisenberg couplings
along the rungs and rails depend on the overlap of the localized electrons,
their ratio can easily be tuned using e.g. strain. The same is true for the
relative spin-orbit interaction strength[83, 84, 85, 86], however, the $D$ and
$\Gamma$ term are additionally restricted by symmetries. These symmetries can
be broken, if not already broken internally, by doping, gate voltages or
structural manipulations like skewing and stretching. The abundance of
experimental handles suggests that a topological phase transition might easily
be switchable, which again expands the potential applications of these systems
in future devices.

Tensor Networks and Dynamical Response.—To relate our model directly to an
experimentally measurable quantity, we compute the dynamic structure factor
(DSF), which is most directly accessed through inelastic neutron
scattering[87, 88] or inverse spin Hall noise spectroscopy (ISHNS)[89].
Because we are investigating the effect of interactions, and since our model
is (quasi) one-dimensional, we use the DRMG algorithm to obtain the exact
many-body ground state of the finite spin ladder[90].

In contrast to electronic systems, the band topology of bosonic quasiparticles
is not a ground state property, so in order to access the low-lying
excitations, we apply time-evolution methods[91, 92, 93, 94, 95, 96] to
compute the odd-part of the time-dependent spin-spin correlation

$\displaystyle
C^{\gamma\gamma^{\prime}}_{ij}(t)=\langle\psi_{0}|\hat{\tilde{S}}^{\gamma}_{i}\hat{U}(t)\hat{\tilde{S}}^{\gamma^{\prime}}_{j}|\psi_{0}\rangle,$
(7)

with
$\hat{\tilde{S}}^{\gamma}_{i}=\hat{S}^{\gamma}_{li}-\hat{S}^{\gamma}_{ri}$ and
the unitary time-evolution operator $\hat{U}(t)$. Already the real-space and
-time correlation function shows signatures of edge modes in the topologically
nontrivial phase region ($h_{y}<|D|$). In Fig. 2 (a) and (b) the real part of
$C^{\gamma\gamma^{\prime}}_{ij}(t)$ for the topological and trivial case at a
fixed time step are plotted. The former shows strong correlations that are
pinned to the edges of the system, whereas they are absent in the latter.

In Fig. 2 (c), we plot the time dependence of $\mathrm{Re}\ C^{xx}_{11}$ at
zero magnetic field for periodic (PBC) and open boundary conditions (OBC). For
PBC there is no topologically protected edge mode at $i=j=1$, and the
excitation decays over time rapidly. In the open system the edge mode has
strongly enhanced time coherence, and we cannot resolve any decay during the
time it takes the excitation to hit the boundary at the other end of the
ladder. Fig. 2 (d) shows the decay envelope of the excitation for increasing
magnetic field for periodic (dashed line) and open (solid line) boundaries.
For finite magnetic fields the boundary correlation starts to decay, and for
stronger fields the difference in lifetime for PBC and OBC - as indicated by
the shaded area - becomes less pronounced.

Figure 3: $S_{k}(\omega)$ for the (a) topological and (b) trivial quantum
paramagnet with $h_{y}/D_{y}=1/2$ and $h_{y}/D_{y}=3/2$ respectively on a
ladder with $L=12$ rungs. The dashed line shows the spectrum of the effective
low-energy models from Ref. [21]. (c) $S_{3\pi/2}(\omega)$ for different
values of $h_{y}/D_{y}$. The first two curves lie in the topological phase
region (see phase diagram at the bottom) with the in-gap modes marked by
dashed lines. At $h_{y}/D_{y}=1$ the gap is closed. The last curve lies in the
trivial region and is gapped, with no mode in-between.

The DSF, which encodes the dynamical susceptibility of the paramagnet is
defined as the Fourier transform of the spin-spin correlations in space and
time

$\displaystyle S^{\gamma\gamma^{\prime}}_{k}(\omega)=\frac{1}{2\pi
L}\int_{-\infty}^{\infty}dt\
e^{i(\omega-\omega_{0})t}\sum_{i,j}e^{i(r_{j}-r_{i})k}C^{\gamma\gamma^{\prime}}_{ij}(t).$
(8)

It is a positive and real quantity and can be calculated using DMRG and time-
evolution methods. Usually, the system size and parameters of the time-
evolution are chosen to minimize finite size effects, and the state is evolved
only for as long as the excitation does not reach the boundary of the system.
In addition, by imposing translational invariance, one of the spatial degrees
of freedom of $C^{\gamma\gamma^{\prime}}_{ij}(t)$ can be fixed [97, 98, 99,
100]. Because we are explicitly studying the excitations in the presence of a
boundary, we require the “full” real-space correlation function, as well as a
“long-enough” time-evolution. $S^{\gamma\gamma^{\prime}}_{k}(\omega)$ obeys
sum rules[99] that we track to diagnose the severity of finite-time and -size
effects, which lead to numerical artifacts (see supplementary material[101]).

We study the field-susceptible part of Eq. (8) which is given by
$S_{k}(\omega)=S^{xx}_{k}(\omega)+S^{zz}_{k}(\omega)$ to which we simply refer
to as DSF from now.

Results.—In Fig. 3 (a) and (b) we plot the DSF of the topological and trivial
paramagnet along with the triplon bands obtained within the harmonic
approximation in Ref. [21]. For small rail and spin-orbit coupling the
harmonic approximation is accurate, because the density of triplons is small
leading to only minor renormalizations from triplon-triplon interactions. In
the topologically nontrivial case (small values of magnetic field) an
additional mode appears between the two bulk bands that we hereafter refer to
as in-gap mode. It is absent in the trivial phase (large values of magnetic
field). Fig. 3 (c) shows the topological phase transition at $k=3\pi/2$ with
the external magnetic field $h_{y}$ as a tuning parameter. At $h_{y}=|D_{y}|$
the gap closes, and further increasing the magnetic field it reopens, but with
no topological in-gap mode. We did finite size scaling with ladders of up to
$L=24$ rungs to confirm that the in-gap mode retains finite spectral weight.

Figure 4: (a) Topological density of states $\rho_{i}^{\text{top}}$ for
$K=0.01J$, $D_{y}=\Gamma_{y}=0.1J$ and different values of $h_{y}$ on a ladder
with $L=8$ rungs. The in-gap mode is localized at the two boundaries of the
system with fractional particle numbers at each termination. For finite
magnetic field the in-gap mode spreads into the bulk and vanishes at the
topological phase transition at $h_{y}/D_{y}=1$. (b) The logarithm of the DSF
for $K=0.01J$, $h_{y}=0.2$ and strong SOC $D_{y}=\Gamma_{y}=1.2J$ on a ladder
with $L=12$ rungs. The quasi-particle peaks are dampened and split. The
boundary mode is flat and marked by text.

To confirm that the in-gap mode is really localized at the boundaries of the
system, we first calculate the local density of states of the lower band

$\displaystyle\rho_{i}=\sum_{k}e^{-ir_{j}k}\int_{0}^{J+\delta}d\omega
S_{k}(\omega),$ (9)

where $\delta$ is a small enough value such that the spectral weight of any
boundary mode, but not that of the upper triplon band is included in the
integral. We then define the “topological density of states” as the difference
of the local densities of states of the nontrivial paramagnet for open and
periodic boundary conditions [21]

$\displaystyle\rho_{i}^{\text{top}}=\rho_{i}^{\text{OBC}}-\rho_{i}^{\text{PBC}}.$
(10)

This quantity isolates the contribution to the spectrum that stems from
introducing a boundary and is shown Fig. 4 (a) for different values of
magnetic field. We see two peaks located at the boundary of the ladder, that
start to spread into the bulk for finite fields. To compute the associated
particle number at one of the edges we integrate $\rho_{i}^{\text{top}}$ for
the system’s left termination

$\displaystyle n_{t}=\int_{1}^{L/2}dr\rho_{i}^{\text{top}},$ (11)

and find that the edge mode has fractional particle number. Our numerical
result $n_{t}\sim 0.43$ deviates slightly from the predicted value of $0.5$
[21] due to spectral leakage caused by the aforementioned finite-size and time
effects. However we expect that in the thermodynamic limit this value would
indeed tend to $0.5$. For larger values of field $n_{t}$ becomes smaller and
vanishes at the phase transition at $h_{y}/D_{y}=1$.

The noninteracting triplon winding number is a topological band property[21],
so the edge modes are naively expected to be stable only as long as there is a
notion of “bands”, or in other words as long as triplons provide a suitable
quasi-particle description of the ladder’s excitation spectrum. We consider
three limits of the model’s parameter space: (i) large magnetic field $h_{y}$,
(ii) large rail coupling $K$ and (iii) large spin-orbit interaction $D_{y}$
and $\Gamma_{y}$. In case (i), the field-polarized phase, the lower triplon
modes condenses at $h_{y}\sim J$ and the system becomes ferromagnetic with
magnons as its low-lying excitations. The band structures derived from linear
spin wave theory are given in the supplemental material[101] and are
topologically trivial. Case (ii) is the limit of two weakly coupled anti-
ferromagnetic Heisenberg chains. The system does not order, even for strong
rail coupling ($K\sim 2J$). The antiferromagnetic order parameter is finite
also for small values of $K$, but approaches zero as the system size $L$ is
increased, as confirmed by finite size scaling [101]. Case (iii) is the most
relevant benchmark of the edge modes’ stability that we can provide with the
available handles in the model, because here the higher-order terms of the
interacting triplon Hamiltonian become important[21]. In this limit there is
significant overlap with the two- and three-triplon continuum that leads to
damping and decay, and the quasi-particle description within the harmonic
approximation breaks down. In contrast to case (i), increasing of $D_{y}$ and
$\Gamma_{y}$ does not lead to condensation of the lower mode, and we did not
find any evidence for a phase transition up to values of
$D_{y}=\Gamma_{y}=2.6J$. In Fig. 4 (b) we plot the logarithm of the DSF for
$D_{y}=\Gamma_{y}=1.2J$. Even though the upper triplon band is strongly damped
and split into many energy levels, we find that the topological edge modes
remain stable.

Conclusion—We have demonstrated that, contrary to what the absence of
experimental evidence for many predicted materials might suggest, topological
bosonic edge modes can be stable in the presence of the full many-body
interaction. Even though correlations have been identified as the reason for
the absence of topological responses in some cases, our work clearly shows
that they do not generically suppress bosonic edge modes, and it prompts the
question of what other mechanism is responsible. A natural extension of this
work is to apply our method to two-dimensional systems wrapped onto a
cylinder. Multi-rail ladders[102, 103], a ferromagnet on a Lieb lattice[104]
or the Shastry-Sutherland model[105] are obvious candidates, but also
topological magnon systems that are potentially more fragile due to their
larger decay phase space.

The absence of signatures of topological edge modes even in systems with good
agreement between theory and experiment for bulk properties suggests that the
responsible decay channel might be a surface effect, i.e., it lies at the
boundary. In a quasi two-dimensional model the edge modes are not completely
localized, but propagate along the boundary of the system. This setup would
further allow the implementation of a defect at the boundary and to
investigate the stability of the topological bosons in its presence.

The typical “workflow” for finding topological edge modes usually consists of
classifying the topology of a translationally invariant model, and then
infering their existence by bulk-boundary correspondence. Even if the symmetry
that protects the band topology is broken weakly, the edge states are still
expected to be present as long as the the perturbation is smaller than the gap
they live in (these approximate symmetries are sometimes called quasi-
symmetries[106, 107]). The chiral symmetry that allows for the definition of a
winding number in the noninteracting case is broken by the rail coupling $K$
[21], and is also clearly broken in the limit of strong SOC, nevertheless our
numerical study shows that the edge modes persist, which invokes the question
of a more suitable bulk classification.

In recent years, considerable effort has been invested in the extension of
topological classifications to interacting systems [108, 109, 110, 111, 112,
113]. These schemes account for the many-body nature of the problem using
either a matrix product state or Green’s function representation, but are not
applicable for the classification of bosonic excitation spectra. Detection of
phase transitions using entanglement have been brought forward, but are - even
though directly related to the dynamical response of a spin system - not
sensitive to topological phase transitions[114]. It is essential to extend
existing or find new methods that are able to capture the topological
properties of e.g. dynamical spin responses, and we hope our work inspires
fruitful research in that direction.

###### Acknowledgements.

We thank A. Nocera and S. M. Winter for helpful discussions. NH acknowledges
financial support from the Max Planck Institute for Solid State Research in
Stuttgart, Germany and the DAAD JSPS summer program 2022. DGJ acknowledges
support from the Department of Atomic Energy, Government of India, under
Project Identification No. RTI 4007. A.P.S. and N.H. acknowledge support by
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – TRR
360 – 492547816. H.K. was supported by JSPS KAKENHI Grant No. JP23H01093 and
MEXT KAKENHI Grant-in-Aid for Transformative Research Areas A “Extreme
Universe” (KAKENHI Grant No. JP21H05191).

## References

  * [1] Onur Hosten and Paul Kwiat.  Observation of the spin hall effect of light via weak measurements.  Science, 319(5864):787–790, 2008.
  * [2] Zheng Wang, Yidong Chong, John D Joannopoulos, and Marin Soljačić.  Observation of unidirectional backscattering-immune topological electromagnetic states.  Nature, 461(7265):772–775, 2009.
  * [3] Philippe Ben-Abdallah.  Photon thermal hall effect.  Physical Review Letters, 116(8):084301, 2016.
  * [4] Frederick Duncan Michael Haldane and Srinivas Raghu.  Possible realization of directional optical waveguides in photonic crystals with broken time-reversal symmetry.  Physical Review Letters, 100(1):013904, 2008.
  * [5] Srinivas Raghu and Frederick Duncan Michael Haldane.  Analogs of quantum-hall-effect edge states in photonic crystals.  Physical Review A, 78(3):033834, 2008.
  * [6] Masaru Onoda, Shuichi Murakami, and Naoto Nagaosa.  Hall effect of light.  Physical review letters, 93(8):083901, 2004.
  * [7] C Strohm, GLJA Rikken, and P Wyder.  Phenomenological evidence for the phonon hall effect.  Physical Review Letters, 95(15):155901, 2005.
  * [8] Yu Kagan and LA Maksimov.  Anomalous hall effect for the phonon heat conductivity in paramagnetic dielectrics.  Physical Review Letters, 100(14):145902, 2008.
  * [9] Kaori Sugii, Masaaki Shimozawa, Daiki Watanabe, Yoshitaka Suzuki, Mario Halim, Motoi Kimata, Yosuke Matsumoto, Satoru Nakatsuji, and Minoru Yamashita.  Thermal hall effect in a phonon-glass ba 3 cusb 2 o 9.  Physical Review Letters, 118(14):145902, 2017.
  * [10] L Sheng, DN Sheng, and CS Ting.  Theory of the phonon hall effect in paramagnetic dielectrics.  Physical Review Letters, 96(15):155901, 2006.
  * [11] Hosho Katsura, Naoto Nagaosa, and Patrick A. Lee.  Theory of the thermal hall effect in quantum magnets.  Phys. Rev. Lett., 104:066403, Feb 2010.
  * [12] Hiroki Kondo, Yutaka Akagi, and Hosho Katsura.  Three-dimensional topological magnon systems.  Physical Review B, 100(14):144401, 2019.
  * [13] Hiroki Kondo and Yutaka Akagi.  Dirac surface states in magnonic analogs of topological crystalline insulators.  Physical Review Letters, 127(17):177201, 2021.
  * [14] Ryo Matsumoto and Shuichi Murakami.  Rotational motion of magnons and the thermal hall effect.  Physical Review B, 84(18):184406, 2011.
  * [15] Satoshi Fujimoto.  Hall effect of spin waves in frustrated magnets.  Physical Review Letters, 103(4):047203, 2009.
  * [16] Ryuichi Shindou, Jun-ichiro Ohe, Ryo Matsumoto, Shuichi Murakami, and Eiji Saitoh.  Chiral spin-wave edge modes in dipolar magnetic thin films.  Physical Review B, 87(17):174402, 2013.
  * [17] Ryuichi Shindou, Ryo Matsumoto, Shuichi Murakami, and Jun-ichiro Ohe.  Topological chiral magnonic edge mode in a magnonic crystal.  Physical Review B, 87(17):174427, 2013.
  * [18] Se Kwon Kim, Héctor Ochoa, Ricardo Zarzuela, and Yaroslav Tserkovnyak.  Realization of the haldane-kane-mele model in a system of localized spins.  Physical Review Letters, 117(22):227201, 2016.
  * [19] Darshan G. Joshi.  Topological excitations in the ferromagnetic kitaev-heisenberg model.  Phys. Rev. B, 98:060405, Aug 2018.
  * [20] Judit Romhányi, Karlo Penc, and Ramachandran Ganesh.  Hall effect of triplons in a dimerized quantum magnet.  Nature communications, 6(1):6805, 2015.
  * [21] Darshan G Joshi and Andreas P Schnyder.  Topological quantum paramagnet in a quantum spin ladder.  Physical Review B, 96(22):220405, 2017.
  * [22] Darshan G. Joshi and Andreas P. Schnyder.  ${\mathbb{z}}_{2}$ topological quantum paramagnet on a honeycomb bilayer.  Phys. Rev. B, 100:020407, Jul 2019.
  * [23] Andrii V Chumak, Vitaliy I Vasyuchka, Alexander A Serga, and Burkard Hillebrands.  Magnon spintronics.  Nature physics, 11(6):453–461, 2015.
  * [24] Z-X Li, Yunshan Cao, and Peng Yan.  Topological insulators and semimetals in classical magnetic systems.  Physics Reports, 915:1–64, 2021.
  * [25] Ryo Matsumoto and Shuichi Murakami.  Theoretical prediction of a rotating magnon wave packet in ferromagnets.  Physical Review Letters, 106(19):197202, 2011.
  * [26] XS Wang, HW Zhang, and XR Wang.  Topological magnonics: A paradigm for spin-wave manipulation and device design.  Physical Review Applied, 9(2):024029, 2018.
  * [27] Igor V Gornyi, Alexander D Mirlin, and Dmitry G Polyakov.  Interacting electrons in disordered wires: Anderson localization and low-t transport.  Physical review letters, 95(20):206603, 2005.
  * [28] Rahul Nandkishore and David A Huse.  Many-body localization and thermalization in quantum statistical mechanics.  Annu. Rev. Condens. Matter Phys., 6(1):15–38, 2015.
  * [29] Ehud Altman and Ronen Vosk.  Universal dynamics and renormalization in many-body-localized systems.  Annu. Rev. Condens. Matter Phys., 6(1):383–409, 2015.
  * [30] Fabien Alet and Nicolas Laflorencie.  Many-body localization: An introduction and selected topics.  Comptes Rendus Physique, 19(6):498–525, 2018.
  * [31] Dmitry A Abanin, Ehud Altman, Immanuel Bloch, and Maksym Serbyn.  Colloquium: Many-body localization, thermalization, and entanglement.  Reviews of Modern Physics, 91(2):021001, 2019.
  * [32] Adam S Jermyn, Roger SK Mong, Jason Alicea, and Paul Fendley.  Stability of zero modes in parafermion chains.  Physical Review B, 90(16):165106, 2014.
  * [33] Jack Kemp, Norman Y Yao, Christopher R Laumann, and Paul Fendley.  Long coherence times for edge spins.  Journal of Statistical Mechanics: Theory and Experiment, 2017(6):063105, 2017.
  * [34] Dominic V Else, Paul Fendley, Jack Kemp, and Chetan Nayak.  Prethermal strong zero modes and topological qubits.  Physical Review X, 7(4):041062, 2017.
  * [35] Maksym Serbyn, Dmitry A Abanin, and Zlatko Papić.  Quantum many-body scars and weak breaking of ergodicity.  Nature Physics, 17(6):675–685, 2021.
  * [36] Anushya Chandran, Thomas Iadecola, Vedika Khemani, and Roderich Moessner.  Quantum many-body scars: A quasiparticle perspective.  Annual Review of Condensed Matter Physics, 14:443–469, 2023.
  * [37] Sanjay Moudgalya, B Andrei Bernevig, and Nicolas Regnault.  Quantum many-body scars and hilbert space fragmentation: a review of exact results.  Reports on Progress in Physics, 85(8):086501, 2022.
  * [38] Daniel E Parker, Romain Vasseur, and Thomas Scaffidi.  Topologically protected long edge coherence times in symmetry-broken phases.  Physical Review Letters, 122(24):240605, 2019.
  * [39] Ivan Amelio and Iacopo Carusotto.  Theory of the coherence of topological lasers.  Physical Review X, 10(4):041060, 2020.
  * [40] Petr Zapletal, Bogdan Galilo, and Andreas Nunnenkamp.  Long-lived elementary excitations and light coherence in topological lasers.  Optica, 7(9):1045–1055, 2020.
  * [41] Yasaman Bahri, Ronen Vosk, Ehud Altman, and Ashvin Vishwanath.  Localization and topology protected quantum coherence at the edge of hot matter.  Nature communications, 6(1):7341, 2015.
  * [42] Rajesh Chaunsali, Haitao Xu, Jinkyu Yang, Panayotis G Kevrekidis, and Georgios Theocharis.  Stability of topological edge states under strong nonlinear effects.  Physical Review B, 103(2):024106, 2021.
  * [43] Daria Smirnova, Daniel Leykam, Yidong Chong, and Yuri Kivshar.  Nonlinear topological photonics.  Applied Physics Reviews, 7(2), 2020.
  * [44] Daniel Leykam and Yi Dong Chong.  Edge solitons in nonlinear-photonic topological insulators.  Physical review letters, 117(14):143901, 2016.
  * [45] Nathan Goldman, Jérôme Beugnon, and Fabrice Gerbier.  Detecting chiral edge states in the hofstadter optical lattice.  Physical review letters, 108(25):255303, 2012.
  * [46] Nathan Goldman, Jean Dalibard, Alexandre Dauphin, Fabrice Gerbier, Maciej Lewenstein, Peter Zoller, and Ian B Spielman.  Direct imaging of topological edge states in cold-atom systems.  Proceedings of the National Academy of Sciences, 110(17):6736–6741, 2013.
  * [47] Yuan Yao, Masahiro Sato, Tetsuya Nakamura, Nobuo Furukawa, and Masaki Oshikawa.  Theory of electron spin resonance in one-dimensional topological insulators with spin-orbit couplings: Detection of edge states.  Physical Review B, 96(20):205424, 2017.
  * [48] António T Costa, Daniel Lourenço R Santos, Nuno MR Peres, and Joaquín Fernández-Rossier.  Topological magnons in cri3 monolayers: an itinerant fermion description.  2D Materials, 7(4):045031, 2020.
  * [49] Shuyi Li and Andriy H Nevidomskyy.  Topological weyl magnons and thermal hall effect in layered honeycomb ferromagnets.  Physical Review B, 104(10):104419, 2021.
  * [50] Kyusung Hwang, Nandini Trivedi, and Mohit Randeria.  Topological magnons with nodal-line and triple-point degeneracies: implications for thermal hall effect in pyrochlore iridates.  Physical Review Letters, 125(4):047203, 2020.
  * [51] Pontus Laurell and Gregory A Fiete.  Magnon thermal hall effect in kagome antiferromagnets with dzyaloshinskii-moriya interactions.  Physical Review B, 98(9):094419, 2018.
  * [52] Alexander Mook, Jürgen Henk, and Ingrid Mertig.  Magnon hall effect and topology in kagome lattices: A theoretical investigation.  Physical Review B, 89(13):134409, 2014.
  * [53] SA Owerre.  Topological honeycomb magnon hall effect: A calculation of thermal hall conductivity of magnetic spin excitations.  Journal of Applied Physics, 120(4), 2016.
  * [54] Jung Hoon Han and Hyunyong Lee.  Spin chirality and hall-like transport phenomena of spin excitations.  Journal of the Physical Society of Japan, 86(1):011007, 2017.
  * [55] Fengfeng Zhu, Lichuan Zhang, Xiao Wang, Flaviano José Dos Santos, Junda Song, Thomas Mueller, Karin Schmalzl, Wolfgang F Schmidt, Alexandre Ivanov, Jitae T Park, et al.  Topological magnon insulators in two-dimensional van der waals ferromagnets crsite3 and crgete3: Toward intrinsic gap-tunability.  Science advances, 7(37):eabi7532, 2021.
  * [56] Heda Zhang, Chunqiang Xu, Caitlin Carnahan, Milos Sretenovic, Nishchay Suri, Di Xiao, and Xianglin Ke.  Anomalous thermal hall effect in an insulating van der waals magnet.  Physical Review Letters, 127(24):247202, 2021.
  * [57] Paul A McClarty, F Krüger, Tatiana Guidi, SF Parker, Keith Refson, AW Parker, Dharmalingam Prabhakaran, and Radu Coldea.  Topological triplon modes and bound states in a shastry–sutherland magnet.  Nature Physics, 13(8):736–741, 2017.
  * [58] Peter Czajka, Tong Gao, Max Hirschberger, Paula Lampen-Kelley, Arnab Banerjee, Nicholas Quirk, David G Mandrus, Stephen E Nagler, and N Phuan Ong.  Planar thermal hall effect of topological bosons in the kitaev magnet $\alpha$-rucl3.  Nature Materials, 22(1):36–41, 2023.
  * [59] Max Hirschberger, Robin Chisnell, Young S Lee, and Nai Phuan Ong.  Thermal hall effect of spin excitations in a kagome magnet.  Physical Review Letters, 115(10):106603, 2015.
  * [60] Y Onose, T Ideue, H Katsura, Y Shiomi, N Nagaosa, and Y Tokura.  Observation of the magnon hall effect.  Science, 329(5989):297–299, 2010.
  * [61] S Suetsugu, T Yokoi, K Totsuka, T Ono, I Tanaka, S Kasahara, Y Kasahara, Z Chengchao, H Kageyama, and Y Matsuda.  Intrinsic suppression of the topological thermal hall effect in an exactly solvable quantum magnet.  Physical Review B, 105(2):024415, 2022.
  * [62] Luke Pritchard Cairns, J-Ph Reid, Robin Perry, Dharmalingam Prabhakaran, and Andrew Huxley.  Thermal hall effect of topological triplons in srcu2 (bo3) 2.  In Proceedings of the International Conference on Strongly Correlated Electron Systems (SCES2019), page 011089, 2020.
  * [63] ME Zhitomirsky and AL Chernyshev.  Colloquium: Spontaneous magnon decays.  Reviews of Modern Physics, 85(1):219, 2013.
  * [64] Jonas Habel, Alexander Mook, Josef Willsher, and Johannes Knolle.  Breakdown of chiral edge modes in topological magnon insulators.  arXiv preprint arXiv:2308.03168, 2023.
  * [65] Andreas Thomasen, Karlo Penc, Nic Shannon, and Judit Romhányi.  Fragility of z 2 topological invariant characterizing triplet excitations in a bilayer kagome magnet.  Physical Review B, 104(10):104412, 2021.
  * [66] D Schmidiger, Pierre Bouillot, S Mühlbauer, S Gvasaliya, Corinna Kollath, Thierry Giamarchi, and A Zheludev.  Spectral and thermodynamic properties of a strong-leg quantum spin ladder.  Physical review letters, 108(16):167201, 2012.
  * [67] Elbio Dagotto and TM Rice.  Surprises on the way from one-to two-dimensional quantum magnets: The ladder materials.  Science, 271(5249):618–623, 1996.
  * [68] Subir Sachdev and RN Bhatt.  Bond-operator representation of quantum spins: Mean-field theory of frustrated quantum heisenberg antiferromagnets.  Physical Review B, 41(13):9323, 1990.
  * [69] Darshan G. Joshi, Kris Coester, Kai P. Schmidt, and Matthias Vojta.  Nonlinear bond-operator theory and $1/d$ expansion for coupled-dimer magnets. i. paramagnetic phase.  Phys. Rev. B, 91:094404, Mar 2015.
  * [70] Darshan G. Joshi and Matthias Vojta.  Nonlinear bond-operator theory and $1/d$ expansion for coupled-dimer magnets. ii. antiferromagnetic phase and quantum phase transition.  Phys. Rev. B, 91:094405, Mar 2015.
  * [71] Song Bao, Jinghui Wang, Wei Wang, Zhengwei Cai, Shichao Li, Zhen Ma, Di Wang, Kejing Ran, Zhao-Yang Dong, DL Abernathy, et al.  Discovery of coexisting dirac and triply degenerate magnons in a three-dimensional antiferromagnet.  Nature communications, 9(1):2591, 2018.
  * [72] M Elliot, Paul A McClarty, D Prabhakaran, RD Johnson, HC Walker, P Manuel, and R Coldea.  Order-by-disorder from bond-dependent exchange and intensity signature of nodal quasiparticles in a honeycomb cobaltate.  Nature Communications, 12(1):3936, 2021.
  * [73] Weiliang Yao, Chenyuan Li, Lichen Wang, Shangjie Xue, Yang Dan, Kazuki Iida, Kazuya Kamazawa, Kangkang Li, Chen Fang, and Yuan Li.  Topological spin excitations in a three-dimensional antiferromagnet.  Nature Physics, 14(10):1011–1015, 2018.
  * [74] Bo Yuan, Ilia Khait, Guo-Jiun Shu, FC Chou, MB Stone, JP Clancy, Arun Paramekanti, and Young-June Kim.  Dirac magnons in a honeycomb lattice quantum xy magnet cotio 3.  Physical Review X, 10(1):011062, 2020.
  * [75] L-C Zhang, YA Onykiienko, PM Buhl, YV Tymoshenko, P Čermák, A Schneidewind, JR Stewart, A Henschel, M Schmidt, S Blügel, et al.  Magnonic weyl states in cu 2 oseo 3.  Physical review research, 2(1):013063, 2020.
  * [76] Olivier Mentré, El Mostafa Ketatni, Marie Colmont, Marielle Huve, Francis Abraham, and Vaclav Petricek.  Structural features of the modulated bicu2 (p1-x v x) o6 solid solution; 4-d treatment of x= 0.87 compound and magnetic spin-gap to gapless transition in new cu2+ two-leg ladder systems.  Journal of the American Chemical Society, 128(33):10857–10867, 2006\.
  * [77] Frédéric Mila, Patrice Millet, and Jacques Bonvoisin.  Exchange integrals of vanadates as revealed by magnetic-susceptibility measurements of nav 2 o 5.  Physical Review B, 54(17):11925, 1996.
  * [78] Masahiko Isobe and Yutaka Ueda.  Magnetic susceptibility of quasi-one-dimensional compound $\alpha$’-n a v 2 o 5–possible spin-peierls compound with high critical temperature of 34 k–.  Journal of the Physical Society of Japan, 65(5):1178–1181, 1996\.
  * [79] Z Hiroi and M Takano.  Absence of superconductivity in the doped antiferromagnetic spin-ladder compound (la, sr) cuo2. 5.  Nature, 377(6544):41–43, 1995.
  * [80] Yoshimitsu Kohama, Shuang Wang, Atsuko Uchida, Krunoslav Prsa, Sergei Zvyagin, Yuri Skourski, Ross D McDonald, Luis Balicas, Henrik M Ronnow, Christian Rüegg, et al.  Anisotropic cascade of field-induced phase transitions in the frustrated spin-ladder system bicu 2 po 6.  Physical Review Letters, 109(16):167204, 2012.
  * [81] V Kiryukhin, YJ Kim, KJ Thomas, FC Chou, RW Erwin, Q Huang, MA Kastner, and RJ Birgeneau.  Magnetic properties of the s= 1 2 quasi-one-dimensional antiferromagnet cacu 2 o 3.  Physical Review B, 63(14):144418, 2001.
  * [82] Masaki Azuma, Z Hiroi, M Takano, K Ishida, and Y Kitaoka.  Observation of a spin gap in sr cu 2 o 3 comprising spin-$1/2$ quasi-1d two-leg ladders.  Physical Review Letters, 73(25):3463, 1994.
  * [83] Ya-ping Wang, Chang-wen Zhang, Wei-xiao Ji, Run-wu Zhang, Ping Li, Pei-ji Wang, Miao-juan Ren, Xin-lian Chen, and Min Yuan.  Tunable quantum spin hall effect via strain in two-dimensional arsenene monolayer.  Journal of Physics D: Applied Physics, 49(5):055305, 2016.
  * [84] Chi-Ho Cheung, Huei-Ru Fuh, Ming-Chien Hsu, Yeu-Chung Lin, and Ching-Ray Chang.  Spin orbit coupling gap and indirect gap in strain-tuned topological insulator-antimonene.  Nanoscale Research Letters, 11:1–8, 2016.
  * [85] Lars Winterfeld, Luis A Agapito, Jin Li, Nicholas Kioussis, Peter Blaha, and Yong P Chen.  Strain-induced topological insulator phase transition in hgse.  Physical Review B, 87(7):075143, 2013.
  * [86] Tsai-Jung Liu, Maximilian A Springer, Niclas Heinsdorf, Agnieszka Kuc, Roser Valentí, and Thomas Heine.  Semimetallic square-octagon two-dimensional polymer with high mobility.  Physical Review B, 104(20):205419, 2021.
  * [87] K. Sturm.  Dynamic structure factor: An introduction.  Zeitschrift für Naturforschung A, 48(1-2):233–242, 1993.
  * [88] Léon Van Hove.  Correlations in space and time and born approximation scattering in systems of interacting particles.  Physical Review, 95(1):249, 1954.
  * [89] Darshan G Joshi, Andreas P Schnyder, and So Takei.  Detecting end states of topological quantum paramagnets via spin hall noise spectroscopy.  Physical Review B, 98(6):064401, 2018.
  * [90] Steven R White.  Density matrix formulation for quantum renormalization groups.  Physical Review Letters, 69(19):2863, 1992.
  * [91] Guifré Vidal.  Efficient simulation of one-dimensional quantum many-body systems.  Physical Review Letters, 93(4):040502, 2004.
  * [92] Andrew John Daley, Corinna Kollath, Ulrich Schollwöck, and Guifré Vidal.  Time-dependent density-matrix renormalization-group using adaptive effective hilbert spaces.  Journal of Statistical Mechanics: Theory and Experiment, 2004(04):P04005, 2004.
  * [93] Steven R White and Adrian E Feiguin.  Real-time evolution using the density matrix renormalization group.  Physical review letters, 93(7):076401, 2004.
  * [94] Jutho Haegeman, J Ignacio Cirac, Tobias J Osborne, Iztok Pižorn, Henri Verschelde, and Frank Verstraete.  Time-dependent variational principle for quantum lattices.  Physical Review Letters, 107(7):070601, 2011.
  * [95] Jutho Haegeman, Christian Lubich, Ivan Oseledets, Bart Vandereycken, and Frank Verstraete.  Unifying time evolution and optimization with matrix product states.  Physical Review B, 94(16):165116, 2016.
  * [96] Michael P Zaletel, Roger SK Mong, Christoph Karrasch, Joel E Moore, and Frank Pollmann.  Time-evolving a matrix product state with long-ranged interactions.  Physical Review B, 91(16):165112, 2015.
  * [97] Nicholas E Sherman, Maxime Dupont, and Joel E Moore.  Spectral function of the j 1- j 2 heisenberg model on the triangular lattice.  Physical Review B, 107(16):165146, 2023.
  * [98] Pierre Bouillot, Corinna Kollath, Andreas M Läuchli, Mikhail Zvonarev, Benedikt Thielemann, Christian Rüegg, Edmond Orignac, Roberta Citro, Martin Klanjšek, Claude Berthier, et al.  Statics and dynamics of weakly coupled antiferromagnetic spin-1 2 ladders in a magnetic field.  Physical Review B, 83(5):054407, 2011.
  * [99] Pontus Laurell, Allen Scheie, D Alan Tennant, Satoshi Okamoto, Gonzalo Alvarez, and Elbio Dagotto.  Magnetic excitations, nonclassicality, and quantum wake spin dynamics in the hubbard chain.  Physical Review B, 106(8):085110, 2022.
  * [100] Matthias Gohlke, Roderich Moessner, and Frank Pollmann.  Dynamical and topological properties of the kitaev model in a [111] magnetic field.  Physical Review B, 98(1):014418, 2018.
  * [101] supplemental.
  * [102] Rhine Samajdar, Mathias S Scheurer, Shubhayu Chatterjee, Haoyu Guo, Cenke Xu, and Subir Sachdev.  Enhanced thermal hall effect in the square-lattice néel state.  Nature Physics, 15(12):1290–1294, 2019.
  * [103] Masataka Kawano and Chisa Hotta.  Thermal hall effect and topological edge states in a square-lattice antiferromagnet.  Physical Review B, 99(5):054422, 2019.
  * [104] Xiaodong Cao, Kai Chen, and Dahai He.  Magnon hall effect on the lieb lattice.  Journal of Physics: Condensed Matter, 27(16):166003, 2015.
  * [105] B Sriram Shastry and Bill Sutherland.  Exact ground state of a quantum mechanical antiferromagnet.  Physica B+ C, 108(1-3):1069–1070, 1981.
  * [106] Chunyu Guo, Lunhui Hu, Carsten Putzke, Jonas Diaz, Xiangwei Huang, Kaustuv Manna, Feng-Ren Fan, Chandra Shekhar, Yan Sun, Claudia Felser, et al.  Quasi-symmetry-protected topology in a semi-metal.  Nature physics, 18(7):813–818, 2022.
  * [107] Chao Shang, Owen Ganter, Niclas Heinsdorf, and Stephen M Winter.  Irf 4: From tetrahedral compass model to topological semimetal.  Physical Review B, 107(12):125111, 2023.
  * [108] Victor Gurarie.  Single-particle green’s functions and interacting topological insulators.  Physical Review B, 83(8):085426, 2011.
  * [109] Zhong Wang and Binghai Yan.  Topological hamiltonian as an exact tool for topological invariants.  Journal of Physics: Condensed Matter, 25(15):155601, 2013.
  * [110] Zhong Wang and Shou-Cheng Zhang.  Simplified topological invariants for interacting insulators.  Physical Review X, 2(3):031008, 2012.
  * [111] Mikel Iraola, Niclas Heinsdorf, Apoorv Tiwari, Dominik Lessnich, Thomas Mertz, Francesco Ferrari, Mark H Fischer, Stephen M Winter, Frank Pollmann, Titus Neupert, et al.  Towards a topological quantum chemistry description of correlated systems: The case of the hubbard diamond chain.  Physical Review B, 104(19):195125, 2021.
  * [112] Anton Kapustin and Lev Spodyneiko.  Higher-dimensional generalizations of berry curvature.  Physical Review B, 101(23):235130, 2020.
  * [113] Ken Shiozaki, Niclas Heinsdorf, and Shuhei Ohyama.  Higher berry curvature from matrix product states.  arXiv preprint arXiv:2305.08109, 2023.
  * [114] Philipp Hauke, Markus Heyl, Luca Tagliacozzo, and Peter Zoller.  Measuring multipartite entanglement through dynamic susceptibilities.  Nature Physics, 12(8):778–782, 2016.
  * [115] Johannes Hauschild and Frank Pollmann.  Efficient numerical simulations with tensor networks: Tensor network python (tenpy).  SciPost Physics Lecture Notes, page 005, 2018.
  * [116] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant.  Array programming with NumPy.  Nature, 585(7825):357–362, September 2020.
  * [117] Niclas Heinsdorf.  topwave python package.  to be published, see Supplementary.

## Supplementary Material

##  I DMRG and Time-Evolution

We use DMRG and the $\hat{W}^{\textup{II}}$ time-evolution method as
implemented in the tenpy package [115]. The sampling rate (or time step
$\Delta t$) is chosen such that the fastest frequency of the signal is
resolved, which is the ground state frequency $\omega_{0}$. The resolution of
the Fourier transformed output along the frequency axis - but also the
severity of finite size effects due to reflections at the boundaries - is then
determined by the total evolved time $T=N_{\Delta t}\Delta t$, where
$N_{\Delta t}$ is the number of time steps.

To reduce ringing effects the correlation functions are convoluted with a
hamming window along the time axis, before being Fourier transformed using
Fast Fourier Transform (FFT) as implemented in the numpy package [116]. To
obtain the Fourier transform in real-space, all matrix elements $C_{ij}$ are
collected with respect to their relative distance along the ladder $d=i-j$.

Because the DSF is a response function, it is required to be real and
positive. In principle, this can be used to choose a suitable Fourier
transform, the output of which will always meet said requirements. In this
work, we choose an agnostic, pre-implemented FFT to be able to judge the
severity of finite-size and -time effects by keeping track of the deviations
from these properties as well as the violation of sum rules.

The algorithm that is used to compute the DSF is summarized as follows:

  1. 1.

Compute ground state $|\psi_{0}\rangle$ using (finite) DMRG.

  2. 2.

Compute
$|\psi^{\gamma^{\prime}}_{j}(t=0)\rangle=\hat{\tilde{S}}^{\gamma^{\prime}}_{j}|\psi_{0}\rangle$.

  3. 3.

Time-evolve one time step using $\hat{W}^{\textup{II}}$:

$\displaystyle\hat{U}(\Delta
t)|\psi^{\gamma^{\prime}}_{j}(t=0)\rangle=|\psi^{\gamma^{\prime}}_{j}(t+\Delta
t)\rangle$

  4. 4.

Compute $|\psi^{\gamma\gamma^{\prime}}_{ij}(t+\Delta
t)\rangle=\hat{\tilde{S}}^{\gamma}_{j}|\psi^{\gamma}_{i}(t+\Delta t)\rangle$
for all $i$.

  5. 5.

Compute overlap $C^{\gamma\gamma^{\prime}}_{ij}(t+\Delta
t)=\langle\psi_{0}|\psi^{\gamma\gamma^{\prime}}_{ij}(t+\Delta t)\rangle$

  6. 6.

Repeat steps 3 \- 5 $N_{\Delta t}=T/\Delta t$ times.

  7. 7.

Repeat steps 2 \- 6 for all $j$.

  8. 8.

Fourier transform $C^{\gamma\gamma^{\prime}}_{ij}(t)$ in time and sites.

##  II Case (i): Strong Field and Linear Spin Wave Theory

For strong external magnetic fields $h_{y}/J\sim 1$, there is a transition to
the field polarized phase as can be seen in Fig. S1. This ordered state has
magnons as its low-lying excitations. We perform Linear Spin Wave Theory
(LSWT) as implemented in the topwave package[117], and plot the magnon band
structures for a variety of different parameters in Fig. S2. In the parameter
regime relevant to our model, no band crossing can be enforced.

##  III Case (ii): Strong Rail Coupling

In the limit of strong rail coupling $K$, the ladder model resembles two
weakly coupled anti-ferromagnetic Heisenberg chains. We define the
antiferromagnetic order parameter (for the left rail) as

$\displaystyle\Delta^{\text{AFM}}=\frac{1}{\left(L/2\right)^{2}}\sum_{\gamma=x,y,z}\sum_{i,j}^{L}(-1)^{i+j}\langle\hat{S}^{\gamma}_{li}\hat{S}^{\gamma}_{lj}\rangle/3.$
(12)

The order parameter for the right rail ($\alpha=r)$ is equivalent. For finite
system size and finite SOC, the order parameter is nonzero, but approaches
zero for larger system sizes as shown in Fig. S3. Hence, the system does not
order, as expected for one dimension.

##  IV Case (iii): Strong Spin-Orbit Coupling

Even though SOC gaps out the two triplon modes, the lower mode does not
condense even for large values of $D_{y}$ and $\Gamma_{y}$. In Fig. S4 the
ground state energy and its derivatives are given. There is no sign of a phase
transition. Fig. S5 shows the logarithm of the DSF for strong SOC
($\Gamma_{y}=D_{y}=2.6J$) and no magnetic field. The two triplon modes are
split into many different energy levels, which are completely flat. The lower
modes do not condense, but are instead “repelled” while the upper modes rise
in energy and are strongly damped. The topological boundary mode is stable.

Figure S1: The ground state energy, its derivative and the magnetization along
the external field axis as a function of field strength $h_{y}/J$. At
$h_{y}/J\sim 1$ there is a phase transition from the paramagnetic to the
field-polarized phase. Figure S2: Magnon spectra of the field polarized ladder
for $h_{y}/J=2.5$ obtained by LSWT as implemented in the topwave package[117]
for different values of rail coupling $K$ and SOC. The model has two magnetic
sites per unit cell and hence two magnon bands. Two bands of the same model
are given in the same color. There are no band crossings and the bands are
topologically trivial. Figure S3: (a) Antiferrromagnetic order parameter for
$h_{y}=0$, $D_{y}=\Gamma_{y}=0.1J$ on a ladder with $L=24$ rungs. The order
parameter is small, but finite even for $K=0$. (b) Finite size scaling of the
order parameter at $K=0$. The order parameter approaches zero in the
thermodynamic limit. Figure S4: Ground state energy (a) and its first (b) and
second (c) derivative with respect to SOC $D_{y}$ (in units of $J$) with
$K=0.01J$, $h_{y}=0$ and $\Gamma_{y}=D_{y}$ on a ladder with $L=12$ rungs. The
curves are smooth and there is no evidence of a phase transition. Figure S5:
The logarithm of the DSF for $K=0.01J$, $h_{y}=0$ and $\Gamma_{y}=D_{y}=2.6J$
on a ladder with $L=12$ rungs. There are many flat energy levels with the ones
at higher energies being damped strongly. The boundary mode remains present
and is marked.


$, which together with Assumption <ref> gives
\sup \limits_{\theta \in \Theta}\big\langle  F(\hat{\theta}), \hat{\theta}-\theta\big\rangle  \leqslant \frac{1}{\beta}
\big\langle  \hat{\theta}-\theta_{k}, \hat{\theta}-\theta\big\rangle \leqslant \frac{2D}{\beta} \left\|\theta_{k}-\hat{\theta}\right\|.

[Zhang et al., 2021]
If Assumption $\ref{assumption0}$ hold, then for any $\theta=\left(\theta_{i}, \theta_{-i}\right)$ and  $\hat{\theta}^{i}=\left(\theta_{i}^{\prime}, \theta_{-i}\right)
\in \Theta$, we have
$J_{i}^{(\theta_{i}^{\prime}, \theta_{-i})}(t) -J_{i}^{(\theta_{i}, \theta_{-i})}(t)
=\frac{1}{1-\gamma} \int_{\mathcal{S}\times \mathcal{A}_{i}} d_{\rho_{t}}^{\hat{\theta}^{i}}(s^{\prime}) \pi_{\theta_{i}^{\prime}}\left(a_{i} \mid s^{\prime}\right) \\
\times \Big(\int_{\mathcal{A}_{-i}} \pi_{\theta_{-i}}\left(a_{-i} \mid s^{\prime}\right) A_{i}^{\theta}\left(s^{\prime}, a,t+1\right) \mathrm{d} a_{-i}\Big)\mathrm{d} a_{i} \mathrm{d} s^{\prime}$, $\forall \ i \in \mathcal{N}$.

Define proximal mapping $P_{\Theta_{i}}(x-y): \mathbb{R}^{d_{i}} \times \mathbb{R}^{d_{i}} \rightarrow \Theta_{i}$ as

P_{\Theta_{i}}(x-y)=\underset{z \in \Theta_{i}}{\operatorname{argmin}}\left\{2\langle y, z\rangle+\left\|x-z\right\|^{2}\right\}, \quad \forall \ x ,\\\ y \in \mathbb{R}^{d_{i}}.

(i) For any $x$, $y$, $z \in \mathbb{R}^{d_{i}}$,

\frac{\left\|P_{\Theta_{i}}(x-y)- z\right\|^{2}}{2}
\leqslant \frac{\left\|x- z\right\|^{2}}{2}+\left\langle y, z-P_{\Theta_{i}}(x-y)\right\rangle-\frac{\left\|x-P_{\Theta_{i}}(x-y))\right\|^{2}}{2};

(ii) For any $x_{1}$, $z \in \mathbb{R}^{d_{i}}$ and $y_{1}$, $y_{2} \in \mathbb{R}^{d_{i}}$,

\begin{equation}\nonumber
\begin{aligned}
& \frac{\left\|P_{\Theta_{i}}(x_{1}-y_{2})- z\right\|^{2}}{2}\\
\leqslant &  \frac{\left\|x_{1}- z\right\|^{2}}{2}-\left\langle y_{2}, P_{\Theta_{i}}(x_{1}-y_{1})-z)\right\rangle+
\left\langle y_{2}-y_{1}, P_{\Theta_{i}}(x_{1}-y_{1})-P_{\Theta_{i}}(x_{1}-y_{2})\right\rangle\\
\leqslant &\frac{\left\|x_{1}- z\right\|^{2}}{2}-\left\langle y_{2}, P_{\Theta_{i}}(x_{1}-y_{1})-z)\right\rangle+\frac{\left\|y_{2}-y_{1}\right\|^{2}}{2}-\frac{\left\|P_{\Theta_{i}}(x_{1}-y_{1})- x_{1}\right\|^{2}}{2};
\end{aligned}
\end{equation}

(iii) Non-expansion  property:

\begin{equation}\nonumber
\left\|P_{\Theta_{i}}(x_{1}-y_{1})-P_{\Theta_{i}}(x_{1}-y_{2})\right\| \leqslant \left\|y_{1}-y_{2}\right\|.
\end{equation}

[Liu et al., 2021] If $G$ is a $\sqrt{2L^{2}+\frac{2}{\beta^{2}}}$-Lipschitz continous mapping and strong-
ly monotone and let $w_{*}$ be the solution of SVI$(G, K)$, then for any $\widehat{z} \in K$, by constructing  $\overline{z}=\operatorname{P}_{K}(\widehat{z}-\widetilde{\eta} G(\widehat{z}))$, where $\widetilde{\eta}=\frac{1}{\sqrt{2} \sqrt{2L^{2}+\frac{2}{\beta^{2}}}}$, we have $\sup \limits_{z \in K} G(\overline{z})^{\top}(\overline{z}-z) \leqslant 2D \sqrt{2L^{2}+\frac{2}{\beta^{2}}}(2+\sqrt{2})
\left\|\widehat{z}-w_{*}\right\|$, where $\sup \limits_{z, z^{\prime} \in K}\left\|z-z^{\prime}\right\| \leqslant 2D$ and $L$ and $\beta$ are constants.

[Pinelis, 1994]  (Concentration inequality) If $X_{1}, X_{2}, ..., X_{N} \in \mathbb{R}^{d}$
denote a vector-valued martingale difference sequence satisfying $\big\|X_{n}\big\| \leqslant V$ and $\mathbb{E}\big[X_{n} \big| X_{1}, ..., X_{n-1}\big]\\=\mathbf{0}$, $\forall \ n \in \{1,\ldots,N\}$, then for any $\delta \in(0,1]$, we have

P\left\{\left\|\sum_{n=1}^{N} X_{n}\right\|^{2}>2 \log (2 / \delta) V^{2} N\right\} \leqslant \delta.

[Chung, 1954] If $\{u_{k}, k\geqslant 1\}$ is a real sequence satisfying

u_{k+1} \leqslant \left(1-\frac{c}{k^{s}}\right) u_{k}+\frac{d}{k^{t}}, \quad 0<s<1, \ s<t, \ c>0, \ d>0,

u_{k} \leqslant \frac{d}{c} \frac{1}{k^{t-s}}+\phi(k),

where $\phi(k)=o\left(\frac{1}{k^{t-s}}\right)$.

§ APPENDIX B PROOFS OF THEOREM <REF>, THEOREM <REF> AND LEMMA <REF>

Proof of Theorem <ref>:
For any countable set $\widetilde{\Theta_{i}} \subseteq \Theta_{i}$, by the property of conditional expectation, we have

\begin{align}
&\mathbb{E}\left[\sup _{\theta_{i} \in \widetilde{\Theta_{i}}} J_{i}^{(\theta_{i}, \theta_{-i,\tau_K})}(t)-J_{i}^{(\theta_{i,\tau_K}, \theta_{-i,\tau_K})}(t)\biggm|\theta_{k},k=1,...,K\right] \notag\\
=&\mathbb{E}\left\{\mathbb{E}\left[\sup _{\theta_{i} \in \widetilde{\Theta_{i}}}J_{i}^{(\theta_{i}, \theta_{-i,\tau_K})}(t)-J_{i}^{(\theta_{i,\tau_K}, \theta_{-i,\tau_K})}(t)\biggm|\tau_K,\theta_{k},k=1,...,K\right]\biggm|\theta_{k},k=1,...,K\right\}\notag\\
=&\sum_{k=1}^{K}\mathbb{E}\left[\sup _{\theta_{i} \in \widetilde{\Theta_{i}}}J_{i}^{(\theta_{i}, \theta_{-i,\tau_K})}(t)-J_{i}^{(\theta_{i,\tau_K}, \theta_{-i,\tau_K})}(t)\biggm|\tau_K=k, \theta_{k},k=1,...,K\right] \notag\\
=&\sum_{k=1}^{K}\mathbb{E}\left[\sup _{\theta_{i} \in \widetilde{\Theta_{i}}}J_{i}^{(\theta_{i}, \theta_{-i,k})}(t)-J_{i}^{(\theta_{i,k}, \theta_{-i,k})}(t)\biggm|\tau_K=k, \theta_{k},k=1,...,K\right]\notag\\
\label{conditionalexpectation}
\end{align}

Noting that $\tau_{K}$ is independent of $\theta_{k}$, it follows that $P\bigg\{\tau_K=k\biggm|\theta_{k},k=1,...,K\bigg\}=P\bigg\{\tau_K=k\bigg\}$.
It's easy to see that

\begin{align}
&\mathbb{E}\left[\sup _{\theta_{i} \in \widetilde{\Theta_{i}}}J_{i}^{(\theta_{i}, \theta_{-i,k})}(t)-J_{i}^{(\theta_{i,k}, \theta_{-i,k})}(t)\biggm|\tau_K=k, \theta_{k},k=1,...,K\right]\notag\\
=&\mathbb{E}\left[\sup _{\theta_{i} \in \widetilde{\Theta_{i}}}J_{i}^{(\theta_{i}, \theta_{-i,k})}(t)-J_{i}^{(\theta_{i,k}, \theta_{-i,k})}(t)\biggm|\theta_{k},k=1,...,K\right].\label{independence}
\end{align}

From the definition of total reward function, we have
$$\mathbb{E}\left[\sup _{\theta_{i} \in \widetilde{\Theta_{i}}}J_{i}^{(\theta_{i}, \theta_{-i,k})}(t)-J_{i}^{(\theta_{i,k}, \theta_{-i,k})}(t)\biggm|\theta_{k},k=1,...,K\right]
=\sup _{\theta_{i} \in \widetilde{\Theta_{i}}}J_{i}^{(\theta_{i}, \theta_{-i,k})}(t)-J_{i}^{(\theta_{i,k}, \theta_{-i,k})}(t)

This together with (<ref>) gives

\begin{equation}\nonumber
\begin{aligned}
&\mathbb{E}\left[\sup _{\theta_{i} \in \widetilde{\Theta_{i}}}J_{i}^{(\theta_{i}, \theta_{-i,k})}(t)-J_{i}^{(\theta_{i,k}, \theta_{-i,k})}(t)\biggm|\tau_K=k, \theta_{k},k=1,...,K\right]\\
=&\sup _{\theta_{i} \in \widetilde{\Theta_{i}}}J_{i}^{(\theta_{i}, \theta_{-i,k})}(t)-J_{i}^{(\theta_{i,k}, \theta_{-i,k})}(t).
\end{aligned}
\end{equation}

Combining the above inequality with (<ref>) gives

\begin{align}
&\mathbb{E}\left[\sup_{\theta_{i}\in \widetilde{\Theta_{i}}}J_{i}^{(\theta_{i}, \theta_{-i,\tau_K})}(t)-J_{i}^{(\theta_{i,\tau_K}, \theta_{-i,\tau_K})}(t)\biggm|\theta_{k},k=1,...,K\right] \notag\\
=&\sum_{k=1}^{K}\left(\sup_{\theta_{i} \in \widetilde{\Theta_{i}}}J_{i}^{(\theta_{i}, \theta_{-i,k})}(t)-J_{i}^{(\theta_{i,k}, \theta_{-i,k})}(t)\right)P\bigg\{\tau_K=k\bigg\} \notag\\
=&\frac{\sum_{k=1}^{K}\gamma_{k}\left(\sup \limits_{\theta_{i} \in \widetilde{\Theta_{i}}}J_{i}^{(\theta_{i}, \theta_{-i,k})}(t)-J_{i}^{(\theta_{i,k}, \theta_{-i,k})}(t)\right)}{\sum_{k=1}^{K}\gamma_{k}}.\notag
\end{align}

Noting that $\left(\left(\pi_{\theta_{i,k}}\right)_{i=1}^{N}\right)_{k=1}^{\infty}$ is a $\gamma_{k}$-weighted $\epsilon_{k}$-Nash equilibrium of the game $\Gamma$, it follows that

\begin{equation}\nonumber
\begin{aligned}
&\sup_{i \in \mathcal{N}}
\mathbb{E}\left[\sup_{\theta_{i}\in \widetilde{\Theta_{i}}}J_{i}^{(\theta_{i}, \theta_{-i,\tau_K})}(t)-J_{i}^{(\theta_{i,\tau_K}, \theta_{-i,\tau_K})}(t)\biggm|\theta_{k},k=1,...,K\right]\\
=& \sup_{i \in \mathcal{N}}\frac{\sum_{k=1}^{K}\gamma_{k}\left(\sup \limits _{\theta_{i} \in \widetilde{\Theta_{i}}}J_{i}^{(\theta_{i}, \theta_{-i,k})}(t)-J_{i}^{(\theta_{i,k}, \theta_{-i,k})}(t)\right)}{\sum_{k=1}^{K}\gamma_{k}}
\leqslant \epsilon_{K},
\end{aligned}
\end{equation}

then the theorem holds.

Proof of Theorem <ref>:
By Assumption <ref> and  Lemma 3.1 in [Hartman and Stampacchia, 1966], we know that SVI$(F,\Theta)$ has a solution. By Assumption <ref> and Assumption <ref>, Lemma <ref> implies the equivalence between the Nash equilibrium problem and SVI$(F,\Theta)$. Therefore, the game has a Nash equilibrium.

Proof of Lemma <ref>:
For agent $i$, by Assumptions
<ref>-<ref>, Lemma $3.2$ in [Zhang et al., 2020] and Proposition <ref>, we have

\begin{equation}\nonumber
\begin{aligned}
\nabla_{\theta_{i}} V_{i}^{\theta}(s,t)=\sum_{l=t}^{\infty} \sum_{\tau=0}^{\infty} \gamma^{l+\tau-t} \int & r_{i}\left(s(l+\tau), a(l+\tau)\right)  \nabla_{\theta_{i}} \log \pi_{\theta_{i}} \left(a_{i}(l) \mid s(l)\right)\\
& \rho_{\theta, t: l+\tau}  \mathrm{d} s_{t+1} ... \mathrm{d} s_{l+\tau} \mathrm{d} a_{t+1}  ... \mathrm{d} a_{l+\tau},
\end{aligned}
\end{equation}

where $
\rho_{\theta, t:l+\tau}=\left[\prod_{u=t}^{l+\tau-1} \rho\left(s(u+1) \mid s(u), a(u)\right)\right] \left[\prod_{u=t}^{l+\tau} \pi_{\theta}\left(a(u) \mid s(u)\right)\right]
is the probability density function of the trajectory $\big(s(t), a(t),..., s(l+\tau), a(l+\tau)
\big)$ and the integral is  over all trajectories. Then, it follows that

\begin{align}\label{differencenabla}
&\left\|\nabla_{{\theta}_{i}} V_{i}^{\theta^{1}}(s,t)-\nabla_{{\theta}_{i}} V_{i}^{\theta^{2}}(s,t)\right\| \notag\\
=&\Bigg\| \sum_{l=t}^{\infty} \sum_{\tau=0}^{\infty} \gamma^{l+\tau-t} \Bigg(\int \left\{\nabla_{\theta_{i}} \log \pi_{\theta^{1}_{i}}\left(a_{i}(l) \mid s(l)\right)-\nabla_{\theta_{i}} \log \pi_{\theta^{2}_{i}}\left(a_{i}(l)  \mid s(l)\right)\right\}\notag\\
&\times r_{i}\left(s(l+\tau), a(l+\tau)\right)\rho_{\theta^{1}, t: l+\tau} +\int r_{i}\left(s(l+\tau), a(l+\tau)\right)  \nabla_{\theta_{i}} \log \pi_{\theta^{2}_{i}}\left(a_{i}(l) \mid s(l)\right)\notag\\
& \times \left(\rho_{\theta^{1}, t: l+\tau}-\rho_{\theta^{2}, t: l+\tau}\right) \Bigg) \mathrm{d} s_{t+1} ... \mathrm{d} s_{l+\tau} \mathrm{d} a_{t+1}  ... \mathrm{d} a_{l+\tau} \Bigg\| \notag\\
\leqslant &\sum_{l=t}^{\infty} \sum_{\tau=t}^{\infty} \gamma^{l+\tau-t} \Big(\int \left\|\nabla_{\theta_{i}} \log \pi_{\theta^{1}_{i}}\left(a_{i}(l)  \mid s(l)\right)-\nabla_{\theta_{i}} \log \pi_{\theta^{2}_{i}}\left(a_{i}(l) \mid s(l)\right)\right\| \notag\\
&\times \left|r_{i}\left(s(l+\tau), a(l+\tau)\right)\right| \rho_{\theta^{1}, t: l+\tau} \mathrm{d} s_{t+1} ... \mathrm{d} s_{l+\tau} \mathrm{d} a_{t+1}  ... \mathrm{d} a_{l+\tau}
+\int\left|r_{i}\left(s(l+\tau), a(l+\tau)\right)\right| \notag\\
&\times \left\|\nabla_{\theta_{i}} \log \pi_{\theta^{2}_{i}}\left(a_{i}(l) \mid s(l)\right)\right\|
\left|\rho_{\theta^{1}, t: l+\tau}-\rho_{\theta^{2}, t: l+\tau}\right|\Big) \mathrm{d} s_{t+1} ... \mathrm{d} s_{l+\tau} \mathrm{d} a_{t+1}  ... \mathrm{d} a_{l+\tau}.
\end{align}

\begin{align}
I_{1}=&\int\left|r_{i}\left(s(l+\tau), a(l+\tau)\right)\right| \left\|\nabla_{\theta_{i}} \log \pi_{\theta^{1}_{i}}\left(a_{i}(l)  \mid s(l)\right)-\nabla_{\theta_{i}} \log \pi_{\theta^{2}_{i}}\left(a_{i}(l) \mid s(l)\right)\right\| \notag\\
&\times \rho_{\theta^{1}, t: l+\tau} \mathrm{d} s_{t+1} ... \mathrm{d} s_{l+\tau} \mathrm{d} a_{t+1}  ... \mathrm{d} a_{l+\tau},\label{I11}\\
I_{2}=&\int\left|r_{i}\left(s(l+\tau), a(l+\tau)\right)\right| \left\|\nabla_{\theta_{i}} \log \pi_{\theta^{2}_{i}}\left(a_{i}(l) \mid s(l)\right)\right\|\notag\\
& \times
\left|\rho_{\theta^{1}, t: l+\tau}-\rho_{\theta^{2}, t: l+\tau}\right| \mathrm{d} s_{t+1} ... \mathrm{d} s_{l+\tau} \mathrm{d} a_{t+1}  ... \mathrm{d} a_{l+\tau}.\label{I22}
\end{align}

For the term $I_{1}$, by  Assumption <ref>, we have

\begin{equation}\label{I1}
I_{1} \leqslant U_{R}  L_{\Theta}  \left\|\theta^{1}_{i}-\theta^{2}_{i}\right\|.
\end{equation}

For the term $I_{2}$, denote $\mathcal{U}_{l+\tau}=\{u: u=t, ..., l+\tau\}$. From the definitions of  $\rho_{\theta^{1}, t: l+\tau}$ and $\rho_{\theta^{2}, t:l+\tau}$, we have

\begin{align}\label{density}
\rho_{\theta^{1}, t: l+\tau}-\rho_{\theta^{2}, t: l+\tau}=&\left[\prod_{u=t}^{l+\tau-1} \rho\left(s(u+1) \mid s(u), a(u)\right)\right] \notag\\
& \times\Bigg[\prod_{u \in \mathcal{U}_{l+\tau}} \pi_{\theta^{1}}\left(a(u) \mid s(u)\right)-
\prod_{u \in \mathcal{U}_{l+\tau}} \pi_{\theta^{2}}\left(a(u) \mid s(u)\right)\Bigg].
\end{align}

From Assumption <ref> and Taylor expansion of $\prod_{u \in \mathcal{U}_{l+\tau}} \pi_{\theta}\left(a(u) \mid s(u)\right)$ near $\theta=\theta^{1}$, there exists some $\widetilde{\theta}=\lambda \theta^{1}+(1-\lambda) \theta^{2} $ and $\lambda \in[0,1]$, such that

\begin{align}\nonumber
&\left|\prod_{u \in \mathcal{U}_{l+\tau}} \pi_{\theta^{1}}\left(a(u) \mid s(u)\right)-\prod_{u \in \mathcal{U}_{l+\tau}} \pi_{\theta^{2}}\left(a(u) \mid s(u)\right)\right| \notag\\
=&\left|\left(\theta^{1}-\theta^{2}\right)^{\top}\left[\sum_{m \in \mathcal{U}_{l+\tau}} \nabla \pi_{\widetilde{\theta}}\left(a(m) \mid s(m)\right) \prod_{u \in \mathcal{U}_{l+\tau}, u \neq m} \pi_{\widetilde{\theta}}\left(a(u) \mid s(u)\right)\right]\right| \notag\\
= &\left|\left(\theta^{1}-\theta^{2}\right)^{\top}\left[\sum_{m \in \mathcal{U}_{l+\tau}} \nabla \log \pi_{\widetilde{\theta}}\left(a(m) \mid s(m)\right) \prod_{u \in \mathcal{U}_{l+\tau}} \pi_{\widetilde{\theta}}\left(a(u) \mid s(u)\right)\right]\right| \notag\\
\leqslant &\left\|\theta^{1}-\theta^{2}\right\| \sum_{m \in \mathcal{U}_{l+\tau}}\left\|\nabla \log \pi_{\widetilde{\theta}}\left(a(m) \mid s(m)\right)\right\| \prod_{u \in \mathcal{U}_{l+\tau}} \pi_{\tilde{\theta}}\left(a(u) \mid s(u)\right) \notag\\
\leqslant &\left\|\theta^{1}-\theta^{2}\right\| (l+\tau-t+1)  B_{\Theta} \sqrt{N}  \prod_{u \in \mathcal{U}_{l+\tau}} \pi_{\widetilde{\theta}}\left(a(u) \mid s(u)\right)\notag ,
\end{align}

where $\nabla \log \pi_{\widetilde{\theta}}\left(a(u) \mid s(u)\right)=(\nabla_{\theta_{i}} \log \pi_{\theta_{i}}(a_{i}(u) \mid s(u))_{i=1}^{N}$.
From (<ref>), (<ref>), the above inequality and Assumption <ref>, it follows that

\begin{align*}
I_{2}  \leqslant &\left\|\theta^{1}-\theta^{2}\right\|  U_{R} \sqrt{N}  B_{\Theta}^{2} \int\left[\prod_{u=t}^{l+\tau-1} \rho\left(s(u+1) \mid s(u), a(u)\right)\right] (l+\tau-t+1)\notag\\
&\times \prod_{u \in \mathcal{U}_{l+\tau}} \pi_{\widetilde{\theta}}\left(a(u) \mid s(u)\right)  \mathrm{d} s_{t+1} ... \mathrm{d} s_{l+\tau} \mathrm{d} a_{t+1}  ... \mathrm{d} a_{l+\tau} \\
=&\left\|\theta^{1}-\theta^{2}\right\|  U_{R} \sqrt{N} B_{\Theta}^{2} (l+\tau-t+1).
\end{align*}

By (<ref>), (<ref>) and the above inequality, we have

\begin{align*}
\left\|\nabla_{{\theta}_{i}} V_{i}^{\theta^{1}}(s,t)-\nabla_{{\theta}_{i}} V_{i}^{\theta^{2}}(s,t)\right\| &\leqslant \sum_{l=t}^{\infty} \sum_{\tau=0}^{\infty}\gamma^{l+\tau-t}\left(I_{1}+I_{2}\right)\\
\leqslant & \sum_{l=t}^{\infty} \sum_{\tau=0}^{\infty} \gamma^{l+\tau-t}  \Bigg(U_{R} L_{\Theta}  \left\|\theta^{1}_{i}-\theta^{2}_{i}\right\|\\
&+(l+\tau-t+1) U_{R}  \sqrt{N}  B_{\Theta}^{2} \left\|\theta^{1}-\theta^{2}\right\|\Bigg)\\
\leqslant & \frac{1}{(1-\gamma)^{2}}  U_{R}  L_{\Theta} \left\|\theta^{1}_{i}-\theta^{2}_{i}\right\|+\frac{1+\gamma}{(1-\gamma)^{3}} U_{R}  \sqrt{N} B_{\Theta}^{2} \left\|\theta^{1}-\theta^{2}\right\| \ .
\end{align*}

This together with Assumption <ref> and Proposition <ref> gives

\begin{align*}
\left\|\nabla_{{\theta}_{i}} J_{i}^{\theta^{1}}(t)-\nabla_{{\theta}_{i}} J_{i}^{\theta^{2}}(t)\right\|
=&\left\|\int_{s \in \mathcal{S}}\left(\nabla_{{\theta}_{i}} V_{i}^{\theta^{1}}(s,t)-\nabla_{{\theta}_{i}} V_{i}^{\theta^{2}}(s,t)\right)\rho_{t}(s) \mathrm{d}s\right\|\\
\leqslant & \int_{s \in \mathcal{S}}
\left\|\nabla_{{\theta}_{i}} V_{i}^{\theta^{1}}(s,t)-\nabla_{{\theta}_{i}} V_{i}^{\theta^{2}}(s,t)\right\| \rho_{t}(s) \mathrm{d}s\\
\leqslant &\frac{1}{(1-\gamma)^{2}}  U_{R}  L_{\Theta} \left\|\theta^{1}_{i}-\theta^{2}_{i}\right\|+\frac{1+\gamma}{(1-\gamma)^{3}} U_{R}  \sqrt{N} B_{\Theta}^{2} \left\|\theta^{1}-\theta^{2}\right\| \ .
\end{align*}

Hence, we have

\begin{align*}
&\sum_{i=1}^{N} \left\|\nabla_{{\theta}_{i}} J_{i}^{\theta^{1}}(t)-\nabla_{{\theta}_{i}} J_{i}^{\theta^{2}}(t)\right\|^{2}\\
\leqslant & 2\sum_{i=1}^{N}\left(\frac{ U_{R} \cdot L_{\Theta}}{(1-\gamma)^{2}} \right)^{2} \left\|\theta^{1}_{i}-\theta^{2}_{i}\right\|^{2} +2 \left(\frac{1+\gamma}{(1-\gamma)^{3}}  U_{R} \sqrt{N}  B_{\Theta}^{2}\right)^{2} \left\|\theta^{1}-\theta^{2}\right\|^{2}\\
\leqslant &2\left(\frac{ U_{R}  L_{\Theta}}{(1-\gamma)^{2}} \right)^{2} \left\|\theta^{1}-\theta^{2}\right\|^{2} +2 \left(\frac{1+\gamma}{(1-\gamma)^{3}}  U_{R} \sqrt{N} B_{\Theta}^{2}\right)^{2} \left\|\theta^{1}-\theta^{2}\right\|^{2}\\
=&\left[2\left(\frac{ U_{R}  L_{\Theta}}{(1-\gamma)^{2}} \right)^{2}+2 \left(\frac{1+\gamma}{(1-\gamma)^{3}}  U_{R}  \sqrt{N} B_{\Theta}^{2}\right)^{2}\right]
\left\|\theta^{1}-\theta^{2}\right\|^{2},
\end{align*}

that is,

$$\left\|F(\theta^{1})-F(\theta^{2})\right\| \leqslant  L\left\|\theta^{1}-\theta^{2}\right\|.$$

[Agarwal et al., 2020]
A. Agarwal, S. M. Kakade, J. D. Lee, and G. Mahajan.
Optimality and approximation with policy gradient methods in Markov
decision processes.
In Proceedings of the 33rd Conference on Learning Theory,
pages 64–66, Graz, Austria, 2020.

[Baxter and Bartlett, 2001]
J. Baxter and P. L. Bartlett.
Infinite-horizon policy-gradient estimation.
Journal of Artificial Intelligence Research, 15:0
319–350, 2001.

[Chen et al., 2021]
T. Y. Chen, K. Zhang, G. B. Giannakis, and T. Başar.
Communication efficient policy gradient methods for distributed
reinforcement learning.
IEEE Transactions on Control of Network Systems, forthcoming,

[Chung, 1954]
K. L. Chung.
On a stochastic approximation method.
Annals of Mathematical Statistics, 250 (3):0
463–483, 1954.

[Daskalakis et al., 2006]
C. Daskalakis, A. Mehta, and C. Papadimitriou.
A note on approximate Nash equilibria.
In Proceedings of the 2nd International Workshop on Internet
and Network Economics, pages 297–306, Patras, Greece, 2006.

[Daskalakis et al., 2020]
C. Daskalakis, D. J. Foster, and N. Golowich.
Independent policy gradient methods for competitive reinforcement
Advances in Neural Information Processing Systems,
33:0 5527–5540, 2020.

[Facchinei and Pang, 2003]
F. Facchinei and J. S. Pang.
Finite-Dimensional Variational Inequalities and Complementarity
Problems, Vol. I.
Springer-Verlag, New York, USA, 2003.

[Filar and Vrieze, 2012]
J. A. Filar and K. Vrieze.
Competitive Markov Decision Processes.
Springer Science and Business Media, Berlin, German, 2012.

[Filar et al., 1991]
J. A. Filar, T. A. Schultz, F. Thuijsman, and O. J. Vrieze.
Nonlinear programming and stationary equilibria in stochastic games.
Mathematical Programming, 500 (1):0 227–237,

[Greenwald et al., 2003]
A. Greenwald, K. Hall, and R. Serrano.
Correlated Q-learning.
In Proceedings of the 20th International Conference on
International Conference on Machine Learning, volume 20, pages 242–249,
Washington, USA, 2003.

[Hartman and Stampacchia, 1966]
P. Hartman and G. Stampacchia.
On some nonlinear elliptic differential-functional equations.
Acta Mathematica, 115:0 271–310, 1966.

[Hoffman and Karp, 1966]
A. Hoffman and R. Karp.
On nonterminating stochastic games.
Management Science, 120 (5):0 359–370, 1966.

[Hsieh et al., 2019]
Y. G. Hsieh, J. Malick F. Iutzeler, and P. Mertikopoulos.
On the convergence of single-call stochastic extra-gradient methods.
In Proceedings of the 33rd International Conference on Neural
Information Processing Systems, pages 6938–6948, Virtual, 2019.

[Hu and Wellman, 1998]
J. Hu and M. P. Wellman.
Multiagent reinforcement learning: theoretical framework and an
In Proceedings of the 15th International Conference on Machine
Learning, volume 98, pages 242–250, Wisconson, USA, 1998.

[Hu and Wellman, 2003]
J. Hu and M. P. Wellman.
Nash Q-learning for general-sum stochastic games.
Journal of Machine Learning Research, 4:0 1039–1069,

[Kinderlehrer and Stampacchia, 1980]
D. Kinderlehrer and G. Stampacchia.
An Introduction to Variational Inequalities and Their
Academic Press, New York, USA, 1980.

[Konda and Tsitsiklis, 1999]
V. Konda and J. Tsitsiklis.
Actor-Critic Algorithms.
Advances in Neural Information Processing Systems, pages
1008–1014, 1999.

[Koshal et al., 2010]
J. Koshal, A. Nedić, and U. V. Shanbhag.
Single timescale regularized stochastic approximation schemes for
monotone Nash games under uncertainty.
In Proceedings of the 49th IEEE Conference on Decision and
Control, page 231–236, Georgia, USA, 2010.

[Koshal et al., 2013]
J. Koshal, A. Nedić, and U. V. Shanbhag.
Regularized iterative stochastic approximation methods for stochastic
variational inequality problems.
IEEE Transactions on Automatic Control, 580
(3):0 594–609, 2013.

[Leonardos et al., 2021]
S. Leonardos, W. Overman, I. Panageas, and G. Piliouras.
Global convergence of multi-agent policy gradient in Markov
potential games.
arXiv preprint arXiv:2106.01969, 2021.

[Littman, 1994]
M. L. Littman.
Markov games as a framework for multi-agent reinforcement learning.
In Proceedings of the 11th International Conference on Machine
Learning, pages 157–163, New Jersey, USA, 1994.

[Littman, 2001]
M. L. Littman.
Friend-or-Foe Q-learning in general-sum games.
In Proceedings of the 18th International Conference on Machine
Learning, volume 1, pages 322–328, Massachusetts, USA, 2001.

[Liu et al., 2021]
M. Liu, H. Rafique, Q. Lin, and T. Yang.
First-order convergence theory for weakly-convex-weakly-concave
min-max problems.
Journal of Machine Learning Research, 220
(169):0 1–34, 2021.

[Liu et al., 2020]
Y. Liu, K. Zhang, T. Başar, and W. Yin.
An improved analysis of (variance-reduced) policy gradient and
natural policy gradient methods.
Advances in Neural Information Processing Systems,
33:0 7624–7636, 2020.

[Lu et al., 2021]
S. Lu, K. Zhang, T. Chen, T. Başar, and L. Horesh.
Decentralized policy gradient descent ascent for safe multi-agent
reinforcement learning.
In Proceedings of the 35th AAAI Conference on Artificial
Intelligence, volume 35, pages 8767–8775, Virtual, 2021.

[Mao and Başar, 2022]
W. C. Mao and T. Başar.
Provably efficient reinforcement learning in decentralized generalsum
Markov games.
Dynamic Games and Applications, pages 1–22, 2022.

[Minty, 1962]
G. J. Minty.
Monotone (nonlinear) operators in Hilbert space.
Duke Mathematical Journal, 290 (3):0 341–346,

[Nash, 1951]
J. F. Nash.
Non-cooperative games.
Annals of Mathematics, 540 (2):0 286–295,

[Perolat et al., 2018]
J. Perolat, B. Piot, and O. Pietquin.
Actor-Critic fictitious play in simultaneous move multistage games.
In Proceedings of the 21st International Conference on
Artificial Intelligence and Statistics, pages 919–928, Playa Blanca,
Lanzarote, 2018.

[Pinelis, 1994]
I. Pinelis.
Optimum bounds for the distributions of martingales in Banach
The Annals of Probability, 220 (4):0
1679–1706, 1994.

[Prasad et al., 2015]
H. L. Prasad, L. A. Prashanth, and S. Bhatnagar.
Two-timescale algorithms for learning Nash equilibria in
general-sum stochastic games.
In Proceedings of the 2015 International Conference on
Autonomous Agents and Multiagent Systems, pages 1371–1379, Istanbul,
Turkey, 2015.

[Shapley, 1953]
L. S. Shapley.
Stochastic games.
Proceedings of the National Academy of Sciences of the United
States of America, pages 1095–1100, 1953.

[Song et al., 2020]
C. Song, Z. Zhou, Y. Jiang Y. Zhou, and Y. Ma.
Optimistic dual extrapolation for coherent non-monotone variational
Advances in Neural Information Processing Systems,
33:0 14303–14314, 2020.

[Sutton et al., 1999]
R. S. Sutton, S. P. Singh D. A. McAllester, and Y. Mansour.
Policy gradient methods for reinforcement learning with function
Advances in Neural Information Processing Systems, pages
1057–1063, 1999.

[Watkins and Dayan, 1992]
C. J. Watkins and P. Dayan.
Machine Learning, 80 (3):0 279–292, 1992.

[Wei et al., 2021]
C. Y. Wei, C. W. Lee, M. X. Zhang, and H. P. Luo.
Last-iterate convergence of decentralized optimistic gradient
descent/ascent in infinite-horizon competitive Markov games.
In Proceedings of the 34th Conference on Learning Theory,
volume 134, pages 4259–4299, Colorado, USA, 2021.

[Yang et al., 2018]
Z. Yang, M. Hong K. Zhang, and T. Başar.
A finite sample analysis of the Actor-Critic algorithm.
In Proceedings of the 57th IEEE Conference on Decision and
Control, Florida, USA, 2018.

[Zhang et al., 2020]
K. Zhang, A. Koppel, H. Zhu, and T. Başar.
Global convergence of policy gradient methods to (almost) locally
optimal policies.
SIAM Journal on Control and Optimization, 580
(6):0 3586–3612, 2020.

[Zhang et al., 2021]
R. Zhang, Z. Ren, and N. Li.
Gradient play in stochastic games: stationary points, convergence,
and sample complexity.
arXiv preprint arXiv: 2106.00198, 2021.

[Zhao et al., 2022]
Y. Zhao, Y. Tian, J. D. Lee, and S. S. Du.
Provably efficient policy gradient methods for two-player zero-sum
Markov games.
In Proceedings of the 25th International Conference on
Artificial Intelligence and Statistics, volume 151, pages 2736–2761,
Virtual, 2022.