#  Observation of magnon bound states in the long-range, anisotropic
Heisenberg model

Florian Kranzl  The first two coauthors contributed equally. Institute for
Quantum Optics and Quantum Information, Austrian Academy of Sciences,
Technikerstraße 21a, 6020 Innsbruck, Austria  Institut für Experimentalphysik,
Universität Innsbruck, Technikerstraße 25, 6020 Innsbruck, Austria  Stefan
Birnkammer  The first two coauthors contributed equally. Department of
Physics, Technical University of Munich, 85748 Garching, Germany  Munich
Center for Quantum Science and Technology (MCQST), Schellingstr. 4, 80799
München, Germany  Manoj K. Joshi  Institute for Quantum Optics and Quantum
Information, Austrian Academy of Sciences, Technikerstraße 21a, 6020
Innsbruck, Austria  Alvise Bastianello  Department of Physics, Technical
University of Munich, 85748 Garching, Germany  Munich Center for Quantum
Science and Technology (MCQST), Schellingstr. 4, 80799 München, Germany
Rainer Blatt  Institute for Quantum Optics and Quantum Information, Austrian
Academy of Sciences, Technikerstraße 21a, 6020 Innsbruck, Austria  Institut
für Experimentalphysik, Universität Innsbruck, Technikerstraße 25, 6020
Innsbruck, Austria  Michael Knap  Department of Physics, Technical University
of Munich, 85748 Garching, Germany  Munich Center for Quantum Science and
Technology (MCQST), Schellingstr. 4, 80799 München, Germany  Christian F. Roos
Institute for Quantum Optics and Quantum Information, Austrian Academy of
Sciences, Technikerstraße 21a, 6020 Innsbruck, Austria  Institut für
Experimentalphysik, Universität Innsbruck, Technikerstraße 25, 6020 Innsbruck,
Austria

###### Abstract

Over the recent years coherent, time-periodic modulation has been established
as a versatile tool for realizing novel Hamiltonians. Using this approach,
known as Floquet engineering, we experimentally realize a long-ranged,
anisotropic Heisenberg model with tunable interactions in a trapped ion
quantum simulator. We demonstrate that the spectrum of the model contains not
only single magnon excitations but also composite magnon bound states. For the
long-range interactions with the experimentally realized power-law exponent,
the group velocity of magnons is unbounded. Nonetheless, for sufficiently
strong interactions we observe bound states of these unconventional magnons
which possess a non-diverging group velocity. By measuring the configurational
mutual information between two disjoint intervals, we demonstrate the
implications of the bound state formation on the entanglement dynamics of the
system. Our observations provide key insights into the peculiar role of
composite excitations in the non-equilibrium dynamics of quantum many-body
systems.

Recent conceptual and technical progress has enabled the controlled study of
quantum magnetism with quantum simulators and quantum processors. For example,
the dynamics of magnons, which are itinerant flipped spin excitations that
propagate through the system by the superexchange mechanism [1, 2, 3], has
been probed in quantum gas microscopes [4]. Already in the early 1930s it has
been realized that in addition to isolated magnons, bound states of magnons
exist as low-energy excitations of a ferromagnet [1]. For short-range
superexchange interactions, the remarkable light-cone dynamics of magnon bound
states has been investigated [5] and their robustness has recently been
demonstrated as well [6]. Long-range spin interactions significantly modify
this picture. When considering long-range spin flip-flop interactions, as
realized by trapped ions in a strong transverse field [7, 8, 9], magnons can
acquire an unbounded group velocity [10, 11]. Hence, they are rather exotic
quasiparticles. This raises the question of whether interactions between
magnons, which enable the bound-state formation in the short-range model, can
also stabilize bound states in a long-range spin chain. Recent theoretical
work suggested that this is indeed the case [12]. It remains an open challenge
to experimentally realize and control models which allow for the observation
of such bound states.

In this work, we investigate bound-state dynamics in the long-range,
anisotropic Heisenberg model realized in a linear string of ions, where two
electronic states of the ion represent a spin-1/2. Spin-spin interactions can
be engineered by illuminating the ions with a bichromatic laser beam. This
coupling gives rise to an effective long-range Ising interaction
$H_{\text{XX}}=\sum_{i<j}J_{ij}\sigma_{i}^{x}\sigma_{j}^{x}$ [13] with spin-
spin couplings that decay approximately algebraically with distance,
$J_{ij}=J/|i-j|^{\alpha}$. A key challenge is to enrich the set of available
spin-spin interactions, starting from the long-range Ising interactions. By
applying a strong transverse field, obtained by introducing either a third
laser frequency component or by shifting the frequency of the bichromatic
beam, spin-changing processes are penalized. This effectively realizes the
long-range XY model describing flip-flop dynamics of spins [7, 8]. Enhancing
the symmetry of the spin interactions from a discrete Ising to a continuous XY
symmetry, enabled for example the spectroscopy of quasiparticle excitations
[11], the variational simulation of target Hamiltonians with related
symmetries [14], the study of hydrodynamic spin transport [15], and the
observation of spontaneous breaking of a continuous symmetry [16]. However, in
the long-range XY model direct interactions between magnons, which are
necessary for the bound state formation [12], are absent. In order to realize
these interactions between magnons, we resort to Floquet engineering [17], and
implement basis rotations such that the averaged Hamiltonian in the high
frequency limit yields the long-range, anisotropic Heisenberg model

$\hat{H}=\frac{1}{3}\sum_{i<j}J_{ij}(\hat{\sigma}_{i}^{x}\hat{\sigma}^{x}_{j}+\hat{\sigma}_{i}^{y}\hat{\sigma}^{y}_{j}+\Delta\hat{\sigma}_{i}^{z}\hat{\sigma}^{z}_{j})\,,$
(1)

where $\hat{\sigma}^{x,y,z}$ are the Pauli matrices and $\Delta$ determines
the Ising interactions between the magnons; see Fig. 1 (a). Several
implementations of long-range Heisenberg models have been proposed; see for
example Refs. [18, 19, 20].

Figure 1: Floquet engineering of the long-range, anisotropic Heisenberg model
and measured dispersion relations. (a) The low-energy excitations of the
ferromagnetic state in the long-range, anisotropic Heisenberg model consist of
isolated magnons (red) and bound states thereof (blue). These excitations move
and interact with each other by the long-range spin exchange couplings. (b)
The Hamiltonian is experimentally realized from conventional long-range Ising
interactions (illustrated by blocks of XX and YY) by Floquet engineering,
which induces rotational symmetry in the XY plane of the spin and tunable ZZ
interactions. (c) Plane-wave spectroscopy of the single magnon excitation
energy indicates a cusp singularity at low momentum consistent with the
theoretical prediction (red dashed line). Inset: The group velocity $v(k)$
diverges in the long-wavelength limit. (d) Plane-wave spectroscopy of two-
magnon states measured for strong interactions $\Delta=3$ shows a well-defined
resonance at large momenta. Measurements agree well with theoretical
computations (red dashed line). Inset: The velocity of the magnon bound state
stays finite even when approaching the region of instability at small momenta,
indicated by the gray area.

The experiments described below are realized in a linear string of twenty
${}^{40}\mathrm{Ca}^{+}$ ions confined in a linear Paul trap. Qubits are
encoded in the Zeeman levels
$\left|\downarrow\right>=\left|S_{1/2},m=+1/2\right>$ and
$\left|\uparrow\right>=\left|D_{5/2},m=+5/2\right>$ and manipulated with a
narrow-linewidth laser at 729 nm. The ion crystal is laser-cooled by Doppler
and polarization gradient cooling and its transverse modes are prepared close
to their ground states by resolved sideband cooling. The effective spin-spin
interaction is characterized by a power-law exponent of $\alpha=1.4$ and a
maximum coupling strength of $J=369$ rad/s. Further experimental details are
provided in the Supplementary Information [21]. For the realization of the
long-range, anisotropic Heisenberg model, we use a Floquet decomposition of
the evolution operator; Fig. 1(b); see also [22]. The anisotropy $\Delta$ is
controlled by setting the relative length of the Floquet substeps. To reduce
the impact of dephasing noise, we implement the basis rotation between the
Floquet substeps as dynamical decoupling pulses. Furthermore, the
inhomogeneous coupling of the Gaussian laser beam, that induces collective
spin rotations, causes rotation errors. These errors are mitigated by
alternating the direction of the basis rotations. The laser pulses of the
individual Floquet steps are shaped to reduce spectral overlap with the radial
motional modes. Related Floquet protocols which enhance the symmetries of the
model by driving have for example recently been developed for Rydberg atoms
[23, 24] and trapped ions [22, 25].

Figure 2: Dynamics of magnon bound states. We prepare two magnons in the
center of the chain and study their evolution for interactions and maximum
time $\Delta=1.0$, $t_{\text{max}}J=5.5$ (left), $\Delta=2.0$,
$t_{\text{max}}J=4.1$ (middle), and $\Delta=3.5$, $t_{\text{max}}J=3.0$
(right), respectively. (a) Selection of three experimental snapshots measured
at time $t=(0,t_{\mathrm{max}}/2,t_{\mathrm{max}})$ from bottom to top. For
small interactions (left), magnons typically spread independently, whereas for
intermediate and strong interactions (middle, right), two magnons are
typically found nearby. (b) The one-magnon projector $P_{j}^{\uparrow}$
spreads quickly for small interactions. It significantly slows down for strong
interactions, for which it mainly captures the bound state dynamics, indicated
by a sharp light cone. (c) The two-magnon bound state (BS) projector
$P_{j,j+1}^{\uparrow\uparrow}$ decays rapidly for weak interaction (left) but
spreads with a light cone for larger interactions (middle, right) confirming
the robust formation of a bound state of two nearby magnons. (d) The
experimental measurements (symbols) for the two-magnon projector are in good
agreement with the theoretical predictions (solid lines). Data is obtained
from snapshots, which have been post-selected for two excitations. After post-
processing we retain $\approx(22\%,20\%,25\%)$ of typically 1500 measured
snapshots. Error bars are obtained from a jackknife analysis of the post-
processed snapshots.

The dispersion of a single magnon excitation does not depend on the spin
interaction $\Delta$, except for an overall constant offset. Using plane-wave
spectroscopy introduced in Ref. [11] we characterize the dispersion
$\varepsilon_{1}(k)$ of a single magnon (see Methods for details); Fig. 1 (c).
For all momenta in the Brillouin zone we measure a resonance, indicating the
existence of a well-defined quasi-particle excitation. We find that the magnon
velocity $v(k)=\partial_{k}\varepsilon_{1}(k)$ diverges in the long-wavelength
limit, which becomes manifest by the cusp of the dispersion at zero momentum.
The single-magnon dispersion can be readily computed analytically
$\varepsilon_{1}(k)=\frac{4J}{3}\sum_{\ell=1}^{\infty}\frac{\cos(k\ell)-\Delta}{\ell^{\alpha}}$
and is in good agreement with the measured data. This solution also shows that
the maximal magnon velocity is unbounded for $1<\alpha\leq 2$. The
experimentally studied case of $\alpha=1.4$ is well within this range. As a
consequence of the diverging group velocity, a local magnon excitation
propagates arbitrarily fast.

Figure 3: Characterizing magnon bound states. (a) For small interactions
($\Delta\\!=\\!1.0$) the bound state (BS) participation vanishes in time,
resembling a random distribution of magnons (see Methods). By contrast, stable
plateaus are reached at stronger interactions establishing the robust
existence of the bound state. (b) The bound state participation at time
$tJ=2.0$ [orange dashed line in (a)] unveils a crossover from a weak-
interaction regime without bound states to a strong-interaction regime with
magnon bound state at an interaction strength of $\Delta\approx 2.2$. For
comparison we show numerical data of the bound-state participation in larger
systems evaluated at late times (gray dashed line). (c) For large
interactions, we determine the velocity of the dominant bound state from the
light cone spreading. Data is obtained from typically 1000 snapshots, which we
post-select for two excitations. We retain $\approx 22\%$ of the snapshots at
the latest times. Error bars are obtained from a jackknife analysis of the
post-processed snapshots.

Given the unconventional properties of the single-magnon excitation, in
particular their unbounded velocity, it is _a priori_ unclear whether the
model supports bound states of two magnons. We measure the spectrum of two-
magnon excitations by initializing the system approximately in a plane wave of
two magnons situated next to each other, which can be achieved by
perturbatively acting with long-range Ising interactions for a short period
and tomographically reconstructing the two-site density matrix after an
evolution with the long-range, anisotropic Heisenberg model (see Methods). The
results of this measurement is shown in Fig. 1 (d) for strong magnon
interaction $\Delta=3$. For low momenta, $k\lesssim\pi/4$, the spectral weight
is distributed over a large energy window. By contrast, for high momenta,
$k\gtrsim\pi/4$ a sharp resonance emerges, which is indicative of the
formation of a bound state.

We compare the measurement signal with the dispersion of the magnon bound
state $\varepsilon_{2}(k)$, that we obtain from solving the two-body problem
numerically, and find good agreement. The following features of the
measurement should be emphasized: On the one hand, the bound state dispersion
law does not stretch across the whole Brillouin zone, but becomes unstable
upon approaching the long-wavelength limit. On the other hand, the group
velocity of the composite two-magnon excitation remains always finite, in
particular also when approaching the unstable regime at low momenta. As a
consequence, the bound state of two magnons should possess a well-defined
light cone. This is in stark contrast with the single magnon case.

In order to establish the well-defined light cone of the magnon bound state,
we prepare an initial state with two flipped spins situated next to each other
in the center of the spin chain. Subsequently, we probe the non-equilibrium
dynamics of this initial state under the evolution of Hamiltonian (1) for
different values of the interaction $\Delta$. The snapshots measured at times
$t=(0,t_{\mathrm{max}}/2,t_{\mathrm{max}})$ already exhibit striking
signatures of the bound state formation; Fig. 2 (a): At weak interactions
$\Delta=1.0$, the two magnons propagate freely and are typically measured at
uncorrelated positions in space. By contrast, when increasing the interactions
to $\Delta=2.0$ or even to $\Delta=3.5$, we find them with high probability
next to each other.

To quantitatively analyze the situation, we measure the single-magnon
$P^{\uparrow}_{j}$ and two-magnon projectors $P^{\uparrow\uparrow}_{j,j+1}$ by
averaging over a collection of typically $1500$ snapshots. The main
contribution to the bound state results from two magnons located nearby, which
is why we evaluate the two-magnon projector on adjacent sites in the spin
chain (see Supplementary Information [21]). As a consequence of the diverging
group velocity, single-magnon excitations propagate faster than linear in time
in the thermodynamic limit. For the accessible system size, this results in a
broad distribution of the one-magnon projector for small interactions
($\Delta\\!=\\!1.0$). The two-magnon projector, by contrast, quickly decays in
time and does not carry substantial weight in the wave function, which
indicates the absence of a magnon bound state; Fig. 2 (left column).

When increasing the interactions to $\Delta=2.0$ or even to $\Delta=3.5$ a
distinct behavior is observed; Fig. 2 (middle and right column, respectively).
First, the one-magnon projector spreads slower than for $\Delta=1.0$. This is
incompatible with having a significant contribution of single unbound magnons
in the wave function, as their dynamics are insensitive to the value of the
interaction $\Delta$. As a consequence, the observed signal can be attributed
to the formation of a magnon bound state. This is further confirmed by
measuring the two-magnon projector, which displays a well-defined light cone
with a velocity that is consistent with the one obtained from the single
magnon projector. To further support that the effective time evolution of the
Floquet protocol describes the long-range, anisotropic Heisenberg model, we
compare time slices of the experimental data to numerical simulations and find
excellent agreement, Fig. 2 (d).

Figure 4: Dynamics of mutual information. We realize two distinct initial
configurations, (a) one with two magnons next to each other and (b) one with
two magnons separated in space. We characterize the spreading of
configurational mutual information between segments A and B each consisting of
three consecutive spins. (c) For weak interactions ($\Delta\\!=\\!0.5$) the
configurational mutual information $I_{c}$ yields similar results for both
initial configurations. (d) For strong interactions ($\Delta\\!=\\!4.5$), a
significant difference in the growth of the mutual information is observed for
the two initial configurations. Insets: Late-time dynamics and comparison of
the measured configurational mutual information $I_{c}$ with the conventional
mutual information $I$ (see text for details). Data is obtained from typically
2200 snapshots, which we post-select for two excitations. We retain $\approx
19\%$ of the snapshots at the latest times. Error bars are obtained from a
jackknife analysis of the post-processed snapshots.

We now analyze the magnon bound state formation for a large range of
interaction strengths $\Delta$, by measuring total two-magnon projector
$\sum_{j}P^{\uparrow\uparrow}_{j,j+1}$. This quantity not only captures the
bound state but also two unbound magnons that accidentally occupy neighboring
lattice sites. To correct for this effect we introduce the bound state
participation, in which we normalize the two-magnon projector such that it
discards contributions arising from a random placement of two magnons on
neighboring sites. The random placement occurs with probability $2/L$, which
is sizeable for the 20-ion chain considered here (see Methods). This bound-
state participation then reaches a finite value at late times, provided a
magnon bound state exists, and decays to zero otherwise.

In agreement with the light-cone measurements, we find a vanishing bound state
participation for $\Delta\\!=\\!1.0$ at late times, and a saturation for
larger values of the interaction $\Delta$; Fig. 3 (a). The bound state
contribution increases continuously with increasing interactions; Fig. 3 (b)
and is in agreement with our numerical computations. The crossover in the
bound state participation becomes increasingly sharp, as we increase the
system size. In the Supplementary Information [21] we study the finite-size
dependence and find a crossover in the spectrum of the model: bound states
exist only for $\Delta\gtrsim 2.2$. This is consistent with recently reported
theoretical results [12]. We extract the light-cone velocity $v_{\mathrm{BS}}$
of the bound state from the propagation front of the two-magnon projector,
which we identify by the half-maximum of the signal; Fig. 3 (c).

In systems with well-defined quasi-particles, the growth of entanglement
following a quantum quench can be understood as follows: The initial state
creates entangled quasi-particle excitations that propagate ballistically
through the the system [26, 27]. In the long-range, anisotropic Heisenberg
model two competing effects arise. On the one hand, single magnons spread
arbitrarily fast, which could lead to entanglement growth that is faster than
linear in time. On the other hand, composite bound states of magnons form,
which strongly modify the excitation content of the state. In order to
elucidate the consequences of the bound-state formation on the entanglement
growth, we study the dynamics of two distinct initial states: one with two
magnons situated next to each other [Fig. 4 (a)] that have a large overlap
with the bound state provided the interaction is large enough, and one with
two magnons that are separated by one lattice site [Fig. 4 (b)]. We then study
the time evolution of this state under the long-range, anisotropic Heisenberg
model and quantify the entanglement between two separated intervals $A$ and
$B$.

One way of quantifying the entanglement between two separated subregions that
are embedded in a larger environment is achieved by the mutual information
$I=S_{A}+S_{B}-S_{A\cup B}$, where $S_{A}$ is the von Neumann entanglement
entropy between the interval $A$ and the rest of the system $\bar{A}$. The
entanglement entropy can be measured tomographically [28, 7] and its Rényi
variants by randomized measurements [29] and beam splitter operations [30].
However, proxies for the mutual information can already be obtained from the
snapshots in one computational basis. Following Ref. [31], the basic idea is
to observe that for systems which conserve the number of excitations $N$, the
total density matrix has a block diagonal form
$\hat{\rho}=\sum_{n=0}^{N}p(n)\hat{\rho}_{A}^{(n)}\otimes\hat{\rho}_{\bar{A}}^{(N-n)}$,
where $p(n)$ is the probability of having $n$ magnons in the subsystem $A$,
and $\hat{\rho}^{(n)}_{A}$ is the (normalized) reduced density matrix of the
subsystem $A$ in the $n$-magnon sector. Then the total entanglement entropy
$S_{A}$ can then be split into a number entropy $S_{N;A}$ and a
configurational entropy $S_{C;A}$

$S_{A}=S_{N;A}+S_{C;A},$ (2)

where $S_{N;A}=-\sum_{n=0}^{N}p(n)\log(p(n))$ and
$S_{C;A}=-\sum_{n=0}^{N}p(n)\text{Tr}(\hat{\rho}^{(n)}_{A}\log\hat{\rho}^{(n)}_{A})$
[32, 33]. Measuring the configurational entropy exactly has an exponential
cost. Hence, we approximate it for short times and small correlations between
$A$ and $\bar{A}$ by the mutual configurational entropy [31, 34]

$S_{C;A}\simeq\sum_{n}p(n)\sum_{\\{A_{n}\\},\\{\bar{A}_{n}\\}}\big{[}p(A_{n}\otimes\bar{A}_{n})-p(A_{n})p(\bar{A}_{n})\big{]}$
(3)

where $\\{A_{n}\\}$ $(\\{\bar{A}_{n}\\})$ is the set of all the possible
magnon configurations in the local basis in the subsystem $A$ $(\bar{A})$.
Using this approximation yields the proxy $I_{c}$ for the mutual information,
which we measure in the experiment.

For weak interactions, bound states cannot form. Hence, the initial state with
two magnons next to each other and the one in which the magnons are separated
by one site should give rise to comparable entanglement dynamics, as in both
cases magnons propagate freely. We find that both initial states give rise to
comparable dynamics of the configurational mutual information $I_{c}$; Fig. 4
(c), confirming the independent spreading of magnons. In the inset, we compare
the theoretically computed configurational mutual information $I_{c}$ with the
mutual information $I$ obtained from the von Neumann entanglement entropy and
find good agreement between the two.

For strong interactions the situation drastically changes. In that case, a
robust magnon bound state forms, which reduces the growth of entanglement, as
the bound state behaves effectively as a single particle. By contrast, when
initializing the dynamics with two magnons separated by one lattice site, the
magnons remain mainly unbound and move freely through the system, which
increases the dynamically accessible number of configurations, hence leading
to higher entanglement; Fig. 4 (d). The insets of Fig. 4 (c,d) confirm that
the observed behavior persists at late times: the mutual information $I$
saturates to a comparable value for weak interactions but is distinct for
strong interactions, demonstrating the impact of magnon bound states on the
entanglement dynamics.

By Floquet engineering, we have realized the long-range, anisotropic
Heisenberg model in a linear crystal of 20 ions. The periodic drive enhances
the discrete Ising interactions realized by sideband modulation to continuous
Heisenberg-type interactions. This enabled us to study the interaction between
magnons created on top of a ferromagnetic state. We experimentally found
robust bound states for strong interactions.

The formation of bound states can have direct implications on the relaxation
dynamics following a quantum quench in which a finite density of excitations
is created. In our system bound states will be stable against few-body
scattering processes for large interactions $\Delta$ due to energy
conservation. We expect that in this regime the relaxation will proceed in
multiple stages: First, the system relaxes to an ensemble in which bound
states act as well-defined quasiparticles. Second, bound states can decay due
to higher order processes giving rise to full thermalization. Other ergodic
systems featuring long-lived composite particles, such as mesons in confined
spin chains [35, 36, 37] or doublons protected by strong interactions in
Hubbard models [38, 39, 40, 41], are expected to undergo a similar prethermal
relaxation dynamics. This multi-stage dynamics should also manifest itself in
the entanglement growth as indicated by our measurements on mutual information
spreading. In future work, it will be interesting to investigate such multi-
stage relaxation dynamics in systems which support long-lived bound states.

Methods

### One-magnon spectroscopy.—

In order to examine the dispersion law of single magnon excitations
$\varepsilon_{1}(k)$ we use plane-wave spectroscopy, introduced in Ref. [11].
We start by creating a plane-wave superposition state of single magnon
excitations with given momentum on top of the ferromagnetic configuration
$\ket{0}=\ket{\downarrow\downarrow...\downarrow\downarrow}$. Plane-wave
initial states can be readily implemented using a sequence of modulated single
qubit rotations mixing both locally accessible states of the computational
spin basis. For our purpose such a rotation can be generated by individual
single-ion rotations around the X-axis,

$\ket{\psi_{1}}=e^{i\sum_{j}\gamma\mathcal{A}^{(k)}_{j}\hat{\sigma}_{j}^{x}}\ket{0}=\ket{0}+i\gamma\sum_{j}\mathcal{A}^{(k)}_{j}\hat{\sigma}_{j}^{x}\ket{0}+\mathcal{O}(\gamma^{2}).$
(4)

Our goal is to excite plane waves of the form
$\ket{k}=\frac{1}{\sqrt{L}}\sum_{j}e^{ikj}\hat{\sigma}_{j}^{x}\ket{0}$. For
the open boundary conditions realized in the experiment the momentum is
quantized as $k=\pi n/(L+1)$ for $n=1,\ldots,L$. To this end, we choose the
rotation angle to create a superposition of counter-propagating waves with
$\mathcal{A}^{(k)}_{j}=\sqrt{\frac{2}{L}}\sin(kj)$. In the measurements, we
chose $\gamma=0.7$ and then postprocess for single magnon excitations. This
implementation scheme allows us to create arbitrary linear combinations of
single magnon states. For example, a superposition of two standing waves with
different momenta $k$ and $q$ can be realized by fixing
$\mathcal{A}^{(k)}_{j}=\sqrt{\frac{2}{L}}(\sin(kj)+\sin(qj))$, from which
relative frequencies between the excitations at momentum $k$ and $q$ are
obtained. This has the major advantage, that measurements are performed within
a single magnon number sector, and hence the experiment operates in the
subspace that is free of decoherence between different magnon number sectors.

Having implemented a suitable initial state, we continue by evolving it under
the long-range, anisotropic Heisenberg Hamiltonian Eq. (1). Because the
dispersion law of a single magnon is expected to be independent of the chosen
interaction $\Delta$, we consider the experimentally simplest case of
vanishing interaction $\Delta=0$. Next, we perform spectroscopic measurements
of the magnon projector $P_{j}^{\uparrow}=\frac{1}{2}(\hat{\sigma}_{j}^{z}+1)$
in the evolved state. The measured signal is

$\displaystyle\bra{\psi_{1}(t)}P_{j}^{\uparrow}\ket{\psi_{1}(t)}=\gamma^{2}\;2\Re\Big{[}e^{-it(\varepsilon_{1}(k)-\varepsilon_{1}(q))}$
$\displaystyle\times\big{(}\bra{q}P_{j}^{\uparrow}\ket{k}+\bra{-q}P_{j}^{\uparrow}\ket{k}\big{)}\Big{]}+\mathcal{O}(\gamma^{3}),$
(5)

where non-oscillating contributions are dropped. We furthermore took advantage
of the symmetry of the dispersion law under sign changes of the momentum. The
energy difference $\varepsilon_{1}(k)-\varepsilon_{1}(q)$ is then obtained by
a Fourier transform. The results of Fig. 1 (c) correspond to a choice of the
reference momentum $q=\pi/(L+1)$ as the standing wave of lowest frequency
supported by the chain.

### Two-magnon spectroscopy.—

A straightforward generalization of the one-magnon spectroscopy to two magnons
requires the coherent manipulation of pairs of spins, which is difficult in
our experimental setting. To solve this challenge we proceed by creating a
coherent superposition of two-magnon states utilizing the short-time evolution
of an Ising Hamiltonian
$\sum_{j<j^{\prime}}J_{jj^{\prime}}\sigma_{j}^{x}\sigma_{j^{\prime}}^{x}$ over
a period $tJ=0.19$. The coupling matrix decays as
$J_{jj^{\prime}}=J/|j^{\prime}-j|^{\alpha}$. As a result we will predominantly
obtain two-magnon states with both magnons situated at neighboring sites.
Configurations with larger separation $d$ of the magnons will be suppressed.
To realize a plane-wave initial state of two close-by magnons we have to
furthermore tune the relative phase between the individual configurations.
This can be achieved by applying a sequence of single qubit rotations with
angles $\\{\phi_{j}\\}$ from the left-most to the right-most ion. As a result,
we realize a two-magnon initial state

$\ket{\psi_{2}}=e^{-i\sum_{j}\frac{\phi_{j}}{2}\hat{\sigma}_{j}^{z}}e^{-i\gamma\sum_{j<j^{\prime}}J_{jj^{\prime}}\hat{\sigma}^{x}_{j}\hat{\sigma}^{x}_{j^{\prime}}}\ket{0}\approx
e^{i\sum_{j}\frac{\phi_{j}}{2}}\\\ \times\Big{(}\ket{0}-i\gamma
J_{12}\sum_{j}e^{-i(\phi_{j}+\phi_{j+1})}\ket{...\uparrow_{j}\uparrow_{j+1}...}\Big{)}+\mathcal{O}(\gamma^{2})\,.$
(6)

By judiciously choosing the phase imprinting $\phi_{j}+\phi_{j+1}=jk$, we can
target the desired plane wave. This state is then evolved with the long-range,
anisotropic Heisenberg Hamiltonian (1). We obtain the energy of the two-magnon
plane-wave configuration by measuring
$\hat{\sigma}^{-}_{j}\hat{\sigma}^{-}_{j+1}$, where
$\hat{\sigma}^{-}_{j}=(\hat{\sigma}^{x}_{j}-i\hat{\sigma}^{y}_{j})/2$. This
observable connects configurations of two magnons located next to each other
with the ferromagnetic state $\ket{0}$

$\displaystyle\bra{\psi_{2}(t)}\hat{\sigma}^{-}_{j}\hat{\sigma}^{-}_{j+1}\ket{\psi_{2}(t)}=-i\gamma\Gamma
J_{12}e^{-ijk}e^{-i\big{(}\varepsilon_{2}(k)-\varepsilon_{0}\big{)}t}$
$\displaystyle+\;\mathcal{O}(\gamma^{2}).$ (7)

The measurement signal hence oscillates with the desired frequency of the two-
magnon excitation. The contrast of the oscillation frequency carries also
information about the probability amplitude $\Gamma$ of the bound state in the
total many-body wave function. Performing a Fourier transformation in time and
averaging over sites 8 to 13 in the center of the chain to reduce
contributions arising from the boundary, we obtain the dispersion law of the
two-magnon bound state with respect to the ferromagnetic configuration
$\varepsilon_{0}-\varepsilon_{2}(k)$, illustrated in Fig. 1 (d).

### Bound state participation.—

The two-magnon projector $P^{\uparrow\uparrow}_{j,j+1}$ is sensitive to bound
states, as they consist of correlated magnons at a short relative distance
(see Supplementary Information [21]). When two unbound magnons are randomly
placed on a chain of $L$ sites, a non-vanishing expectation value for the
projector will arise statistically $\langle
P^{\uparrow\uparrow}_{j,j+1}\rangle_{\text{random}}=2/L$. To account for this,
we properly renormalize the bound state participation shown in Fig. 3 as

$\text{BS participation}=\frac{1}{1-2/L}\Bigg{(}\sum_{j=1}^{L-1}\langle
P_{j,j+1}^{\uparrow\uparrow}(t)\rangle-2/L\Bigg{)}.$ (8)

###### Acknowledgements.

Acknowledgements: We acknowledge support from the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy–EXC–2111–390814868 and DFG grants No. KN1254/1-2,
KN1254/2-1, and from the Austrian Science Fund through the SFB BeyondC
(F7110). Furthermore, The project leading to this application has received
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme (grant agreements No. 741541
and No. 851161), from the Munich Quantum Valley, which is supported by the
Bavarian state government with funds from the Hightech Agenda Bayern Plus, and
from the Institut für Quanteninformation GmbH. Data and Informations
availability: Raw data, data analysis, and simulation codes are available on
Zenodo [42].

Author contributions. The research was devised by MK, and CFR. SB, AB, and MK
developed the theoretical analysis. MJ, FK, RB, and CFR contributed to the
experimental setup. FK and MJ performed the experiments. FK, SB, AB, MK and
CFR analyzed the data and SB carried out numerical simulations. FK, SB, AB,
MK, and CFR wrote the manuscript. All authors contributed to the discussion of
the results and the manuscript. All correspondence should be addressed to MK
<EMAIL_ADDRESS>and CFR (christian.roos@uibk.ac.at).

Competing interests. The authors declare no competing interests.

## References

  * Bethe [1931] H. Bethe, Zur Theorie der Metalle, Zeitschrift für Physik 71, 205 (1931).
  * Auerbach [2012] A. Auerbach,  _Interacting electrons and quantum magnetism_ (Springer Science & Business Media, 2012).
  * Ganahl  _et al._ [2012] M. Ganahl, E. Rabel, F. H. L. Essler, and H. G. Evertz, Observation of Complex Bound States in the Spin-$1/2$ Heisenberg $XXZ$ Chain Using Local Quantum Quenches, Phys. Rev. Lett. 108, 077206 (2012).
  * Fukuhara  _et al._ [2013a] T. Fukuhara, A. Kantian, M. Endres, M. Cheneau, P. Schauß, S. Hild, D. Bellem, U. Schollwöck, T. Giamarchi, C. Gross, I. Bloch, and S. Kuhr, Quantum dynamics of a mobile spin impurity, Nature Physics 9, 235 (2013a).
  * Fukuhara  _et al._ [2013b] T. Fukuhara, P. Schauß, M. Endres, S. Hild, M. Cheneau, I. Bloch, and C. Gross, Microscopic observation of magnon bound states and their dynamics, Nature 502, 76 (2013b).
  * Morvan  _et al._ [2022] A. Morvan, T. I. Andersen, X. Mi, C. Neill, A. Petukhov, K. Kechedzhi, D. Abanin, R. Acharya, F. Arute, K. Arya, _et al._ , Formation of robust bound states of interacting photons, arXiv:2206.05254 (2022).
  * Jurcevic  _et al._ [2014] P. Jurcevic, B. P. Lanyon, P. Hauke, C. Hempel, P. Zoller, R. Blatt, and C. F. Roos, Quasiparticle engineering and entanglement propagation in a quantum many-body system, Nature 511, 202 (2014).
  * Richerme  _et al._ [2014] P. Richerme, Z.-X. Gong, A. Lee, C. Senko, J. Smith, M. Foss-Feig, S. Michalakis, A. V. Gorshkov, and C. Monroe, Non-local propagation of correlations in quantum systems with long-range interactions, Nature 511, 198 (2014).
  * Monroe  _et al._ [2021] C. Monroe, W. C. Campbell, L.-M. Duan, Z.-X. Gong, A. V. Gorshkov, P. W. Hess, R. Islam, K. Kim, N. M. Linke, G. Pagano, P. Richerme, C. Senko, and N. Y. Yao, Programmable quantum simulations of spin systems with trapped ions, Rev. Mod. Phys. 93, 025001 (2021).
  * Hauke and Tagliacozzo [2013] P. Hauke and L. Tagliacozzo, Spread of Correlations in Long-Range Interacting Quantum Systems, Phys. Rev. Lett. 111, 207202 (2013).
  * Jurcevic  _et al._ [2015] P. Jurcevic, P. Hauke, C. Maier, C. Hempel, B. P. Lanyon, R. Blatt, and C. F. Roos, Spectroscopy of Interacting Quasiparticles in Trapped Ions, Phys. Rev. Lett. 115, 100501 (2015).
  * Macrì  _et al._ [2021] T. Macrì, L. Lepori, G. Pagano, M. Lewenstein, and L. Barbiero, Bound state dynamics in the long-range spin-$\frac{1}{2}$ XXZ model, Phys. Rev. B 104, 214309 (2021).
  * Porras and Cirac [2004] D. Porras and J. I. Cirac, Effective Quantum Spin Systems with Trapped Ions, Phys. Rev. Lett. 92, 207901 (2004).
  * Kokail  _et al._ [2019] C. Kokail, C. Maier, R. van Bijnen, T. Brydges, M. K. Joshi, P. Jurcevic, C. A. Muschik, P. Silvi, R. Blatt, C. F. Roos, and P. Zoller, Self-verifying variational quantum simulation of lattice models, Nature 569, 355 (2019).
  * Joshi  _et al._ [2022] M. K. Joshi, F. Kranzl, A. Schuckert, I. Lovas, C. Maier, R. Blatt, M. Knap, and C. F. Roos, Observing emergent hydrodynamics in a long-range quantum magnet, Science 376, 720 (2022).
  * Feng  _et al._ [2022] L. Feng, O. Katz, C. Haack, M. Maghrebi, A. V. Gorshkov, Z. Gong, M. Cetina, and C. Monroe, Continuous Symmetry Breaking in a Trapped-Ion Spin Chain, arXiv:2211.01275 (2022).
  * Bukov  _et al._ [2015] M. Bukov, L. D’Alessio, and A. Polkovnikov, Universal high-frequency behavior of periodically driven systems: from dynamical stabilization to Floquet engineering, Advances in Physics 64, 139 (2015).
  * Arrazola  _et al._ [2016] I. Arrazola, J. S. Pedernales, L. Lamata, and E. Solano, Digital-Analog Quantum Simulation of Spin Models in Trapped Ions, Scientific Reports 6, 30534 (2016).
  * Bermudez  _et al._ [2017] A. Bermudez, L. Tagliacozzo, G. Sierra, and P. Richerme, Long-range Heisenberg models in quasiperiodically driven crystals of trapped ions, Phys. Rev. B 95, 024431 (2017).
  * Birnkammer  _et al._ [2022a] S. Birnkammer, A. Bohrdt, F. Grusdt, and M. Knap, Characterizing topological excitations of a long-range Heisenberg model with trapped ions, Phys. Rev. B 105, L241103 (2022a).
  * sup [2022] Supplementary Material for Observation of magnon bound states in the long-range, anisotropic Heisenberg model (2022).
  * Kranzl  _et al._ [2022] F. Kranzl, A. Lasek, M. K. Joshi, A. Kalev, R. Blatt, C. F. Roos, and N. Yunger Halpern, Experimental observation of thermalisation with noncommuting charges, arXiv:2202.04652 (2022).
  * Geier  _et al._ [2021] S. Geier, N. Thaicharoen, C. Hainaut, T. Franz, A. Salzinger, A. Tebben, D. Grimshandl, G. Zürn, and M. Weidemüller, Floquet Hamiltonian engineering of an isolated many-body spin system, Science 374, 1149–1152 (2021).
  * Scholl  _et al._ [2022] P. Scholl, H. J. Williams, G. Bornet, F. Wallner, D. Barredo, L. Henriet, A. Signoles, C. Hainaut, T. Franz, S. Geier, A. Tebben, A. Salzinger, G. Zürn, T. Lahaye, M. Weidemüller, and A. Browaeys, Microwave Engineering of Programmable $XXZ$ Hamiltonians in Arrays of Rydberg Atoms, PRX Quantum 3, 020303 (2022).
  * Morong  _et al._ [2022] W. Morong, K. S. Collins, A. De, E. Stavropoulos, T. You, and C. Monroe, Engineering dynamically decoupled quantum simulations with trapped ions, arXiv:2209.05509 (2022).
  * Calabrese and Cardy [2005] P. Calabrese and J. Cardy, Evolution of entanglement entropy in one-dimensional systems, Journal of Statistical Mechanics: Theory and Experiment 2005, P04010 (2005).
  * Alba and Calabrese [2017] V. Alba and P. Calabrese, Entanglement and thermodynamics after a quantum quench in integrable systems, Proceedings of the National Academy of Sciences 114, 7947 (2017).
  * Roos  _et al._ [2004] C. F. Roos, G. P. T. Lancaster, M. Riebe, H. Häffner, W. Hänsel, S. Gulde, C. Becher, J. Eschner, F. Schmidt-Kaler, and R. Blatt, Bell States of Atoms with Ultralong Lifetimes and Their Tomographic State Analysis, Phys. Rev. Lett. 92, 220402 (2004).
  * Brydges  _et al._ [2019] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch, C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F. Roos, Probing Rényi entanglement entropy via randomized measurements, Science 364, 260 (2019).
  * Islam  _et al._ [2015] R. Islam, R. Ma, P. M. Preiss, M. E. Tai, A. Lukin, M. Rispoli, and M. Greiner, Measuring entanglement entropy in a quantum many-body system, Nature 528, 77 (2015).
  * Lukin  _et al._ [2019] A. Lukin, M. Rispoli, R. Schittko, M. E. Tai, A. M. Kaufman, S. Choi, V. Khemani, J. Léonard, and M. Greiner, Probing entanglement in a many-body–localized system, Science 364, 256 (2019).
  * Wiseman and Vaccaro [2003] H. M. Wiseman and J. A. Vaccaro, Entanglement of Indistinguishable Particles Shared between Two Parties, Phys. Rev. Lett. 91, 097902 (2003).
  * Goldstein and Sela [2018] M. Goldstein and E. Sela, Symmetry-Resolved Entanglement in Many-Body Systems, Phys. Rev. Lett. 120, 200602 (2018).
  * Wolf  _et al._ [2008] M. M. Wolf, F. Verstraete, M. B. Hastings, and J. I. Cirac, Area Laws in Quantum Systems: Mutual Information and Correlations, Phys. Rev. Lett. 100, 070502 (2008).
  * Kormos  _et al._ [2017] M. Kormos, M. Collura, G. Takács, and P. Calabrese, Real-time confinement following a quantum quench to a non-integrable model, Nature Physics 13, 246 (2017).
  * Lerose  _et al._ [2020] A. Lerose, F. M. Surace, P. P. Mazza, G. Perfetto, M. Collura, and A. Gambassi, Quasilocalized dynamics from confinement of quantum excitations, Phys. Rev. B 102, 041118 (2020).
  * Birnkammer  _et al._ [2022b] S. Birnkammer, A. Bastianello, and M. Knap, Prethermalization in confined spin chains, arXiv:2202.12908 (2022b).
  * Winkler  _et al._ [2006] K. Winkler, G. Thalhammer, F. Lang, R. Grimm, J. Hecker Denschlag, A. J. Daley, A. Kantian, H. P. Büchler, and P. Zoller, Repulsively bound atom pairs in an optical lattice, Nature (London) 441, 853 (2006).
  * Strohmaier  _et al._ [2010] N. Strohmaier, D. Greif, R. Jördens, L. Tarruell, H. Moritz, T. Esslinger, R. Sensarma, D. Pekker, E. Altman, and E. Demler, Observation of Elastic Doublon Decay in the Fermi-Hubbard Model, Phys. Rev. Lett. 104, 080401 (2010).
  * Sensarma  _et al._ [2010] R. Sensarma, D. Pekker, E. Altman, E. Demler, N. Strohmaier, D. Greif, R. Jördens, L. Tarruell, H. Moritz, and T. Esslinger, Lifetime of double occupancies in the Fermi-Hubbard model, Phys. Rev. B 82, 224302 (2010).
  * Abanin  _et al._ [2017] D. Abanin, W. D. Roeck, W. W. Ho, and F. Huveneers, A rigorous theory of many-body prethermalization for periodically driven and closed quantum systems, Commun. Math. Phys. 354, 809 (2017).
  * zen [2022] Observation of magnon bound states in the long-range, anisotropic Heisenberg model, Zenodo 10.5281/zenodo.7357735 (2022).

Supplementary Information
Observation of magnon bound states in the long-range, anisotropic Heisenberg
model

##  1 Experimental platform

In our setup, linear crystals of twenty ${}^{40}\mathrm{Ca}^{+}$ ions are
trapped in a linear Paul trap with a low axial trapping frequency of 92 kHz
and radial trapping frequencies of 2.926 MHz and 2.894 MHz, respectively. The
qubit is encoded in the Zeeman levels
$\left|\downarrow\right>=\left|S_{1/2},m=+1/2\right>$ and
$\left|\uparrow\right>=\left|D_{5/2},m=+5/2\right>$ which can be manipulated
with a narrow-linewidth laser at 729 nm. Entanglement between the qubits is
generated by a bichromatic light field that off-resonantly couples the
collective motional modes from the direction perpendicular to the ion crystal.
The two frequency components of the laser beam are detuned by $\pm 40$ kHz
from the blue and the red sidebands of the radial center-of-mass mode of the
ion crystal, giving rise to a long-range effective Ising coupling
$H_{\text{XX}}=\sum_{i<j}J_{ij}\sigma_{i}^{x}\sigma_{j}^{x}$, that decays
approximately as $J_{ij}=J/|i-j|^{\alpha}$ with a maximum coupling strength of
$J=369$ rad/s and a power-law exponent of $\alpha=1.4$.

##  2 Floquet Hamiltonian realized with dynamical decoupling

The Floquet sequence uses dynamical decoupling to reduce the impact of
dephasing errors caused by variations in the ambient magnetic field.
Furthermore, global rotations are applied in an alternating fashion to reduce
the impact of systematic over- and under-rotations, respectively, caused by
inhomogeneous illumination by the global laser field. A numerical simulation
of the Floquet sequence shows the robustness against constant detuning errors
(Fig. S1). The dynamically decoupled sequence $U$ maintains a high fidelity
over a much larger detuning range than the simple Floquet sequence $U^{(0)}$.

As a consequence of the dynamical decoupling pulses, the basis of the quantum
state is rotated with each Floquet step. If the number $n$ of Floquet steps is
not an integer multiple of eight (periodicity of Floquet sequence), the basis
has to be rotated accordingly. The final rotation $R_{\mathrm{f}}$ is
(operation on the right applied first)
$R_{\mathrm{f},n=1}=R_{-x}(\pi/2)R_{y}(\pi/2)$,
$R_{\mathrm{f},n=2}=R_{-y}(\pi)R_{x}(\pi/2)$,
$R_{\mathrm{f},n=3}=R_{-x}(\pi/2)R_{-y}(\pi/2)$,
$R_{\mathrm{f},n=4}=R_{-x}(\pi/2)$,
$R_{\mathrm{f},n=5}=R_{-x}(\pi/2)R_{-y}(\pi/2)$,
$R_{\mathrm{f},n=6}=R_{-y}(\pi)R_{x}(\pi/2)$,
$R_{\mathrm{f},n=7}=R_{-x}(\pi/2)R_{y}(\pi/2)$, and
$R_{\mathrm{f},n=8}=R_{-x}(\pi/2)$. Here, $R_{\alpha}(\theta)$ denotes a
rotation around the axis $\alpha$ by an angle $\theta$.

For the experimental implementation of the Floquet sequence, an important
source of error is the resonant coupling of the entangling light field to the
radial motional modes of the ion crystal. As each Trotter substep is realized
by a finite-length light pulse, the light’s spectrum of the Floquet sequence
is Fourier-broadened, with distinct peaks spanning over several tens of
kilohertz. To reduce overlap with the radial motional modes, we use pulse
shaping for each Trotter substep (Blackman window, rise and fall time 18 µs,
respectively).

Figure S1: Floquet protocol under dephasing errors. (a) Floquet Hamiltonian
realized with dynamical decoupling. The series of Floquet steps (dashed boxes)
repeats after $n=8$ steps. The final rotation $R_{\mathrm{f}}$ depends on the
number of Trotter steps. Here, $R_{\alpha}(\theta)$ denotes a rotation around
the axis $\alpha$ by an angle $\theta$. (b) Floquet Hamiltonian without
dynamical decoupling. (c) We evaluate the robustness of the Floquet protocol
under a dephasing error $\delta/2\sum_{i}\sigma_{z}^{(i)}$, where $\delta$
denotes the detuning between laser and qubit. The numerical simulation is
carried out for a system size of 10 qubits, magnon interaction $\Delta=3.5$
and a power-law spin-spin coupling $J_{i,j}=J/\left|i-j\right|^{\alpha}$ with
$J=369$ rad/s and $\alpha=1.4$. The maximum evolution time is $tJ\approx 3.3$,
divided into 32 Trotter steps.

##  3 Numerical Simulations with Exact Diagonalization

The numerical results presented in this work are obtained using exact
diagonalization (ED) for the studied spin Hamiltonian $H$. We explicitly
implement the magnon number conservation of the model and perform numerically
exact simulations for comparatively large systems up to $L=600$ sites.

Figure S2: Theoretically characterizing bound states. (a) We show the
$L_{4}$-norm of the two-magnon eigenstates as a function of momentum $k$ and
interaction $\Delta$ for comparatively large systems of $L=600$. A finite
$L_{4}$-norm indicates the existence of a bound state, which arise in some
part of the Brillouin zone for $\Delta\gtrsim 2.2$. The momentum range in
which bound states exist increases with $\Delta$. Inset: $L_{4}$-norm for the
experimentally studied system size of $L=20$. (b) Examples of the eigenstates
for large interaction $\Delta=3.0$ (upper panel) and intermediate interaction
$\Delta=2.0$ (lower panel). When bound states exist (upper panel), the
dominant contributions to the wave functions results from neighboring sites in
the spin chain, as indicated by sharp peaks of
$|\psi(k,d)|^{2}/|\psi(k,1)|^{2}$.

##  4 Characterization of two-magnon bound states

The results presented in the main text indicated a non-trivial dependence of
bound states on the interaction $\Delta$ and momentum $k$. We now provide a
theoretical analysis of the eigenstates of the long-range, anisotropic
Heisenberg model in the two-magnon sector for different interaction strength
$\Delta$ and momenta $k$. A convenient strategy to detect localized
wavefunctions is by looking at their spatial support using the $L_{4}$ norm. A
finite $L_{4}$ norm in the thermodynamic limit indicates a bound state. Let
$\psi(k,d)$ be an eigenstate in the two-magnon sector with momentum $k$, and
$d$ characterizes the relative distance between the flipped spins. Then the
$L_{4}$ norm is defined as
$\lVert\psi(k)\rVert_{4}=\sum_{d=1}^{L-1}|\psi(k,d)|^{4}$.

In Fig. S2 (a) we show results of the $L_{4}$-norm for the state with highest
energy, representing the probable candidate for a bound state in the case of
repulsive magnon interactions ($\Delta>0$). For large systems ($L=600$) we
find well-separated regions with finite $L_{4}$-norm for $\Delta>2.2$. When
increasing the interaction strength, the regime in which the bound states
exist starts to extend over a larger momentum range but it always vanishes at
the peculiar point $k=0$, i.e., at the point where the single-magnon velocity
diverges. Simulations carried out with the experimentally realized system size
of $L=20$ yield qualitatively similar results, but a non-negligible
contribution persists for smaller interaction $\Delta$; see Fig. S2 (a)
(inset).

We study the decay of the bound state wave function with relative distance
between the magnons in Fig. S2 (b). Within the region featuring stable bound
states ($\Delta=3$, upper panel) we find that the bound state decays as a
power law with distance in space (inset). The bound state therefore mainly
occupies nearest neighbor sites, which motivates the choice of the two-site
projector of magnons on adjacent sites, which we have used in the main text.
For smaller interactions ($\Delta=2$, lower panel) strong tails emerge in the
wave function that do not decay to zero. Such wave functions give rise to a
vanishing $L_{4}$ norm and represent an extended state.