††thanks: Corresponding author.††thanks: Present address: CGnal S.r.l. Corso
Venezia 43, 20121 Milano, Italy

#  Theory of robust quantum many-body scars in long-range interacting systems

Alessio Lerose  0000-0003-1555-5327<EMAIL_ADDRESS>Department of
Theoretical Physics, University of Geneva, Quai Ernest-Ansermet 30, 1205
Geneva, Switzerland  Tommaso Parolini  The Abdus Salam International Centre
for Theoretical Physics (ICTP), Strada Costiera 11, 34151 Trieste, Italy
Scuola Internazionale di Studi Superiori Avanzati, Via Bonomea, 265, 34136
Trieste, Italy  Rosario Fazio  0000-0002-7793-179X The Abdus Salam
International Centre for Theoretical Physics (ICTP), Strada Costiera 11, 34151
Trieste, Italy  Dipartimento di Fisica, Università di Napoli “Federico II”,
Monte S. Angelo, I-80126 Napoli, Italy  Dmitry A. Abanin  Department of
Theoretical Physics, University of Geneva, Quai Ernest-Ansermet 30, 1205
Geneva, Switzerland  Department of Physics, Princeton University, Princeton,
New Jersey 08544, USA  Silvia Pappalardi  0000-0001-6931-8736 Institut für
Theoretische Physik, Universität zu Köln, Zülpicher Straße 77, 50937 Köln,
Germany

###### Abstract

Quantum many-body scars (QMBS) are exceptional energy eigenstates of quantum
many-body systems associated with violations of thermalization for special
non-equilibrium initial states. Their various systematic constructions require
fine-tuning of local Hamiltonian parameters. In this work we demonstrate that
the setting of _long-range_ interacting quantum spin systems generically hosts
_robust_ QMBS. We analyze spectral properties upon raising the power-law decay
exponent $\alpha$ of spin-spin interactions from the solvable permutationally-
symmetric limit $\alpha=0$. First, we numerically establish that despite
spectral signatures of chaos appear for infinitesimal $\alpha$, the towers of
$\alpha=0$ energy eigenstates with large collective spin are smoothly deformed
as $\alpha$ is increased, and exhibit characteristic QMBS features. To
elucidate the nature and fate of these states in larger systems, we introduce
an analytical approach based on mapping the spin Hamiltonian onto a
relativistic quantum rotor non-linearly coupled to an extensive set of bosonic
modes. We exactly solve for the eigenstates of this interacting impurity
model, and show their self-consistent localization in large-spin sectors of
the original Hamiltonian for $0<\alpha<d$. Our theory unveils the stability
mechanism of such QMBS for arbitrary system size and predicts instances of its
breakdown e.g. near dynamical critical points or in presence of semiclassical
chaos, which we verify numerically in long-range quantum Ising chains. As a
byproduct, we find a predictive criterion for presence or absence of heating
under periodic driving for $0<\alpha<d$, beyond existing Floquet-
prethermalization theorems. Broader perspectives of this work range from
independent applications of the technical toolbox developed here to informing
experimental routes to metrologically useful multipartite entanglement.

##  I Introduction

While Feynman’s vision of a universal quantum simulator fuelled early
developments of synthetic matter platforms Feynman  _et al._ (2018), it
quickly became clear that these setups allow to investigation a separate set
of fundamental questions on quantum matter out of equilibrium, hardly
conceivable in traditional condensed-matter setups. A primary question
concerns the mechanism of thermalization in isolated quantum many-particle
systems – a central tenet of macroscopic thermodynamics Polkovnikov  _et al._
(2011). The possibility of violating this behavior has recently attracted
great interest, as it underpins the ongoing quest for design and coherent
control of non-thermal states of matter. In this context several mechanisms
have been discovered and characterized, including various instances of many-
body localization phenomena Abanin  _et al._ (2019); De Roeck and Huveneers
(2014); Yao  _et al._ (2016); Schiulaz  _et al._ (2015); Michailidis  _et al._
(2018); Smith  _et al._ (2017); Brenes  _et al._ (2018); van Nieuwenburg  _et
al._ (2019); Schulz  _et al._ (2019); James  _et al._ (2019); Lerose  _et al._
(2020), prethermalization arising from weak integrability breaking Bertini
_et al._ (2015); Mallayya  _et al._ (2019); Abanin  _et al._ (2017); De Roeck
and Verreet (2019), and Hilbert space fragmentation Sala  _et al._ (2020);
Khemani  _et al._ (2020).

Certain quantum many-body systems may fail to thermalize only when initialized
in a specific class of simple far-from-equilibrium states, whereas all other
initial states give ordinary thermalization dynamics. This phenomenon has been
associated with the existence of anomalous highly excited energy eigenstates
known as _quantum many-body scars_ (QMBS) — a term coined in Ref. Turner  _et
al._ (2018a) by analogy to wavefunction scarring in single-particle
semiclassical machanics Heller (1984). While distinct mechanisms fall under
the umbrella of QMBS Serbyn  _et al._ (2021); Moudgalya  _et al._ (2022);
Chandran  _et al._ (2023), their common hallmark is a vanishing fraction of
eigenstates violating the strong eigenstate thermalization hypothesis (ETH)
Kim  _et al._ (2014); D’Alessio  _et al._ (2016), i.e. exhibiting non-thermal
values of local observables and anomalously low entanglement, embedded in an
otherwise ETH-compatible spectrum. The discovery of long-lived coherent
oscillations in Rydberg-atom quantum simulations of far-from-equilibrium
dynamics Bernien  _et al._ (2017) triggered a huge theoretical interest in
QMBS Turner  _et al._ (2018b); Shiraishi and Mori (2017); Ho  _et al._ (2019);
Khemani  _et al._ (2019a); Choi  _et al._ (2019); Iadecola  _et al._ (2019);
Schecter and Iadecola (2019); Mark  _et al._ (2020); Turner  _et al._ (2021);
Moudgalya and Motrunich (2022); Moudgalya  _et al._ (2018); Omiya and Müller
(2023); Gotta  _et al._ (2023); Hummel  _et al._ (2023); Surace  _et al._
(2023); Logarić  _et al._ (2023), and it is even more remarkable considering
that all the work thus far points to fragility of QMBS to pertubations Lin
_et al._ (2020); Surace  _et al._ (2021); Serbyn  _et al._ (2021); Moudgalya
_et al._ (2022); Chandran  _et al._ (2023). Identifying Hamiltonians with
_robust_ QMBS within an experimentally sensible subclass of interactions would
be an important breakthrough for both theory and experiments.

With very few exceptions Bull  _et al._ (2022); Desaules  _et al._ (2022a),
the vast majority of work on QMBS concerned local Hamiltonians. _Long-range_
interactions represent however a promising avenue of investigation: on the
experimental side, they are naturally present in several analog quantum-
simulation platforms such as trapped ions Blatt and Roos (2012); Britton  _et
al._ (2012); Richerme  _et al._ (2014), dipolar gases Chomaz  _et al._ (2022),
polar molecules Yan  _et al._ (2013), cold atoms in cavities Baumann  _et al._
(2010), and even solid-state spinful defects Kucsko  _et al._ (2018); on the
theory side, they are known to give rise to anomalous dynamical phenomena
Kastner (2011); Mori (2019); Neyenhuis  _et al._ (2017); Lerose  _et al._
(2018, 2019a); Defenu (2021); Lerose  _et al._ (2019b); Liu  _et al._ (2019);
Verdel  _et al._ (2020); Lerose and Pappalardi (2020); Defenu  _et al._
(2023). This behavior can be traced back to a peculiar form of integrability
breaking. In the limit of all-to-all infinite-range interactions the full
permutational symmetry allows to reduce the quantum many-body dynamics to the
semiclassical dynamics of few collective degrees of freedom. As the range of
interactions is decreased, all other degrees of freedom associated with
gradually shorter wavelengths get “activated” and interact with the collective
ones. This tunable decoupling gives rise a scenario that interpolates between
few-body and many-body physics. As a consequence, non-equilibrium behavior
shows resilience over long prethermal stages of dynamics Žunkovič  _et al._
(2018); Halimeh  _et al._ (2017); Santos  _et al._ (2016a); Liu  _et al._
(2019); Collura  _et al._ (2022), whose duration increases with the
interaction range and grows unbounded in the thermodynamic limit when
interactions decay slower than $1/r^{d}$ ($d=$ system dimensionality) Lerose
_et al._ (2019c); Mori (2019); Lerose and Pappalardi (2020). The occurrence of
thermalization at longer times remains, however, an open question.

In this paper we demonstrate the occurrence of _robust QMBS_ arising from
long-range interactions in isotropic quantum spin lattices. This finding
implies _thermalization breakdown at arbitrarily long times_ in this large
class of quantum many-body systems. As key point of our analysis, we unveil an
emergent exact solvability of eigenstates originating from high-permutational-
symmetry sectors in these systems. We show that such eigenstates get smoothly
deformed upon changing Hamiltonian parameters as long as the decay of spin-
spin interactions remains sufficiently slow. The energy eigenvalues associated
with such QMBS form _unequally_ spaced regular towers, labelled by the
collective spin quantum numbers and by the occupation numbers of quasiparticle
excitations emerging from our analytical construction. Together with the
numerical observation of quantum chaos signatures in level statistics for
arbitrary decaying interactions, our findings establish _QMBS as a generic
feature of long-range interacting spin lattices_. Our theory gives stringent
quantitative criteria for the stability of QMBS, which we numerically
illustrate in the variable-range quantum Ising chain.

##  II Overview

The results of this paper apply to general quantum spin lattices with spin-
spin interactions decaying algebraically as $1/r^{\alpha}$ with the distance
$r$. For the sake of concreteness we work with a variable-range Ising-like
quantum spin chain, introduced in Sec. III. The following Sections present the
main contributions of this work:

Figure 1:  Graphical summary of the theory of robust quantum many-body scars
developed in this work. (Left panel) The starting point is the solvable
semiclassical spectrum of the mean-field limit of interacting spin lattices
($\alpha=0$), labelled by the exact quantum numbers $S=sL^{d}-\delta_{S}$
(collective spin size) and $N$ (quantized classical action). Spin depletion
$\delta_{S}$ is associated with increasing degeneracy, which can be resolved
in terms of magnon excitations with non-zero momentum. (Central panel) We
express the finite-range spin Hamiltonian with $\alpha\neq 0$ in terms of
these degrees of freedom. This yields a reformulation as an effective quantum
impurity model: a relativistic quantum rotor parametrically coupled to an
ensemble of bosons. In the dilute sector $\delta_{S}\ll L/2$ of a given energy
shell the effective rotor-magnon Hamiltonian simplifies to the extent that we
can find a (non-trivial) exact solution for its spectrum and eigenstates.
(Right panel) Self-consistency of this solution is satisfied only for
$0<\alpha<d$, in agreement with our numerical data analysis for small system
size. Our theory unveils the mechanism of eigenstate localization underlying
the robust quantum many-body scars. It also predicts and classifies instances
of its breakdown, which we verify numerically.  Figure 2:  Relation between
classical mean-field trajectories for $\alpha=0$ (left column) and many-body
energy eigenstates for $\alpha\neq 0$ (center and right columns). Classical
integrability gives room to robust quantum many-body scars localized in the
large-spin sectors of Hilbert space, coexisting with a bulk of repelling
energy levels for $0<\alpha<d$ (top row). Classical chaos, instead, prevents
quantum many-body scarring for arbitrarily small $\alpha>0$ (bottom row). The
figure illustrates this dichotomy with data for the quantum Ising model, Eq.
(4) with $\gamma=1$. We report: classical collective-spin trajectories
$\langle\hat{\mathbb{S}}(t)\rangle/L$ for $\alpha=0$ and $L\to\infty$ (left);
exact energy spectra resolved in terms of the average collective spin size
$\langle|\hat{\mathbb{S}}|\rangle$ for $\alpha=0.5$ and $L=18$ (center);
weights on the unperturbed $\alpha=0$ basis $\\{|E^{(0)}_{n}\rangle\\}$ of the
representative eigenstate $|E^{*}\rangle$ highlighted in red in the center
column (right); with time-independent Hamiltonian (top) or subject to strong
periodic driving (bottom) — in the latter case we report stroboscopic
trajectories and Floquet spectra. Quantum many-body scarring is here
highlighted by wavefunction concentration in large collective-spin size
sectors (details are given further below in the paper).

  * •

Performing an extensive study of numerical energy spectra of variable-range
quantum Ising chains, we show that random-matrix-like level repulsion appears
for infinitesimal $\alpha$ along with resilient anomalous eigenstates,
characterized by large collective spin size, low entanglement, and high
overlap with product states (Sec. IV). We identify these eigenstates as
candidate QMBS-like states, smoothly deformed from the large-spin eigenstates
of the permutationally-symmetric (solvable) mean-field limit $\alpha=0$. Using
a refined numerical approach which involves a projection of the Hamiltonian
onto large-spin subspaces (App. B), we find evidence for robustness of the
QMBS for arbitrary $L$ when $0<\alpha\lesssim d$, and for their eventual
hybridization when $\alpha\gg d$.

  * •

We formulate a self-consistent analytical theory of eigenstate localization in
Hilbert space that fully backs our numerical observations (Sec. V). First, we
introduce a procedure to map an interacting spin lattice onto a collective
spin coupled to an ensemble of bosonic modes, physically representing magnon
excitations with non-vanishing momenta. The coupling strength is controlled by
$\alpha$ and $d$ and it is vanishingly small in the mean-field limit. In a
Hilbert space shell with fixed energy density and large collective spin, the
exact spin-boson description reduces to a simpler but still interacting
quantum impurity model: a relativistic quantum rotor parametrically coupled to
a quadratic bosonic bath. We find an exact solution to this effective
Hamiltonian, which allows us to compute the degree of eigenstate
delocalization in Hilbert space and hence to check self-consistency. We find
that the eigenstates of this rotor-magnon Hamiltonian are self-consistent QMBS
of the original spin Hamiltonian only for $0<\alpha<d$. Figure 1 visually
summarizes the steps of our eigenstate localization theory, which relies on
two fundamental ingredients:

    1. 1.

sufficiently slow interaction decay $0<\alpha<d$;

    2. 2.

classical integrability of the mean-field limit.

These QMBS thus possess, by construction, a clear fingerprint from
semiclassical periodic orbit – sometimes considered a defining feature of
QMBS.

  * •

We predict and classify instabilities of the QMBS, arising from violations of
the two conditions above. Eigenstate delocalization most notably occurs in
presence of discrete-symmetry breaking or of mean-field chaotic dynamics. We
accurately test these predictions for the long-range quantum Ising chain,
respectively in the ordered phase or in presence of periodic kicks (Sec. VI).
Figure 2 reports representative examples.

The results of this work have a number of broader consequences, which we
briefly elaborate on in the conclusive Sec. VII.

##  III The model

In this paper we consider a $d$-dimensional lattice of $V=L^{d}$ quantum spins
governed by the Hamiltonian

$\hat{H}_{\alpha}=-\sum_{\mathbf{i},\mathbf{j}}J_{\mathbf{i},\mathbf{j}}\bigg{(}\frac{1+\gamma}{2}\hat{\sigma}^{x}_{\mathbf{i}}\hat{\sigma}^{x}_{\mathbf{j}}+\frac{1-\gamma}{2}\hat{\sigma}^{y}_{\mathbf{i}}\hat{\sigma}^{y}_{\mathbf{j}}+\frac{\Delta}{2}\hat{\sigma}^{z}_{\mathbf{i}}\hat{\sigma}^{z}_{\mathbf{j}}\bigg{)}-h\sum_{\mathbf{j}}\hat{\sigma}^{z}_{\mathbf{j}}$
(1)

where $\gamma$ and $\Delta$ parametrize the XY and XXZ anisotropies,
respectively. In this equation
$\hat{\sigma}^{x,y,z}_{\mathbf{r}}=\hat{s}^{x,y,z}_{\mathbf{j}}/s$ are
rescaled spin-$s$ operators acting on the spin at site
$\mathbf{j}=(j_{1},\dots,j_{d})$ of the lattice, where $j_{a}=1,\dots,L$ for
$a=1,\dots,d$. For spins-$1/2$ these are the standard Pauli matrices.

We will consider interactions depending only on the distance between spins,
with algebraic decay:

$J_{\mathbf{i},\mathbf{j}}\equiv
J_{\mathbf{i}-\mathbf{j}}=\frac{J}{\lvert\lvert\mathbf{i}-\mathbf{j}\rvert\rvert^{\alpha}}.$
(2)

where
$\lvert\lvert\mathbf{i}-\mathbf{j}\rvert\rvert\equiv\sqrt{\sum_{a=1}^{d}[\mathrm{Min}(|i_{a}-j_{a}|,L-|i_{a}-j_{a}|)]^{2}}$
realizes periodic boundary conditions. The constant $J$ is chosen to make the
mean-field energy independent of the decay exponent $\alpha$ (Kac
normalization Kac and Thompson (1969)),

$J=\frac{J_{0}}{2\mathcal{N}_{\alpha,L}},\qquad\mathcal{N}_{\alpha,L}=\frac{1}{2}\sum_{\mathbf{r}\neq\mathbf{0}}\frac{1}{\lvert\lvert\mathbf{r}\rvert\rvert^{\alpha}}\,.$
(3)

Such prescription is necessary to have a well defined thermodynamic limit for
$\alpha\leq d$: The system-size divergent rescaling factor
$\mathcal{N}_{\alpha,L}\sim L^{d-\alpha}$ ensures that the energy per spin
remains finite.

For definiteness, we will assume ferromagnetic interactions $J_{0}>0$
throughout, and we set $s=1/2$ 111For $s=1/2$ the terms
$\mathbf{i}=\mathbf{j}$ produce an inconsequential additive constant
$E=\sum_{\mathbf{j}}J_{\mathbf{i},\mathbf{j}}(1+\Delta/2)$, as Pauli matrices
square to $1$; we will thus set $J_{\mathbf{i},\mathbf{i}}=0$.. Increasing the
range of interactions (i.e. decreasing $\alpha$) weakens spatial fluctuations,
leading the system toward its mean-field limit, analogously to increasing the
system dimensionality Dutta and Bhattacharjee (2001); Defenu  _et al._ (2017).
It is thus convenient to perform our theoretical analysis with fixed $d$ and
study the properties of the model upon varying $\alpha$. As the effects of
spatial fluctuations are strongest in one dimension, we will work with a
variable-range XY quantum spin chain

$\hat{H}_{\alpha}=-\frac{J_{0}}{\mathcal{N}_{\alpha,L}}\sum_{j=1}^{L}\sum_{r=1}^{L/2}\frac{1}{r^{\alpha}}\bigg{(}\frac{1+\gamma}{2}\hat{\sigma}^{x}_{j}\hat{\sigma}^{x}_{j+r}+\frac{1-\gamma}{2}\hat{\sigma}^{y}_{j}\hat{\sigma}^{y}_{j+r}\bigg{)}-h\sum_{j=1}^{L}\hat{\sigma}^{z}_{j}$
(4)

(periodic boundary conditions are understood). We anticipate that our results
do not rely at all on any of the model restrictions above, as we will
explicitly show in Sec. V.8.

For $\alpha=\infty$ the Hamiltonian (4) reduces to the XY quantum spin chain
with nearest-neighbor interactions. This model is exactly solvable by mapping
to a quadratic fermionic chain via Jordan-Wigner transformation Lieb  _et al._
(1961). For finite $\alpha$ integrability is broken by couplings beyond the
nearest neighbors, as the associated Jordan-Wigner string turns into multi-
fermion interactions. As $\alpha$ is decreased, longer-range interactions
cooperate to suppress spatial fluctuations, resulting in a qualitative
enhancement of the system’s ability to order in _excited_ states for
$\alpha<2$ Dyson (1969).

The tendency to sustain collective spin alignment becomes increasingly
prominent as $\alpha$ is decreased. This is conveniently understood by viewing
a long-range interacting system as a “perturbation” to the infinite-range
limit with all-to-all interactions ($\alpha=0$). For $\alpha\to 0$ Eq. (4)
reduces to a Hamiltonian describing dynamics of a single collective spin,
completely decoupled from all the other degrees of freedom — realizing the
mean-field description of the system. The effect of spatially modulated
interactions $\alpha\neq 0$ is then to dynamically couple this collective
degree of freedom to all the other finite-wavelength modes describing spatial
fluctuations, resulting in complex quantum many-body dynamics.

For later purpose we show how to make this picture explicit. We express the
finite-range Hamiltonian $\hat{H}_{\alpha}$ as a perturbation to
$\hat{H}_{\alpha=0}$ by rewriting it in momentum space. Let us define the
Fourier transform the spin operators (recall
$\hat{s}^{\mu}_{j}=\hat{\sigma}^{\mu}_{j}/2$ for $\mu=x,y,z$)

$\tilde{S}^{\mu}_{k}=\sum_{j=1}^{L}e^{ikj}\hat{s}^{\mu}_{j}\,,$ (5)

with

$k\equiv k_{\ell}=2\pi\ell/L,\quad\ell=0,\pm 1,\pm 2,\dots,\pm\lfloor
L/2\rfloor$ (6)

(for $L$ even $k_{\pm L/2}$ coincide);
$\tilde{\mathbb{S}}_{k=0}\equiv\hat{\mathbb{S}}=\sum_{j}\hat{\mathbb{s}}_{j}$
is the system’s collective spin. In terms of momentum-space operators the
variable-range XY spin chain reads

$\hat{H}_{\alpha}=-\frac{J_{0}}{L}\sum_{k}f_{k}(\alpha)\Big{[}\big{(}\tilde{S}^{+}_{k}\tilde{S}^{-}_{-k}+\tilde{S}^{-}_{k}\tilde{S}^{+}_{-k}\big{)}+\gamma\big{(}\tilde{S}^{+}_{k}\tilde{S}^{+}_{-k}+\tilde{S}^{-}_{k}\tilde{S}^{-}_{-k}\big{)}\Big{]}-2h\tilde{S}^{z}_{k=0}\,,$
(7)

where $\tilde{S}^{\pm}_{k}=\tilde{S}^{x}_{k}\pm i\tilde{S}^{y}_{k}$. (Note
that in this expression the various $k$-modes are _not_ dynamically decoupled,
as
$[\tilde{S}^{\alpha}_{k},\tilde{S}^{\beta}_{q}]=i\epsilon^{\alpha\beta\gamma}\tilde{S}^{\gamma}_{k+q}$.)
In Eq. (7) we defined the function

$f_{k}(\alpha)=\frac{1}{\mathcal{N}_{\alpha,L}}\sum_{r=1}^{L/2}\frac{\cos(kr)}{r^{\alpha}}=\frac{\sum_{r=1}^{L/2}\frac{\cos(kr)}{r^{\alpha}}}{\sum_{r=1}^{L/2}\frac{1}{r^{\alpha}}}.$
(8)

By construction, $f_{k=0}(\alpha)=1$.

Figure 3: Left panel: Couplings $f_{k}$ in Eq. (8) for a range of values of
$\alpha$; here we set $L=1000$. Right panel: convergence of the discrete
sequence $f_{\ell}(\alpha)$ for $\alpha=0.5$ as $L\to\infty$; the dashed line
highlights the asymptotics $\sim\ell^{-(1-\alpha)}$ from Eq. (11).

We now identify the variable-range perturbation by singling out the $k=0$ part
of the Hamiltonian:

$\hat{H}_{\alpha}=\hat{H}_{\alpha=0}(\hat{\mathbb{S}})+\hat{V}_{\alpha}\,,\qquad\text{with
}\;\hat{V}_{\alpha}=-\frac{J_{0}}{L}\sum_{k\neq
0}f_{k}(\alpha)\Big{[}\big{(}\tilde{S}^{+}_{k}\tilde{S}^{-}_{-k}+\tilde{S}^{-}_{k}\tilde{S}^{+}_{-k}\big{)}+\gamma\big{(}\tilde{S}^{+}_{k}\tilde{S}^{+}_{-k}+\tilde{S}^{-}_{k}\tilde{S}^{-}_{-k}\big{)}\Big{]}.$
(9)

Note that the term $\hat{V}_{\alpha}$ is in general an extensive perturbation
222The factor $1/L$ is due to the convention chosen for the spin Fourier
transform in Eq. (5)..

The quantity $f_{k}(\alpha)$ controls how strongly spin fluctuation modes with
various wavelengths come into play, thus determining the physical properties
of $\hat{H}_{\alpha}$. Its behavior upon varying $\alpha$ is summarized in
Fig. 3. Its shape shrinks from $f_{k}(\alpha\to\infty)=\cos k$ to

$f_{k}(\alpha\to 0)=\delta_{k,0}\,,$ (10)

becoming increasingly singular at $k=0$ as $\alpha$ is decreased:

$\left\\{\begin{split}f_{k}(\alpha)&\sim 1-c(\alpha)k^{2}&\text{for
}&\alpha>3;\\\ f_{k}(\alpha)&\sim 1-c(\alpha)|k|^{\alpha-1}&\text{for
}&1<\alpha<3;\\\ f_{k_{\ell}}(\alpha)\equiv f_{\ell}(\alpha)&\sim
c(\alpha)|\ell|^{-(1-\alpha)}&\text{for }&\alpha<1.\\\ \end{split}\right.$
(11)

These and other properties of $f_{k}(\alpha)$ are derived in App. A. The last
case is the most interesting for this paper: For $\alpha<1$, $f_{k}(\alpha)$
squeezes onto the vertical axis as $L\to\infty$; upon zooming near $k=0$ one
finds that the values $f_{2\pi\ell/L}$ converge as $L\to\infty$ to _a discrete
sequence of finite values,_ denoted $f_{\ell}$, which fall off as
$|\ell|^{-(1-\alpha)}$, cf. the right panel of Fig. 3 Lerose and Pappalardi
(2020); Defenu (2021). Thus, for small $0<\alpha<1$, only modes with extensive
wavelengths $k_{\ell}\propto 1/L$ do non-trivially participate in dynamics and
interact; physically, in this regime interactions vary so slowly with the
spatial distance that the system behaves as a permutationally invariant system
over finite length scales, and it is thus unable to resolve finite
wavelengths. As $\alpha$ is increased further, all $k\neq 0$ get eventually
activated.

The splitting in Eq. (9) thus gives us an intuitive physical picture of
dynamics for $\alpha>0$ in terms of a collective spin weakly coupled to
finite-wavelength spin excitations. This representation can be expected to
provide a starting point for a meaningful perturbation theory when $f_{k\neq
0}(\alpha)$ is small. The precise notion of smallness will be clarified in
this paper.

##  IV Numerical analysis

We analyze _spectral_ properties of highly excited energy eigenstates

$\hat{H}_{\alpha}|E_{n}\rangle=E_{n}|E_{n}\rangle$ (12)

for finite $0<\alpha<\infty$. While dynamical relaxation to thermal
equilibrium is known to become slow as $\alpha$ is decreased (see
Introduction), it is generally believed that thermalization is ultimately
attained at long times for any $0<\alpha<\infty$, as the system is non-
integrable. Accordingly, the energy spectrum is believed to satisfy the
eigenstate thermalization hypothesis (ETH) Deutsch (1991); Srednicki (1994).
In this Section, we report extensive numerical simulations which show that the
actual scenario may be subtler. While we find level statistics compatible with
ETH down to very small values of $\alpha\gtrsim 0.001$ (in agreement with Ref.
Russomanno  _et al._ (2021), see also Ref. Fratus and Srednicki (2016)), we
will show that a relevant subset of eigenstates exhibits remarkable
deviations, characterized by anomalous QMBS-like properties such as large
collective spin size, low entanglement, and large overlap with product states.
Such deviations are increasingly pronounced as $\alpha$ is decreased. Our
numerical data suggest the permanence of such anomalous properties for
arbitrary $L$ at least for $\alpha\lesssim 1$.

For definiteness, in this Section, we set $\gamma=1$ in Eq. (4) (quantum Ising
chain). We further set $J_{0}=h=1$ so as to avoid signatures of the thermal
phase transition (or of the excited-state quantum phase transition Cejnar  _et
al._ (2021)) in highly excited states. In Sec. VI.2 below we will comment on
how these phenomena affect our results.

The Hamiltonian (4) has translation, spatial reflection, and spin-flip
symmetries, with associated quantum numbers $K=0,2\pi/L,\dots,2\pi(L-1)/L$,
$I=\pm 1$, and $P_{z}=\pm 1$ respectively. In our exact diagonalization (ED)
results, we fix for definiteness $K=0$, $I=1$, and $P_{z}=1$.

###  IV.1 Level repulsion

We study the energy eigenvalue statistics of the model by probing the level
spacings $s_{n}=E_{n}-E_{n-1}$ via the ratio between nearby energy gaps
Oganesyan and Huse (2007)

$r_{n}=\frac{\mathrm{Min}\\{s_{n},s_{n+1}\\}}{\mathrm{Max}\\{s_{n},s_{n+1}\\}}\
.$ (13)

Generic non-integrable models are expected to obey Wigner-Dyson level
statistics Bohigas  _et al._ (1984), characterized by a distribution of the
ratios $p_{\text{WD}}(r)=\frac{27}{4}\,r\,\frac{1+r}{(1+r+r^{2})^{5/2}}$ and
average level spacing ratio $\braket{r}_{\text{WD}}=0.5295$. On the other
hand, integrable Hamiltonians are characterized by Poisson statistics Berry
and Tabor (1977), with $p_{\text{P}}(r)=\frac{2}{(1+r)^{2}}$ and average
$\braket{r}_{\text{P}}=0.386$.

Figure 4: Level repulsion for the long-range Hamiltonian. The red line
indicates the Wigner–Dyson spectral statistics, while the blue is the Poisson
one. (Top) Distribution of the level spacing ratio $p(r)$ for $L=18$ for
$\alpha=0.01$ (left) and $\alpha=0.7$ (right). (Bottom) Average level spacing
ratio as a function of the system size $L$ for different
$\alpha=0.001,\,0.01,\,0.2,\,0.5,\,0.8,\,1.5,\,3$ and $h=1$.

Our data indicate that the mean level spacing ratio is compatible with the
Wigner–Dyson prediction for all $\alpha>0$. The numerical distribution of
$p(r)$ shows a very good agreement with the Wigner–Dyson prediction (dotted
red line) and clear disagreement from the Poisson one (blue line), as
displayed in Fig. 4ab for $\alpha=0.01$ (left) and $\alpha=0.7$ (right) for a
system of $L=18$ spins. This is well confirmed by the trend of $\braket{r}$
with system size. Figure 4c shows convergence to the Wigner–Dyson value
$\braket{r}_{\text{WD}}$ upon increasing system size $L$. This holds for
$\alpha$ as small as $0.001$, signalling high sensitivity of spectral
statistics to the interaction range. These results complement those of Ref.
Russomanno  _et al._ (2021).

Figure 5: QMBS in the long-range quantum Ising chain. We report numerical
energy spectra of the model in Eq. (4) with $L=18$, $\gamma=1$ and $h=1$,
restricted to the translation-invariant, parity and spin-flip symmetric
sector. Scatter plot of the observables $\langle\delta_{S}\rangle$ (row a),
half-chain entanglement entropy $S_{L/2}(E_{n})$ (row b), and overlap with the
polarized state $|\langle\uparrow|E_{n}\rangle|^{2}$ (row c) of each
eigenstate $|E_{n}\rangle$ as a function of the energy density $E_{n}/L$.
Columns (1-4) correspond to increasing values of $\alpha=0$, $0.2$, $0.5$,
$0.8$.

###  IV.2 Anomalous eigenstates

In spite of the level repulsion, we find that the spectrum of the long-range
quantum Ising Hamiltonian has a distinctive structure, featuring a
thermodynamically negligible fraction of atypical non-thermal eigenstates
running througout the many-body energy spectrum.

In Fig. 5 we plot, for all energy eigenstates $\\{\ket{E_{n}}\\}$:
a) the average depletion of the total spin
$\hat{\mathbf{S}}^{2}=|\hat{\mathbf{S}}|(|\hat{\mathbf{S}}|+1)$ from its
maximal value, i.e.

$\displaystyle\langle\hat{\delta}_{S}\rangle_{E_{n}}=\frac{L}{2}-\braket{}{\hat{\mathbf{S}}}{}_{E_{n}}\
,$ (14)

b) the half-chain entanglement entropy:

$S_{L/2}=\text{Tr}(\hat{\rho}_{L/2}\ln\hat{\rho}_{L/2}),\quad\hat{\rho}_{L/2}=\text{Tr}_{L/2+1,\dots,L}|E_{n}\rangle\langle
E_{n}|\ ,$ (15)

c) the overlap with the fully polarized state

$|\langle E_{n}|\Uparrow\rangle|^{2}\ ,$ (16)

where
$|\bm{\Uparrow}\,\rangle\equiv|\uparrow_{1}\uparrow_{2}\dots\uparrow_{L}\,\rangle$.
The anomalous structure of the many-body energy spectrum is highlighted by
either of these three quantities in Fig. 5.

In the limit $\alpha=0$ the Hamiltonian conserves the collective spin size
$[\hat{H}_{\alpha=0},\hat{\mathbf{S}}^{2}]=0$ and hence $\delta_{S}=L/2-S$ is
an exact quantum number. The energy spectrum forms unequally spaced towers
(_Peres lattice_ Peres (1984)) labelled by $\delta_{S}$ – the horizontal rows
in Fig. 5a1. The degeneracy of the eigenspaces increases exponentially with
$\delta_{S}$, i.e., as the collective spin size gets smaller. The so-called
_Dicke manifold_ Dicke (1954) spanned by the $L+1$ non-degenerate eigenstates
with maximum collective spin size $\delta_{S}=0$ is the permutationally
symmetric subspace. This tower of eigenstates corresponds to the large orange
dots in the figure. Their bipartite entanglement entropy is known to scale
logarithmically with system size, i.e., $S_{L/2}\sim\ln L$, away from the
edges of the spectrum Latorre  _et al._ (2005). The fully polarized state
$|\Uparrow\rangle$ has non-vanishing overlap with them only. Note that some
states are missing in the plot due to our symmetry restrictions 333Let us
recall that we diagonalize the Hamiltonian in the symmetry sector $I=1$
(symmetric under reflection about the middle of the spectrum) and $k=0$
(translation invariant). For this reason, there are only $L/2+1$ collective
states $\ket{N}$ rather than $L+1$, and we do not access states with
$\langle\delta_{S}\rangle=1$, which possess only non-zero momentum (cf. the
discussion in Sec. V.2 below)..

As $\alpha$ is increased from $0$ the degeneracies are lifted. While the
splitting of the exponentially degenerate sectors having
$\delta_{S}=\mathcal{O}(L)$ dominate the level spacing statistics and induces
chaotic features (see discussion above), scatter plots of the quantities
(14)-(16) in Fig. 5 indicate a strong deviation from the standard ETH scenario
for $\alpha>0$. Deviation is strongest for the eigenstates originating from
sectors with small $\delta_{S}\ll L/2$ (i.e. large collective spin), as these
towers of states are remarkably resilient to the finite range of interactions.
This is shown in columns 2–4 of Fig. 5 corresponding to $\alpha=0.2$, $0.5$,
$0.7$, respectively, obtained by full ED of a chain of $L=18$ spins. The
states with smallest $\braket{\hat{\delta}_{S}}$ (large orange dots) are still
characterized by entanglement entropy much lower than the rest of the
eigenstates and by very large overlap with the polarized initial state
$\ket{\Uparrow}$ (Fig. 5c) 444We note the presence of another set of
eigenstates with anomalous small entanglement entropy, see Fig. 5b2b3b4. These
possess a defined total spin, i.e. $\braket{\delta_{S}}=2n$ with
$n=1,2,\dotsc,L/4$ and can be identified with the top of each stripe.
Understanding these states goes beyond the scope of this work..

To better characterize the stability properties of these large-spin
eigenstates with small $\braket{\,\hat{\delta}_{S}\,}$ for larger systems, we
resorted to a refined numerical method drawing on Refs. Protopopov  _et al._
(2017, 2020). This consists in _restricting the Hilbert space to a subspace
$\mathcal{H}_{N_{\text{{shell}}}}$ with large collective spin size_, where
$N_{\rm shell}$ is the number of kept sectors $\delta_{S}=0,\dots,N_{\rm
shell}-1$. For instance, ${N_{\text{{shell}}}}=1$ corresponds to the Dicke
manifold $\delta_{S}=0$, while the full Hilbert space is retrieved for
${N_{\text{{shell}}}=L/2+1}$. Constructing an orthonormal basis in
$\mathcal{H}_{N_{\text{{shell}}}}$ carrying the quantum number $\delta_{S}$ is
not an elementary task, as eigenstates of $\hat{\mathbb{S}}^{2}$ are typically
highly entangled. We used the construction of Ref. Protopopov  _et al._ (2017,
2020) based on regular-Cayley-tree fusing rules. The crucial advantage of this
method is that the matrix elements of the basic spin operators
$\hat{\sigma}^{\mu}_{j}$ within $\mathcal{H}_{N_{\text{{shell}}}}$ can be
computed _analytically_. Thus, as the dimension of
$\mathcal{H}_{N_{\text{{shell}}}}$ scales only as
$L^{N_{\text{{shell}}}}/(N_{\text{{shell}}}-1)!$, for not too large
$N_{\text{{shell}}}$ this technique allowed us to construct a sparse matrix
representation of the truncated Hamiltonian
$\hat{H}_{\alpha}\rvert_{\mathcal{H}_{N_{\text{{shell}}}}}$ for system sizes
$L\gg 18$ much larger than reachable by full ED – see Fig. 16 in App. B for
precise estimates. The validity of this truncation is a posteriori justified
by studying the convergence of the various observables [listed as a), b), c)
above] upon increasing the cutoff ${N_{\text{shell}}}$. Here convergence was
achieved for $L>18$ when $\alpha\lesssim 1$, as reported in App. B. There we
also describe the details of our method, which may be of independent interest.

Our numerical data suggest that the large-spin eigenstates exhibit a _sub-
extensive growth of the spin depletion $\langle\hat{\delta}_{S}\rangle$_ and a
_logarithmic scaling of the half-chain entanglement entropy_ $S_{L/2}$ with
system size $L$ for $\alpha\lesssim 1$. In Fig. 6a, we plot
$\langle\hat{\delta}_{S}\rangle$ as a function of $L$ for various
$0<\alpha\leq 1$, restricting to the central [$(L/4+1)$-th in order of energy
in the considered symmetry sector] large-spin eigenstate, with closest energy
density $E_{n}/L\simeq 0$ to the infinite-temperature value. We compute these
quantities up to $L=40$. Data for increasing $N_{\text{{shell}}}$ are
represented by dots with decreasing size. Convergence is excellent even with
very small $N_{\text{{shell}}}$ for $\alpha\lesssim 0.5$, while it is slower
for larger $\alpha$. Results strongly suggest a sublinear scaling of
$\langle\hat{\delta}_{S}\rangle$ with $L$, consistent with saturation to a
finite value for small $\alpha$. Data for the entanglement entropy, shown in
Fig. 6b, are consistent with a logarithmic scaling with system size $L$. Note
that for each value of $L$ we used the largest $N_{\text{{shell}}}$ for which
observables are converged. This constrained us to smaller $L$ for larger
$\alpha$.

Figure 6:  Observables for the large-spin eigenstate closest to zero energy
for $\alpha=0.1,0.2,0.5,0.8,1$ in orange, red, purple, blue and green
respectively. (Top) Spin depletion $\langle\hat{\delta}_{S}\rangle$ (14) as a
function of system size $L$ for increasing $N_{\text{{shell}}}=1,3,4$
represented by dots with decreasing sizes. The dashed lines is given by the
theory prediction in Eq. (73) below. (Bottom) Half-chain entanglement entropy
scaling $S_{L/2}$ with $L$. Only data converged with $N_{\text{shell}}$ are
reported.

These results must be contrasted to the behavior found for $\alpha\gg 1$,
where large-spin eigenstates exhibit a tendency to delocalize through the
many-body spectrum. This is illustrated in Fig. 7. In panel (a) we report the
eigenstates’ spin depletion $\langle\hat{\delta}_{S}\rangle_{E_{n}}$ for
$\alpha=3$ and $L=16$. While largest-spin eigenstates can still be singled out
(large orange dots), the spin depletion is much larger than for $0<\alpha<1$,
cf. Fig. 5 row a. Analysis for various sizes $L$ suggests, furthermore, that
the separation of these eigenstates from the thermal bulk is only a finite-
size effect (Fig. 7b). The spin depletion
$\langle\hat{\delta}_{S}\rangle_{E_{n}}$ averaged over the $L/2+1$ eigenstates
with smallest value increases approximately linearly with $L$ for $\alpha>1$,
in contrast with the strong suppression found for $\alpha<1$.

In summary, we have shown that the numerical energy spectrum of the long-range
quantum Ising chain exhibits non-thermal features compatible with a QMBS
scenario. The system has a vanishing fraction of eigenstates characterized by
a) the depletion from maximal collective spin size
$\langle\hat{\delta}_{S}\rangle=L/2-\langle S\rangle$ scaling sub-extensively
with $L$, b) logarithmic scaling of the half-chain entanglement entropy
$S_{L/2}\sim\ln L$ [cf. Fig. 6] and c) large overlap with fully polarized
states [cf. Fig. 5]. Our numerical data suggest that these non-thermal
features may persist for arbitrarily large system size when $\alpha\lesssim
1$. Below we will back this numerical evidence with an analytical construction
of these _collective QMBS_ , which confirms the suggested stability and
delimits their robustness.

We note that in the limit $\gamma=0$ of Eq. (4) these QMBS reduce to the well-
known Anderson tower of states Anderson (1952); Tasaki (2019). In this limit,
however, their stability is a priori guaranteed by the $U(1)$ symmetry, as
they are the ground states in each symmetry sector.

Figure 7:  Eigenstates’ spin depletion $\langle\hat{\delta}_{S}\rangle$ in the
quantum Ising chain for shorter interaction range $\alpha>1$. (a) Scatter plot
of the exact energy spectrum for $\alpha=3$ with $L=16$ and $h=J_{0}=1$. (b)
Average spin depletion $\langle\delta_{S}\rangle_{E_{n}}$ of the $L/2+1$
eigenstates with smallest value [cf. large orange dots in panel (a)] as a
function of the system size $L$ for increasing $\alpha=0.5,1.5,3$.

###  IV.3 Initial state dependence in dynamics

The dynamic counterpart of the spectral structure of QMBS found above is given
by out-of-equilibrium dynamics depending strongly on the initial condition.
For completeness we report evidence of this behavior, exemplified by Fig. 8,
where we contrast initialization in a polarized state (blue) to a random
product state (red), subject to evolution governed by a long-range quantum
Ising Hamiltonian, (4) with $\gamma=1$ and $\alpha=0.7$.

Polarized states, such as $\ket{\Uparrow}=\ket{\uparrow\dots\uparrow}$, have a
large overlap with the QMBS [cf. Fig. 5]. Consequently, they display anomalous
dynamics of local observables such as the single-site spin polarization
$\braket{\hat{\sigma}_{1}^{z}(t)}$, see Fig. 8a, exhibiting strong revivals
over time scales that grow with $L$. This is accompanied by slow logarithmic
growth of the entanglement entropy, which eventually saturates to a value
scaling only logarithmically with the volume, see Fig. 8b. This phenomenology,
detected numerically in Refs. Schachenmayer  _et al._ (2013); Buyskikh  _et
al._ (2016), was explained via a semiclassical theory Lerose and Pappalardi
(2020) valid at long but finite times. (Qualitatively different behavior
occurs in proximity to semiclassical instabilities Lerose and Pappalardi
(2020); Mori (2019), as we discuss in Sec. VI.2.)

We contrast this behavior with dynamics initialized in a random product state
such as
$|\psi_{0}\rangle={|\uparrow_{1}\rangle}|\mathbf{s}_{2}\rangle\dots|\mathbf{s}_{L}\rangle$
where the first spin points along $z$ while other spins point in random
uncorrelated directions on the Bloch sphere,
$|\mathbf{s}_{i}\rangle=\cos(\theta_{i}/2)|\uparrow\rangle_{i}+e^{i\phi_{i}}\sin(\theta_{i}/2)|\downarrow\rangle_{i}$
with $\theta_{i}\in[0,\pi)$ and $\phi_{i}\in[0,2\pi)$. In this case, the
initial collective spin size is typically low and consequently, dynamics are
compatible with the standard thermalization scenario. Here the local spin
polarization attains the thermal value after a few oscillations, without
revivals; the entanglement entropy displays rapid (linear) growth in time and
volume-law saturation, see Fig. 8.

Figure 8: Dynamics of the long-range quantum Ising chain [Eq. (4) with
$\gamma=1$] with $\alpha=0.7$ from a fully polarized state $|\Uparrow\rangle$
(blue lines) and from a random product state – defined in the text – (red
lines) for system size $L=10,14,18$ in increasing color intensity. (a) Single-
site spin polarization dynamics
$\langle\psi(t)|\hat{\sigma}^{z}_{1}|\psi(t)\rangle$. (b) Entanglement entropy
growth $S_{L/2}(t)$.

##  V Quantum rotor-magnon theory

In this Section we analyze spectral properties at finite $\alpha$ using an
analytical approach we develop. This allows us to formulate a predictive
theory for the stability of non-thermal eigenstates observed in numerics.

Taking the exact solution of the $\alpha=0$ model as starting point (Secs.
V.1), we interpret the spectral degeneracy in terms of magnon excitations
(Sec. V.2). Hence, we provide basic intuition for the persistence of a few
non-thermal eigenstates for small finite $\alpha$ (Sec. V.3). Building on this
intuition we systematically map the original spin Hamiltonian to that of a
collective spin interacting with $L-1$ bosonic modes (Sec. V.4), and we
demonstrate that this mapping takes a particularly simple and physically
suggestive form in the large-spin sectors (Sec. V.5). Remarkably, we show that
the resulting effective model is an interacting impurity model, for which we
find an exact solution (Sec. V.6), and we compute the spreading of energy
eigenstates in Hilbert space as $\alpha$ is increased (Sec. V.7). These
results show that _non-thermal eigenstates are stable for arbitrary system
size $L$ when $0<\alpha<1$ in absence of semi-classical chaos_. We conclude
this Section by expounding the generality of our theory (Sec. V.8). We make
reference to Fig. 1 for a visual summary of the main steps of this Section.

We note that a related rotor-magnon description of the _low-energy_ spectrum
has been recently developed by Roscilde et al. Roscilde  _et al._ (2023a, b).

###  V.1 Spectrum and eigenstates in the infinite-range limit ($\alpha=0$)

We review the exact spectrum of the Hamiltonian (4) with infinite interaction
range ($\alpha=0$) Lipkin  _et al._ (1965); Bapst and Semerjian (2012); Mazza
and Fabrizio (2012); Sciolla and Biroli (2011); Ribeiro  _et al._ (2008),
which plays a central role in the following theory.

The Hamiltonian $\hat{H}_{\alpha=0}$ is invariant under all spin permutations.
It reduces to a function of the collective spin components
$\hat{S}^{x,y,z}=\frac{1}{2}\sum_{j=1}^{L}\hat{\sigma}^{x,y,z}_{j}$ :

$\hat{H}_{\alpha=0}=-\frac{4J_{0}}{L}\bigg{[}\frac{1+\gamma}{2}\big{(}\hat{S}^{x}\big{)}^{2}+\frac{1-\gamma}{2}\big{(}\hat{S}^{y}\big{)}^{2}\bigg{]}-2h\hat{S}^{z}\,.$
(17)

This Hamiltonian manifestly conserves the collective spin size
$\hat{\mathbb{S}}^{2}=(\hat{S}^{x})^{2}+(\hat{S}^{y})^{2}+(\hat{S}^{z})^{2}$.
The eigenvalues of $\hat{\mathbb{S}}^{2}$ can be written as $S(S+1)$, with

$S=L/2-\delta_{S},\qquad\delta_{S}=0,1,2,\dots,\bigg{\lfloor}\frac{L}{2}\bigg{\rfloor}.$
(18)

The good quantum number $\delta_{S}$ quantifies the discrepancy of the
collective spin size from its maximal value (cf. Fig.5a1).

Hilbert space sector associated with $\delta_{S}$ contains $d_{L,\delta_{S}}$
identical copies of a spin-$S$ representation of $SU(2)$, where

$d_{L,\delta_{S}}=\binom{L}{\delta_{S}}-\binom{L}{\delta_{S}-1}.$ (19)

In each such $(2S+1)$-dimensional space the Hamiltonian acts as Eq. (17)
thought as the Hamiltonian of a single spin of size $S$. We denote by

$E_{\delta_{S},N},\qquad N=0,1,\dots,L-2\delta_{S},$ (20)

the “tower” of $L-2\delta_{S}+1=2S+1$ single-spin energy eigenvalues in the
block $\delta_{S}$. We label the eigenstates as

$|\delta_{S},\kappa\,;N\rangle$ (21)

where the additional “silent” quantum number $\kappa$ parametrizes the
$d_{L,\delta_{S}}$-fold degeneracy. Thus, in this expression, the pair
$(\delta_{S},\kappa)$ identifies a $SU(2)$ tower and $N$ identifies the
eigenstate within that tower.

In all the towers with large collective spin $S=\rho L/2$ ($0<\rho\leq 1$) the
single-spin Hamiltonian approaches a semiclassical limit as $L\to\infty$. To
see this, observe that a general infinite-range spin-$1/2$ Hamiltonian can be
written as

$\hat{H}_{\alpha=0}=S\,\mathcal{H}(\hat{\mathbb{s}};\rho),\qquad\hat{\mathbb{s}}\equiv\bigg{(}\frac{\hat{S}^{x}}{S},\frac{\hat{S}^{y}}{S},\frac{\hat{S}^{z}}{S}\bigg{)}$
(22)

where $\mathcal{H}$ is a classical smooth function and $\rho$ enters as a
parameter. [For our model,
$\mathcal{H}(x,y,z;\rho)=-2J_{0}\rho\Big{(}\frac{1+\gamma}{2}x^{2}+\frac{1-\gamma}{2}y^{2}\Big{)}-2hz$.]
The Schrödinger equation reads

$i\frac{1}{S}\partial_{t}|\psi\rangle=\mathcal{H}(\hat{\mathbb{s}};\rho)|\psi\rangle,$
(23)

where $|\psi\rangle$ is a state in the $(2S+1)$-dimensional tower, and the
rescaled spin variables satisfy

$[\hat{s}^{\alpha},\hat{s}^{\beta}]=i\frac{1}{S}\epsilon^{\alpha\beta\gamma}\hat{s}^{\gamma}.$
(24)

Thus, the system has an effective $\hbar_{\rm{eff}}=1/S$, and the
thermodynamic limit coincides with the classical limit. The leading behavior
of spectrum, eigenstates, and dynamics as $L\to\infty$ can be computed by
standard semiclassical tools applied to $\mathcal{H}$, i.e. Bohr-Sommerfeld
quantization rule, Wenzel-Kramers-Brillouin (WKB) and truncated Wigner
approximation, respectively.

For the purposes of this paper, the following considerations suffice. The
single-spin classical Hamiltonian $\mathcal{H}$ is (trivially) Liouville-
integrable, and it can be recast to action-angle variables. For a given value
of the spin length $\rho$ we parametrize the classical phase space with
canonical variables $q,p$ [for example $\mathbb{s}=\;\big{(}\sqrt{1-p^{2}}\cos
q,\sqrt{1-p^{2}}\sin q,p\big{)}$]. Hence, we can rewrite

$\mathcal{H}(\mathbb{s};\rho)\mapsto\mathcal{H}(n;\rho)$ (25)

where $n=\frac{1}{2\pi}\oint_{E}p\,dq$ is the classical action variable,
corresponding to the phase-space area enclosed by a trajectory. Parametrizing
trajectories by their energy $E$, the representation (25) can be found by
inverting the equation $n=n(E;\rho)$. By construction, $\mathcal{H}$ does not
depend on the angle $\varphi$ canonically conjugated to $n$.

Figure 9: Sketch of the unperturbed spectrum of $\hat{H}_{\alpha=0}$, cf. Eq.
(26). Towers of energy levels associated with distinct $SU(2)$ irreps appear
as rows; level thickness graphically denotes the degeneracy
$d_{L,\delta_{S}}$. The inset highlights the locally regular 2d lattice at a
given energy density and $\delta_{S}\ll L/2$, characterized by the two
frequencies $\omega$ and $\bar{\omega}$, cf. Eq. (27).

According to the Bohr-Sommerfeld quantization rule the quantum energy levels
can be found by quantizing the classical action in integer multiples of
$\hbar_{\rm eff}$. We thus set $n=N/S=N/(\rho L/2)$, where $N=0,1,2,\dots,\rho
L$ is the integer quantum number labelling energy levels within the considered
tower [cf. Eq. (20)]. The dimensionless quantum operator $\hat{N}$ associated
with the classical action is simply defined by $\hat{N}|N\rangle=N|N\rangle$,
and the quantum operators $e^{\pm i\hat{\varphi}}$ associated with the
classical conjugated angle act as raising and lowering operators on the ladder
of eigenstates. On the other hand, the collective spin size $\rho$ is also
similarly quantized, $\rho=S/(L/2)=1-2\delta_{S}/L$. Taking $\delta_{S}\ll
L/2$ we can consistently neglect $\mathcal{O}(1/L^{2})$ corrections in $n$ and
obtain a formula for the semiclassical spectrum to leading order in
$\hbar_{\rm{eff}}=2/L$:

$E_{\delta_{S},N}\sim\frac{L}{2}\,\mathcal{H}\bigg{(}\frac{2N}{L};1-\frac{2\delta_{S}}{L}\bigg{)}.$
(26)

This is illustrated in Fig. 9.

The level spacing $E_{\delta_{S},N}-E_{\delta_{S},{N-1}}$ within a tower has a
semiclassical interpretation as $\hbar\omega(n;\rho)$, where
$\omega=\partial_{n}\mathcal{H}$ is the frequency of the classical orbit at
the corresponding energy. This quantity varies smoothly as $N$ runs through a
tower (at least away from isolated phase-space separatrices where $\omega=0$,
corresponding to singularities of the action-angle representation of
$\mathcal{H}$). To make the local structure of the semiclassical spectrum
manifest, we focus on a given energy shell and set $N=n\,L/2+\delta_{N}$, with
$\delta_{N}$ integer. Expanding in $1/L$, we obtain

$\begin{split}E_{\delta_{S},\delta_{N}}\;\sim\;&\frac{L}{2}\,\mathcal{H}\bigg{(}n+\frac{2\delta_{N}}{L};1-\frac{2\delta_{S}}{L}\bigg{)}\\\
\sim\;&L\,\mathcal{E}(n)\;+\omega(n)\,\delta_{N}+\bar{\omega}(n)\,\delta_{S}\;+\mathcal{O}\bigg{(}\frac{1}{L}\bigg{)}\,.\end{split}$
(27)

Thus, the energy spectrum around the energy density
$\mathcal{E}(n)=\mathcal{H}(n;1)/2$ is built by combining integer multiples of
two (positive) fundamental frequencies, $\omega\equiv\partial_{n}\mathcal{H}$
and $\bar{\omega}\equiv-\partial_{\rho}\mathcal{H}$, which vary smoothly with
$n$. This is illustrated in the inset of Fig. 9.

For later convenience, we note that the $\alpha=0$ eigenstates can be
expressed in the language of _spin-coherent states_. Here we synthetically
report the necessary information on this formalism; a minimalistic self-
contained review is in App. C. Considering the $SU(2)$ tower identified by
$(\delta_{S},\kappa)$, we define the spin-coherent state with spherical angle
$\bm{\Omega}\equiv(\theta,\phi)$ as the state with maximal spin projection in
the direction
$\vec{n}(\bm{\Omega})=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$
within the tower:

$\vec{n}(\bm{\Omega})\cdot\hat{\vec{S}}\,|\delta_{S},\kappa;\bm{\Omega}\rangle=\bigg{(}\frac{L}{2}-\delta_{S}\bigg{)}|\delta_{S},\kappa;\bm{\Omega}\rangle.$
(28)

The freedom in the choice of relative phases of spin-coherent states can be
fixed by choosing a reference one — usually in the $\vec{z}$ direction
$(\theta=0)$, i.e. $|\delta_{S},\kappa;\bm{\Uparrow}\rangle$ — and a
particular rotation protocol — usually the rotation by $\theta$ in the plane
generated by $\vec{z}$ and $\vec{n}(\bm{\Omega})$: denoting
$\vec{r}(\bm{\Omega})\equiv\frac{\vec{z}\times\vec{n}(\bm{\Omega})}{|\vec{z}\times\vec{n}(\bm{\Omega})|}$,
one has

$|\delta_{S},\kappa;\bm{\Omega}\rangle=\hat{U}(\bm{\Omega})|\delta_{S},\kappa;\bm{\Uparrow}\rangle,\quad\hat{U}(\bm{\Omega})\equiv
e^{-i\theta\,\vec{r}(\bm{\Omega})\cdot\hat{\vec{S}}}\,.$ (29)

For every choice of $\bm{\Omega}$ we denote the ladder of eigenstates of the
corresponding collective spin projection as

$\vec{n}(\bm{\Omega})\cdot\hat{\vec{S}}\,|\delta_{S},\kappa;M_{\bm{\Omega}}\rangle=\bigg{(}\frac{L}{2}-\delta_{S}-M\bigg{)}|\delta_{S},\kappa;M_{\bm{\Omega}}\rangle$
(30)

for $M=0,1,\dots,L-2\delta_{S}$. On varying $\vec{n}(\bm{\Omega})$ on the unit
sphere, spin-coherent states span the full tower (“overcomplete basis”): one
can show that

$\begin{split}\mathbb{1}_{\delta_{S},\kappa}=&\sum_{M=0}^{L-2\delta_{S}}|\delta_{S},\kappa;M_{\bm{\Uparrow}}\rangle\langle\delta_{S},\kappa;M_{\bm{\Uparrow}}|\\\
=&\frac{L-2\delta_{S}+1}{4\pi}\int
d\bm{\Omega}\,|\delta_{S},\kappa;\bm{\Omega}\rangle\langle\delta_{S},\kappa;\bm{\Omega}|\end{split}$
(31)

where $\int
d\bm{\Omega}\equiv\int_{0}^{2\pi}d\phi\int_{0}^{\pi}d\theta\sin\theta$. Thus,
every state in the tower can be written as a linear combination of spin-
coherent states 555This representation is not unique.. In particular, for the
eigenstates of $H_{\alpha=0}$ we can write

$|\delta_{S},\kappa;N\rangle=\int
d\bm{\Omega}\,\psi_{\delta_{S},N}(\bm{\Omega})|\delta_{S},\kappa;\bm{\Omega}\rangle.$
(32)

Here we introduced the _coherent-state wave function_
$\psi_{\delta_{S},N}(\bm{\Omega})\propto\langle\delta_{S},\kappa;\bm{\Omega}|\delta_{S},\kappa;N\rangle$
of the unperturbed eigenstate, which is efficiently accessible to both
semiclassical and numerical computations (see App. C).

In this work we will consider the $\alpha=0$ spectrum and eigenbasis as known
input data, which can be computed either analytically (via a large-$L$
asymptotic expansion) or numerically (via exact diagonalization of the single-
spin problem for $S\leq L/2\approx 10^{5}$). We briefly mention in passing
that the considered family of models is equivalent in the limit $\alpha=0$ to
the Lipkin-Meshkov-Glick model of nuclear physics Lipkin  _et al._ (1965),
which is actually Bethe-ansatz solvable Ortiz  _et al._ (2005); however, this
solution is not practically useful for large $L$, and semiclassical/numerical
techniques give much easier access to the relevant information.

###  V.2 Bosonic labelling of $\alpha=0$ eigenstates

For $\alpha=0$ the many-body Hilbert space fragments into many single-spin
Hilbert spaces. Each such $SU(2)$ irreducible representation is associated
with a regular energy tower, identified by the collective spin size $S$ and by
an additional quantum number $\kappa$ [cf. Eq. (21)] that parametrizes the
degeneracy of the towers with equal $S$. As soon as interactions acquire a
non-trivial spatial structure, e.g. for $\alpha\neq 0$, different towers get
coupled and the degeneracy gets lifted. To understand how the perturbation
affects the eigenstates, we first need to make sense of $\kappa$ in terms of
additional physical degrees of freedom, which are “inert” for $\alpha=0$ and
get activated as $\alpha\neq 0$.

In this Subsection we show how _the degeneracy of $\hat{H}_{\alpha=0}$ can be
resolved in terms of bosonic excitations (physically interpreted as magnons or
quantized spin waves) with non-vanishing momenta._ This way, each $\alpha=0$
eigenstate can be identified by the collective spin state – specified by the
quantized classical action $N=0,1,\dots,L-2\delta_{S}$ – and by the magnon
content of the state – specified e.g. by $\delta_{S}$ non-vanishing momenta.
This allows us to replace the abstract labelling in Eq. (21) with

$|\\{k_{1},\dots,k_{\delta_{S}}\\};N\rangle,\qquad
k_{1},\dots,k_{\delta_{S}}\neq 0\,,$ (33)

or, equivalently, the corresponding string of bosonic occupation numbers
$\\{n_{k}\\}_{k\neq 0}$, with $\delta_{S}=\sum_{k\neq 0}n_{k}$. In the
following Subsections we will leverage this description of unperturbed
eigenstates to reformulate the finite-range interacting spin Hamiltonian
$\hat{H}_{\alpha\neq 0}$ in terms of a large spin coupled to bosons.

For simplicity of presentation we derive the bosonic labelling by working with
the simplest collective Hamiltonian $\hat{H}_{\alpha=0}\mapsto\hat{S}^{z}$. We
construct a complete common eigenbasis of $\hat{\mathbb{S}}^{2}$ and
$\hat{S}^{z}$, labelled by a magnon occupation string $\\{n_{k}\\}_{k\neq 0}$
and the $z$-magnetization depletion $M$ (where
$S=L/2-\delta_{S}=L/2-\sum_{k\neq 0}n_{k}$ and $S^{z}=S-M$). The eigenstates
of the desired $\hat{H}_{\alpha=0}$ are linear combinations of the eigenstates
of $\hat{S}^{z}$ within individual $SU(2)$ towers, i.e. of states with
different $M$ but fixed magnon content, as we describe at the end.

Figure 10:  Illustration of the construction of a bosonic labeling. Here,
eigenspaces of $\hat{\mathbb{S}}^{2}$ correspond to vertical towers, with
degeneracy increasing leftwards. Individual towers are internally spanned by
collective operators (here $\tilde{S}^{-}_{k=0}$). We label the towers with
$k\neq 0$ magnon excitations created on top of the spin-coherent vacuum. The
quasi-orthogonality of dilute magnon states (blue-shaded area) is then
transferred through the entire towers, to large-spin states that cannot be
easily described with standard Holstein-Primakoff bosons (red-shaded area). We
use this bosonic labeling to resolve the degenerate eigenstates of an
arbitrary $\hat{H}_{\alpha=0}$, and to express the matrix elements of an
arbitrary $\hat{H}_{\alpha\neq 0}$.

We start from the spin-coherent state
$|\bm{\Uparrow}\,\rangle\equiv|\uparrow_{1}\uparrow_{2}\dots\uparrow_{L}\,\rangle$
and act with lowering operators $S^{-}_{k}$ [cf. Eq. (5)], generalizing the
procedure leading to Clebsch-Gordan coefficients Landau and Lifsits (1981). To
follow the construction it is helpful to make reference to the illustration in
Fig. 10:

  * •

The parent state $|\bm{\Uparrow}\,\rangle$ has collective quantum numbers
$S=L/2$ (i.e. $\delta_{S}=0$) and $S^{z}=L/2$. By repeatedly applying the
collective lowering operator $\tilde{S}^{-}_{k=0}$ to it, we generate the
entire $\delta_{S}=0$ tower. We normalize such states and denote them
$|\emptyset;M_{\bm{\Uparrow}}\rangle\propto(\tilde{S}^{-}_{k=0})^{M}{|\emptyset;\bm{\Uparrow}\,\rangle}$
for $M=0,1,\dots,L$. (These permutationally-invariant states are known as
_Dicke states_ Dicke (1954).)

  * •

To construct the towers composing the $S=L/2-1$ subspace (i.e. $\delta_{S}=1$)
we must act on $|\bm{\Uparrow}\,\rangle$ with spin lowering operators to
generate $L-1$ parent states with maximal $S^{z}=L/2-1$ and orthogonal to the
$\delta_{S}=0$ state $|\emptyset;1_{\bm{\Uparrow}}\,\rangle$ previously
constructed. A natural choice for this purpose is to act with
$\tilde{S}^{-}_{k_{1}}$, $k_{1}\neq 0$. This directly yields orthogonal states
$|\\{k_{1}\\};\bm{\Uparrow}\rangle\propto\tilde{S}^{-}_{k_{1}}|\bm{\Uparrow}\,\rangle$
with $\delta_{S}=1$, which we properly normalize (no need to orthogonalize).
By repeatedly applying the collective lowering operator $\tilde{S}^{-}_{k=0}$
to each of them, we generate the entire $\delta_{S}=1$ towers. We normalize
such states and denote them
$|\\{k_{1}\\};M_{\bm{\Uparrow}}\rangle\propto(\tilde{S}^{-}_{k=0})^{M}|\\{k_{1}\\};\bm{\Uparrow}\,\rangle$
for $M=0,1,\dots,L-2$.

  * •

We can then continue along this way: To construct the towers composing the
$S=L/2-2$ subspace (i.e. $\delta_{S}=2$) we act with $\tilde{S}^{-}_{k_{2}}$
on $|\\{k_{1}\\};\bm{\Uparrow}\rangle$ (with $k_{1},k_{2}\neq 0$) and
orthogonalize with respect to the states with $\delta_{S}=0,1$,
$\delta_{S}+M_{\bm{\Uparrow}}=2$ previously constructed, to generate
$L(L-3)/2$ new parent states with maximal $S^{z}=L/2-2$. The states
$\tilde{S}^{-}_{k_{2}}\tilde{S}^{-}_{k_{1}}|\bm{\Uparrow}\,\rangle$ do not
directly result in an orthogonal basis of the new subspace, as they have a
small $\mathcal{O}(1/L)$ overlap between themselves and with $\delta_{S}=0,1$
towers. In this case, we define the normalized states
$|\\{k_{1},k_{2}\\};\bm{\Uparrow}\rangle\propto\hat{P}_{\delta_{S}=2}\tilde{S}^{-}_{k_{2}}\tilde{S}^{-}_{k_{1}}|{\bm{\Uparrow}}\,\rangle$,
where $k_{1},k_{2}\neq 0$ and $\hat{P}_{\delta_{S}}$ projects on the
corresponding eigenspace of $\hat{S}^{2}$ (i.e., removes the components along
the previously constructed states). These states span the full eigenspace with
$\delta_{S}=2$ and $S^{z}=L/2-2$. They can be used to generate the entire
$\delta_{S}=2$ towers by repeatedly applying the collective lowering operator
$\tilde{S}^{-}_{k=0}$ to each of them. We normalize such states and denote
them
$|\\{k_{1},k_{2}\\};M_{\bm{\Uparrow}}\rangle\propto(\tilde{S}^{-}_{k=0})^{M}|\\{k_{1},k_{2}\\};\bm{\Uparrow}\,\rangle$
for $M=0,1,\dots,L-4$.

  * •

This procedure can be iterated to break down the full multi-spin Hilbert space
into simultaneous eigenspaces of $\hat{S}^{2}$ and $\hat{S}^{z}$: For
arbitrary $\delta_{S}$ we define the normalized states
$|\\{k_{1},\dots,k_{\delta_{S}}\\};\bm{\Uparrow}\rangle\propto\hat{P}_{\delta_{S}}\tilde{S}^{-}_{k_{\delta_{S}}}\dots\tilde{S}^{-}_{k_{1}}|\bm{\Uparrow}\,\rangle$
where $k_{1},\dots,k_{\delta_{S}}\neq 0$. Such states span the full eigenspace
with fixed ${\delta_{S}}$ and maximal $S^{z}=L/2-{\delta_{S}}$. By repeatedly
applying the collective lowering operator $\tilde{S}^{-}_{k=0}$ to each of
them we generate the entire towers that compose the ${\delta_{S}}$ eigenspace.
We normalize such states and denote them

$\quad\;\;|\\{k_{1},\dots,k_{\delta_{S}}\\};M_{\bm{\Uparrow}}\rangle\\\
\propto(\tilde{S}^{-}_{k=0})^{M}|\\{k_{1},\dots,k_{\delta_{S}}\\};\bm{\Uparrow}\rangle$
(34)

for $M=0,1,\dots,L-2\delta_{S}$. (Note the slight redundancy of notation
$|\\{k_{1},\dots,k_{\delta_{S}}\\};0_{\bm{\Uparrow}}\rangle\equiv|\\{k_{1},\dots,k_{\delta_{S}}\\};\bm{\Uparrow}\rangle$.)

Using this procedure, we constructed states defined by a set of _non-
vanishing_ momenta $k_{1},\dots,k_{\delta_{S}}\neq 0$ and an integer
$M=0,1,\dots,L-2\delta_{S}$. These states are exact simultaneous eigenstates
of $\hat{\mathbb{S}}^{2}$ and $\hat{S}^{z}$, and they span the entire
$2^{L}$-dimensional spin space. This set of states is however massively
redundant, as there is no restriction on the set of nonzero momenta. Indeed,
for $\delta_{S}=\mathcal{O}(L)$, two such states will typically have a finite
overlap. The states defined by Eq. (34) are only really useful when
$\delta_{S}\ll L/2$. In this “dilute” regime one can check the quasi-
orthonormality property:

$\langle\\{k^{\prime}_{1},\dots,k^{\prime}_{\delta_{S}^{\prime}}\\};M^{\prime}_{\bm{\Uparrow}}\,|\,\\{k_{1},\dots,k_{\delta_{S}}\\};M_{\bm{\Uparrow}}\rangle=\delta_{\delta_{S},\delta_{S}^{\prime}}\delta_{M,M^{\prime}}\Big{[}\delta_{\\{k_{1},\dots,k_{\delta_{S}}\\},\\{k^{\prime}_{1},\dots,k^{\prime}_{\delta_{S}}\\}}+\mathcal{O}\big{(}\delta_{S}/L\big{)}\Big{]}.$
(35)

This property can be traced back to the correspondence between _nearly
polarized_ spin states with large $S^{z}$ and dilute bosonic states, realized
by the Holstein-Primakoff transformation. To see this, let us define the
mapping

$\hat{s}^{-}_{j}\mapsto\hat{b}^{\dagger}_{j}(1-\hat{b}^{\dagger}_{j}\hat{b}_{j}),\;\hat{s}^{+}_{j}\mapsto(1-\hat{b}^{\dagger}_{j}\hat{b}_{j})\hat{b}_{j},\;\hat{s}^{z}_{j}\mapsto\frac{1}{2}-\hat{b}^{\dagger}_{j}\hat{b}_{j},$
(36)

where $\hat{b}_{j},\hat{b}^{\dagger}_{j}$ are canonical bosonic annihilation
and creation operators. [Note that this is an equivalent simpler version of
the standard exact Holstein-Primakoff transformation for $s=1/2$, see e.g.
Ref. Vogl  _et al._ (2020).] The mapping (36) should be understood as an
_embedding_ of the two-dimensional Hilbert space of a spin-$1/2$ in the
infinite-dimensional Hilbert space of a bosonic mode, with
$|\uparrow\,\rangle$ and $|\downarrow\,\rangle$ mapped onto
$|0\rangle\equiv|\emptyset\rangle$ and $|1\rangle\equiv
b^{\dagger}|\emptyset\rangle$ respectively. It is useful to define bosonic
Fock states

$|\\{n_{k}\\}\rangle\equiv\prod_{k}\frac{\big{(}\tilde{b}^{\dagger}_{k}\big{)}^{n_{k}}}{\sqrt{n_{k}!}}|\emptyset\rangle\,,$
(37)

where $\tilde{b}^{\dagger}_{k}=(1/\sqrt{L})\sum_{j}e^{ikj}b^{\dagger}_{j}$
creates a boson with momentum $k$ and $\\{n_{k}\\}$ is a string of occupation
numbers. In App. D we show that, in the dilute regime, the _maximally
polarized_ (i.e. $M=0$) spin states
$|\\{k_{1},\dots,k_{\delta_{S}}\\};\bm{\Uparrow}\rangle$ defined above are
close to the corresponding bosonic Fock states with unoccupied $k=0$ mode:

$\langle\\{n_{k}\\}|\\{k_{1},\dots,k_{\delta_{S}}\\};\bm{\Uparrow}\rangle=\delta_{\\{n_{k}\\},\\{k_{1},\dots,k_{\delta_{S}}\\}}+\mathcal{O}\big{(}\delta_{S}/L\big{)}\,.$
(38)

(Here, with a little abuse of notation, the Kronecker delta on the right-hand
side forces the set of momenta counted with their multiplicities to be that
identified by the string of occupation numbers.) Since Fock states are
orthonormal, Eq. (38) immediately implies Eq. (35) for $M=M^{\prime}=0$.
Hence, crucially, the quasi-orthonormality of the parent maximally polarized
states is transferred down throughout the respective $SU(2)$ towers, to states
$|\\{k_{1},\dots,k_{\delta_{S}}\\};M_{\bm{\Uparrow}}\rangle$ with arbitrary
$M$.

As illustrated in Fig. 10, the construction above non-trivially extends the
standard Holstein-Primakoff mapping between spins and bosons. Indeed, when not
only $\delta_{S}\ll L/2$ but also $M\ll L$, the correspondence between spin
states and bosonic states is one-to-one:

$\langle\\{n_{k}\\}|\\{k_{1},\dots,k_{\delta_{S}}\\};M_{\bm{\Uparrow}}\rangle=\\\
\delta_{n_{k=0},M}\,\delta_{\\{n_{k}\\}_{k\neq
0},\\{k_{1},\dots,k_{\delta_{S}}\\}}+\mathcal{O}\big{(}(M+\delta_{S})/L\big{)}.$
(39)

In other words, when the spin state is close to fully polarized ($L/2-S^{z}\ll
L$), the quantum number $M$ can simply be identified with the occupation of
the $k=0$ bosonic mode. In this corner of the multi-spin Hilbert space (blue-
shaded area in Fig. 10), the correspondence between nearly fully polarized
spin states and dilute bosonic Fock states is complete.

Equation (35) shows however that spin states
$|\\{k_{1},\dots,k_{\delta_{S}}\\};M_{\bm{\Uparrow}}\rangle$ with
$\delta_{S}\ll L/2$ but _arbitrarily large_ $M$ can also be accurately
labelled by their magnon content $\\{n_{k}\\}_{k\neq 0}$ (red-shaded area in
Fig. 10). Loosely speaking, Eq. (35) allows us to factorize the Hilbert space
of $L$ spins as direct product of a large spin and a Fock space of $L-1$
bosonic modes, which is accurate in the subspace $\delta_{S}\ll L/2$ with
dilute bosons. It is tempting to picture the state
$|\\{k_{1},\dots,k_{\delta_{S}}\\};M_{\bm{\Uparrow}}\rangle$ as a macroscopic
“condensate” of $M$ collective magnon excitations with $k=0$ accompanied by
dilute magnons with non-vanishing momenta $k_{1},\dots,k_{\delta_{S}}\neq 0$,
as suggested by the expression in Eq. (34). One must keep in mind, however,
that for extensive $M=\mathcal{O}(L)$ the total population
$M+\delta_{S}=\mathcal{O}(L)$ of Holstein-Primakoff bosons is far from dilute,
making the naive replacement of $\tilde{S}_{k=0}^{-}$ (the actual creation
operator of collective spin excitations) by $\tilde{b}^{\dagger}_{k=0}$ (the
Holstein-Primakoff creation operator) highly inaccurate. The spin state
$|\\{k_{1},\dots,k_{\delta_{S}}\\};M_{\bm{\Uparrow}}\rangle$ actually has
exponentially-small-in-$M$ overlap with its candidate bosonic partner
$|n_{0}=M,\\{n_{k}\\}_{k\neq 0}\rangle$. [In fact, it is straightforward to
check that only an exponentially small fraction of the latter bosonic Fock
state belongs to the physical spin space!]

The set of non-vanishing momenta $\\{n_{k}\\}_{k\neq 0}$ determines
$\delta_{S}=\sum_{k\neq 0}n_{k}$ and the silent quantum number $\kappa$
appearing in Eq. (21), with accuracy $\mathcal{O}\big{(}\delta_{S}/L\big{)}$.
In the construction above, states within each $SU(2)$ tower are identified by
the eigenvalue of the collective operator $\hat{S}^{z}$. As anticipated above,
the role of $\hat{S}^{z}$ can however be played as well by the quantized
action operator $\hat{N}$ associated with an arbitrary collective Hamiltonian
$\hat{H}_{\alpha=0}$. Correspondingly, the role of $\tilde{S}_{k=0}^{-}$ is
played by $e^{i\hat{\varphi}}$. This way, states within each tower are
identified by $N=0,1,\dots,L-2\delta_{S}$,

$|\\{k_{1},\dots,k_{\delta_{S}}\\};N\rangle$ (40)

in place of the collective $z$-magnetization depletion
$M=0,1,\dots,L-2\delta_{S}$. This labelling is equivalent to Eq. (21) upon
identifying
$(\delta_{S},\kappa)\leftrightarrow\\{k_{1},\dots,k_{\delta_{S}}\\}\leftrightarrow\\{n_{k}\\}_{k\neq
0}$. This is the conclusion we aimed for.

###  V.3 Eigenstate localization for $0<\alpha\ll 1$

The fully connected Hamiltonian $\hat{H}_{\alpha=0}$ features an extensive
tower of “ground-state-like” wavefunctions — the $\delta_{S}=0$ states —
running throughout the many-body energy spectrum. From a thermodynamic
viewpoint, such states are anomalous highly-excited states with large overlap
with product states, low entanglement, and regular energy level spacing. Each
such state finds exponentially many other states of non-collective character —
i.e. with $\delta_{S}=\mathcal{O}(L)$ — at equal energy density. As $\alpha>0$
states with different collective spin length get dynamically connected, and
hybridization of the “few” anomalous states with small $\delta_{S}$ with the
“many” normal states with $\delta_{S}=\mathcal{O}(L)$ becomes a priori
favorable. As discussed in Sec. IV, the numerical results indicate a different
scenario, showing a weak or suppressed tendency towards hybridization.

We argue here — and demonstrate below via explicit analytical arguments — that
hybridization of non-thermal eigenstates is hindered for small $\alpha$ by a
combination of two crucial ingredients:

  * i)

the suppression of matrix elements $\propto f_{k}(\alpha)$, ascribed to the
long range of the interactions;

  * ii)

the absence of energy resonances in the subspaces $\delta_{S}\ll L/2$,
ascribed to the classical integrability of the mean-field limit.

Figure 11: Sketch of the intuitive argument for eigenstate localization in the
main text. Left panel: Transitions generated by $\hat{V}_{\alpha}$ are
quasilocal on the unperturbed spectrum graph. The shading of the arrows
denotes the amount of hybridization generated at the perturbative order for
which the corresponding transition comes into play. Right panel: For
$\alpha<1$, we demonstrate that an exact eigenstate originating from the
$\delta_{S}\ll L/2$ subspace stays localized in a local neighbourhood.

This mechanism is pictorially illustrated in Fig. 11. The graph therein shows
the energy levels of $\hat{H}_{\alpha=0}$ arranged in a smoothly varying 2d
lattice spanned by the integers quantum numbers $(\delta_{S},\delta_{N})$, cf.
Eq. (27), connected by shaded arrows representing matrix elements for
$\alpha\neq 0$. As the interaction operator flips at most two spins, each
matrix element connects states with $\delta_{S}$ differing by at most $2$ —
and we will interpret the change in $\delta_{S}$ as creation and annihilation
of magnons 666For a Hamiltonian with $p$-spin interactions, matrix elements
connect states with $\delta_{S}$ differing at most by $p$.. At the same time,
the dependency of matrix elements on the collective spin variables is rooted
in semiclassics, resulting in analytic functions of $n=N/L$ and $\varphi$
(which will be rigorously shown below). These remarks show that transitions
generated by the perturbation $\alpha\neq 0$ are _quasi-local_ on the
$(\delta_{S},\delta_{N})$-lattice. Furthermore, since the two fundamental
frequencies $\omega,\bar{\omega}$ that locally define the lattice smoothly
vary with $n$, for a generic choice of energy shell they are typically _non-
resonant_ (i.e. incommensurate). This condition prevents energy denominators
to become arbitrarily small at any given order in perturbation theory
Giorgilli (2022). The simultaneous suppression of matrix elements for small
$\alpha$ thus prevents the occurrence of energy resonances to all orders,
hindering the delocalization of perturbed eigenstates over the graph in Fig.
9.

Remarkably, the basic mechanism illustrated above finds a natural counterpart
in an emergent _exact solvability_ of the perturbed energy eigenstates
originating from the subspaces $\delta_{S}\ll L/2$. The rest of this Section
is devoted to this demonstration, which constitutes the main theoretical
result of this paper. Combined with the solid observation of level repulsion
at arbitrarily small $\alpha>0$ reported in Sec. IV, this result proves the
coexistence of robust QMBS and a quantum-chaotic bulk many-body spectrum in
long-range interacting systems, which is the main message of this paper.

###  V.4 Matrix elements of $\hat{H}_{\alpha}$

In Sec. V.1 we reviewed how $\hat{H}_{\alpha=0}$ dynamically activates a
single collective degree of freedom, and in Sec. V.2 we showed how the
remaining $L-1$ dynamically frozen degrees of freedom can be conveniently
expressed as bosonic modes — physically associated with magnon or spin-wave
excitations with non-vanishing momenta. It is natural to ask how the variable-
range perturbation $\hat{V}_{\alpha}$ fits into this framework. Considering
the unperturbed spectrum as in Fig. 11, transitions between different towers
generated by $\hat{V}_{\alpha}$ (shaded arrows) correspond to creation or
destruction of few magnons simultaneously with jumps of the collective spin
quantum number $N$. The goal of this and the next Subsection is to derive an
explicit representation of the Hamiltonian $\hat{H}_{\alpha}$ as a large spin
interacting with an ensemble of bosonic modes. We will infer this
representation by computing matrix elements.

Starting from the expression of $\hat{H}_{\alpha}$ in Eq. (7), we insert on
both sides resolutions of the identity with unperturbed eigenstates (21). The
infinite-range Hamiltonian is by construction diagonal,

$\hat{H}_{\alpha=0}=\sum_{\delta_{S},\kappa}\sum_{N}E_{\delta_{S},N}|\delta_{S},\kappa;N\rangle\langle\delta_{S},\kappa;N|,$
(41)

with $E_{\delta_{S},N}$ given by Eq. (26). The variable-range perturbation can
be written as $\hat{V}_{\alpha}=\frac{1}{L}\sum_{k\neq
0}f_{k}(\alpha)\hat{V}_{k}$, with

$\hat{V}_{k}=-{J_{0}}\big{[}(\tilde{S}^{+}_{k}\tilde{S}^{-}_{-k}+\tilde{S}^{-}_{k}\tilde{S}^{+}_{-k})+\gamma(\tilde{S}^{+}_{k}\tilde{S}^{+}_{-k}+\tilde{S}^{-}_{k}\tilde{S}^{-}_{-k})\big{]}\,.$
(42)

We need to compute the matrix elements

$\langle\delta_{S}^{\prime},\kappa^{\prime};N^{\prime}|\hat{V}_{k}|\delta_{S},\kappa;N\rangle.$
(43)

In principle we have an efficient recipe to apply $\hat{V}_{k}$ to maximally
polarized states of the form $|\delta_{S},\kappa;\bm{\Uparrow}\rangle$, by
expressing both in terms of standard Holstein-Primakoff bosons. This strategy
directly works for low-energy states, e.g. $|\delta_{S},\kappa;0\rangle$,
which are nearly-polarized in the direction of the minimum of the classical
Hamiltonian $\mathcal{H}$. However, the “condensate” state
$|\delta_{S},\kappa;N\rangle$ is constructed as
$e^{iN\hat{\varphi}}|\delta_{S},\kappa;0\rangle$, where $e^{i\hat{\varphi}}$
is a function of collective spin operators $\hat{S}^{x,y,z}$ (known from the
$\alpha=0$ solution). We would then need to commute $N$ times the operators
$\tilde{S}^{\pm}_{k}$ appearing in $\hat{V}_{\alpha}$ with
$e^{i\hat{\varphi}}$ using the commutation relations
$[\hat{S}^{\alpha},\tilde{S}^{\beta}_{k}]=i\epsilon^{\alpha\beta\gamma}\tilde{S}^{\gamma}_{k}$.
While a priori well defined, this procedure is impractical, as the expression
of $\hat{\varphi}$ is generally cumbersome. We will thus follow a different
route using the formalism of spin-coherent states, see Sec. V.1 and App. C.

The advantage of decomposing unperturbed eigenstates as linear superpositions
of spin-coherent states, is that we know how to apply $\hat{V}_{k}$ to the
latter: Upon rotating the spin operators into a reference frame aligned with
the coherent-state polarization, we can use the standard Holstein-Primakoff
mapping to bosonize spin fluctuations. In other words, the complication of
dealing with a macroscopically populated collective spin mode is transferred
to computing integrals involving coherent-state wave functions.

To see how this machinery works in practice, let us develop Eq. (43) in spin-
coherent states:

$\int d\bm{\Omega^{\prime}}\int
d\bm{\Omega}\,\psi^{*}_{\delta_{S}^{\prime},N^{\prime}}(\bm{\Omega^{\prime}})\psi_{\delta_{S},N}(\bm{\Omega})\langle\delta_{S}^{\prime},\kappa^{\prime};\bm{\Omega}^{\prime}|\hat{V}_{k}|\delta_{S},\kappa;\bm{\Omega}\rangle.$
(44)

As $\hat{V}_{k}$ flips at most two spins, the quantum numbers $S$ and $S^{z}$
can only change by at most $2$ units. Thus, its action on the maximally
polarized state $|\delta_{S},\kappa;\bm{\Omega}\rangle$ results in a linear
combination of states of the form
$|\delta_{S}^{{}^{\prime\prime}},\kappa^{{}^{\prime\prime}};M_{\bm{\Omega}}\rangle$
with $|\delta_{S}^{{}^{\prime\prime}}-\delta_{S}|\leq 2$ and
$M_{\bm{\Omega}}\leq 4$ 777Recall $\delta_{S}=L/2-S$ and $M=S-S^{z}$.. To take
the overlap with the bra state
$\langle\delta_{S}^{\prime},\kappa^{\prime};\bm{\Omega}^{\prime}|$, we need to
decompose the latter on bra eigenstates of the collective spin projection
along $\vec{n}(\bm{\Omega})$, i.e.

$\langle\delta_{S}^{\prime},\kappa^{\prime};\bm{\Omega}^{\prime}|=\sum_{M=0}^{4}g^{*}_{\delta_{S}^{\prime},M}(\bm{\Omega}^{\prime},\bm{\Omega})\;\langle\delta_{S}^{\prime},\kappa^{\prime};M_{\bm{\Omega}}|\
,$ (45)

where we already implemented the restriction $M\leq 4$ arising from applying
$\hat{V}_{k}$ to the ket. The coefficients
$g_{\delta_{S},M}(\bm{\Omega}^{\prime},\bm{\Omega})$ are given by the standard
single-spin coherent-state wave functions of spin-$S=L/2-\delta_{S}$ states
with projection $S-M$ (see App. C). Explicitly,

$g_{\delta_{S},M}(\bm{\Omega}^{\prime},\bm{\Omega})=\binom{L-2\delta_{S}}{M}^{1/2}k(\bm{\Omega}^{\prime},\bm{\Omega})^{M}\langle\delta_{S};\bm{\Omega}|\delta_{S};\bm{\Omega}^{\prime}\rangle;$
(46)

where we set
$k(\bm{\Omega}^{\prime},\bm{\Omega})=\tan\big{(}\alpha(\bm{\Omega}^{\prime},\bm{\Omega})/2\big{)}e^{i\beta(\bm{\Omega}^{\prime},\bm{\Omega})}$,
and $\alpha$ and $\beta$ are the polar and azimuthal angle of
$\vec{n}(\bm{\Omega}^{\prime})$ in the rotated frame adapted to
$\vec{n}(\bm{\Omega})$.

Now we are left with matrix elements of the form

$\langle\delta_{S}^{\prime},\kappa^{\prime};M_{\bm{\Omega}}|\hat{V}_{k}|\delta_{S},\kappa;\bm{\Omega}\rangle.$
(47)

To evaluate this expression it is convenient to use Eq. (29)

$\langle\delta_{S}^{\prime},\kappa^{\prime};M_{\bm{\Uparrow}}|\hat{U}^{\dagger}(\bm{\Omega})\hat{V}_{k}\hat{U}(\bm{\Omega})|\delta_{S},\kappa;\bm{\Uparrow}\rangle\,;$
(48)

and transfer the rotations to spin operators
$\hat{U}^{\dagger}(\bm{\Omega})\tilde{S}^{\pm}_{\pm
k}\hat{U}(\bm{\Omega})\equiv\tilde{S}^{\pm}_{\pm k}(\bm{\Omega})$ appearing in
$\hat{V}_{k}$,

$\tilde{S}^{\pm}_{k}(\bm{\Omega})=R_{\pm,+}(\bm{\Omega})\tilde{S}^{+}_{k}+R_{\pm,-}(\bm{\Omega})\tilde{S}^{-}_{k}+R_{\pm,z}(\bm{\Omega})\tilde{S}_{k}^{z}$
(49)

where

$\begin{split}R_{+,+}(\bm{\Omega})=&R_{-,-}^{*}(\bm{\Omega})=-\sin^{2}(\theta/2)e^{2i\phi},\\\
R_{+,-}(\bm{\Omega})=&R_{-,+}^{*}(\bm{\Omega})=+\cos^{2}(\theta/2),\\\
R_{+,z}(\bm{\Omega})=&R_{-,z}^{*}(\bm{\Omega})=+\sin(\theta)e^{i\phi}.\\\
\end{split}$ (50)

Plugging in these transformations we obtain

$\hat{U}^{\dagger}(\bm{\Omega})\hat{V}_{k}\hat{U}(\bm{\Omega})=\sum_{\mu,\nu=+,-,z}{J}^{\mu,\nu}(\bm{\Omega})\,\tilde{S}^{\mu}_{k}\tilde{S}^{\nu}_{-k}\,,$
(51)

where we defined (omitting $\bm{\Omega}$-dependency for simplicity)

${J}^{\mu,\nu}\equiv-
J_{0}\Big{[}\big{(}R_{+,\mu}R_{-,\nu}+R_{-,\mu}R_{+,\nu}\big{)}\\\
+\gamma\big{(}R_{+,\mu}R_{+,\nu}+R_{-,\mu}R_{-,\nu}\big{)}\Big{]}\,.$ (52)

Note the symmetry property ${J}^{\mu,\nu}\equiv{J}^{\nu,\mu}$.

Here comes the crucial point: We have now isolated the amplitudes of the
collective spin’s and the magnons’ transitions in the expression of the matrix
element. Putting everything together and performing the integrals over the
spherical angles we obtain

$\langle\delta_{S}^{\prime},\kappa^{\prime};N^{\prime}|\hat{V}_{\alpha}|\delta_{S},\kappa;N\rangle=\sum_{k\neq
0}f_{k}(\alpha)\sum_{\mu,\nu=+,-,z}\sum_{M=0}^{4}\;\binom{L-2\delta_{S}}{M}^{1/2}\;\mathcal{J}^{\mu,\nu}_{\delta_{S}^{\prime},N^{\prime};\delta_{S},N;(M)}\;\langle\delta_{S}^{\prime},\kappa^{\prime};M_{\bm{\Uparrow}}|\frac{\tilde{S}^{\mu}_{k}\tilde{S}^{\nu}_{-k}}{L}|\delta_{S},\kappa;\bm{\Uparrow}\rangle\
,$ (53)

where we defined the coefficients

$\mathcal{J}^{\mu,\nu}_{\delta_{S}^{\prime},N^{\prime};\delta_{S},N;(M)}\;\equiv\;\int
d\bm{\Omega^{\prime}}\int
d\bm{\Omega}\;\psi^{*}_{\delta_{S}^{\prime},N^{\prime}}(\bm{\Omega^{\prime}})\psi_{\delta_{S},N}(\bm{\Omega})\,k(\bm{\Omega}^{\prime},\bm{\Omega})^{M}\,{J}^{\mu,\nu}(\bm{\Omega})\,\langle\delta_{S};\bm{\Omega}|\delta_{S};\bm{\Omega}^{\prime}\rangle\
.$ (54)

This is an _exact_ expression for the matrix elements of the finite-range
perturbation $\hat{V}_{\alpha}$ between _arbitrary_ eigenstates of an
infinite-range Hamiltonian $\hat{H}_{\alpha=0}$ (no relation between the two
has to be assumed). The computation of the matrix element is reduced to a
finite sum of matrix elements of two-spin operators between maximally (ket)
and a nearly-maximally (bra) polarized states in the corresponding $SU(2)$
towers, with coefficients computed from the single-spin $\alpha=0$ solution,
i.e. by taking integrals of simple trigonometric polynomials and the two
coherent-state wave functions.

###  V.5 Interacting rotor-magnon effective Hamiltonian

The expression (53) is only really useful when $\delta_{S}\ll L/2$. In this
case the remaining matrix element can be efficiently accessed numerically by
constructing the exact few-magnon subspace. Note that this procedure gives an
alternative efficient projection method, based on resolving the unperturbed
spectrum in terms of magnon occupations rather than the decorated trees
exploited above in Sec. IV and reviewed in App. B. A numerical comparison of
the performances of the two methods is left to future work.

In this Subsection we show that the problem in this form is amenable to
remarkable analytical progress. The key observation is that the remaining
matrix element states can be computed via the bosonic representation _in the
dilute approximation_ [in stark contrast with the original matrix element in
Eq. (43)!]. Concerning the operators, we take the Fourier transformation of
the exact Holstein-Primakoff mapping (36),

$\begin{split}\tilde{S}^{-}_{k}&\mapsto
L^{1/2}\;\tilde{b}^{\dagger}_{k}-\frac{1}{L^{1/2}}\sum_{q_{1},q_{2}}\tilde{b}^{\dagger}_{q_{1}}\tilde{b}^{\dagger}_{q_{2}}\tilde{b}_{q_{1}+q_{2}-k},\\\
\tilde{S}^{+}_{k}&\mapsto
L^{1/2}\;\tilde{b}_{-k}-\frac{1}{L^{1/2}}\sum_{q_{1},q_{2}}\tilde{b}^{\dagger}_{q_{1}+q_{2}+k}\tilde{b}_{q_{1}}\tilde{b}_{q_{2}},\\\
\tilde{S}^{z}_{k}&\mapsto\frac{L}{2}\delta_{k,0}-\sum_{q}\tilde{b}^{\dagger}_{q+k}\tilde{b}_{q}\,.\end{split}$
(55)

Concerning the states, following Sec. V.2, for small $\delta_{S}/L$ we can
identify $(\delta_{S},\kappa)\mapsto\\{n_{k}\\}_{k\neq 0}$ and make the
following replacement in Eq. (53):

$|\delta_{S},\kappa;M_{\bm{\Uparrow}}\rangle\quad\mapsto\quad\frac{\big{(}b^{\dagger}_{k=0}\big{)}^{M}}{\sqrt{M!}}\prod_{k\neq
0}\frac{\big{(}\tilde{b}^{\dagger}_{k}\big{)}^{n_{k}}}{\sqrt{n_{k}!}}|\emptyset\rangle$
(56)

[Recall $0\leq M\leq 4\ll L$ here, cf. Eq. (39).] Thus, the leading
contributions to the matrix elements can be determined for each value of $M$
by straightforward substitution of Eqs. (55) and (56) into Eq. (53):

  * •

For $M=0$, the leading contributions in $1/L$ are given by the terms

$\begin{split}\frac{1}{L}\tilde{S}^{\mp}_{k}\tilde{S}^{\pm}_{-k}&\;\mapsto\;\tilde{b}^{\dagger}_{k}\tilde{b}_{k}\,,\;\tilde{b}_{-k}\tilde{b}^{\dagger}_{-k}\,;\\\
\frac{1}{L}\tilde{S}^{\pm}_{k}\tilde{S}^{\pm}_{-k}&\;\mapsto\;\tilde{b}_{-k}\tilde{b}_{k}\,,\;\tilde{b}^{\dagger}_{k}\tilde{b}^{\dagger}_{-k}\,.\end{split}$
(57)

  * •

For $M=1$ contributions of the same leading order are given by cubic terms
with $\tilde{b}_{0}^{\dagger}$ and a zero-momentum pair:

$\begin{split}\frac{1}{L}\tilde{S}^{z}_{k}\tilde{S}^{-}_{-k}&\;\mapsto\;-\frac{1}{L^{1/2}}\tilde{b}^{\dagger}_{0}\,\tilde{b}_{-k}\tilde{b}^{\dagger}_{-k}\,;\\\
\frac{1}{L}\tilde{S}^{-}_{k}\tilde{S}^{z}_{-k}&\;\mapsto\;-\frac{1}{L^{1/2}}\tilde{b}^{\dagger}_{0}\,\tilde{b}^{\dagger}_{k}\tilde{b}_{k}\,;\end{split}$
(58)

indeed, in this case the extra factor $L^{1/2}$ in Eq. (53) compensates the
$1/L^{1/2}$ in the right-hand sides.

  * •

Similarly, for $M=2,3,4$, contributions of the same leading order are given by
terms with $(\tilde{b}_{0}^{\dagger})^{M}$ and a zero-momentum pair. Here the
factors $L^{-M/2}$ on the right-hand sides are exactly compensated by the
extra factors $L^{M/2}$ in Eq. (53). As it turns out, only the following term
with $M=2$ contributes: 888Note that here we used hermiticity
$\langle\delta_{S}^{\prime},\kappa^{\prime};N^{\prime}|\hat{V}_{k}|\delta_{S},\kappa;N\rangle=\langle\delta_{S},\kappa;N|\hat{V}_{k}|\delta_{S}^{\prime},\kappa^{\prime};N^{\prime}\rangle^{*}$
to rule out off-diagonal terms of the form
$(\tilde{b}^{\dagger}_{0})^{M}\,\tilde{b}_{-k}\tilde{b}_{k}$ for $M>0$.

$\frac{1}{L}\tilde{S}^{-}_{k}\tilde{S}^{-}_{-k}\;\mapsto\;-\frac{1}{L}(\tilde{b}^{\dagger}_{0})^{2}\,\big{(}\tilde{b}^{\dagger}_{k}\tilde{b}_{k}+\tilde{b}_{-k}\tilde{b}^{\dagger}_{-k}\big{)}\,.$
(59)

We conclude that the transitions in magnon space are described by a
_quadratic_ bosonic Hamiltonian (up to corrections of higher order in the
magnon density $\delta_{S}/L$).

What about the simultaneous transitions in collective-spin space? Here the
crucial observation is that the coefficients
$\mathcal{J}^{\pm,\pm}_{\delta_{S}^{\prime},N^{\prime};\delta_{S},N;(M)}$ can
be identified with matrix elements of a certain polynomial collective operator
$P(\hat{\mathbb{S}}/S)$ between eigenstates of the $\alpha=0$ problem, as
derived in App. C. It follows that the matrices
$\mathcal{J}^{\mu,\nu}_{\delta_{S}^{\prime},N^{\prime};\delta_{S},N;(M)}$ are
independent of the spin size and only sensitive to
$\delta_{S}-\delta_{S}^{\prime}$ to leading order in $1/L$. Furthermore, these
matrices have characteristic semiclassical structure, which is best
appreciated upon expressing the spin components as analytic functions of the
action and angle variables of the unperturbed problem: The analytic dependency
on the angle $\varphi$ translates into an (at most) exponential decay of the
matrix elements in $N^{\prime}-N$ away from the diagonal, and the analytic
dependency on the action $n$ translates into a smooth variation with
$(N+N^{\prime})/L$ along the diagonal. (These standard general properties can
easily be checked numerically for an arbitrary polynomial collective operator
$P(\hat{\mathbb{S}}/S)$ and an arbitrary single-spin Hamiltonian
$\hat{H}_{\alpha=0}$ eigenbasis.)

Putting everything together, we have found a remarkable result: for
$\delta_{S}\ll L/2$, the matrix elements of the finite-range Hamiltonian
$\hat{H}_{\alpha}$ coincide with those of an effective Hamiltonian describing
a large spin coupled to bosons,

$\hat{H}_{\alpha,\text{eff}}=\sum_{N}\bigg{[}\frac{L}{2}\mathcal{H}\bigg{(}\frac{N}{L/2};1\bigg{)}-\bar{\omega}\bigg{(}\frac{N}{L/2};1\bigg{)}\sum_{k\neq
0}\tilde{b}_{k}^{\dagger}\tilde{b}_{k}\bigg{]}|N\rangle\langle N|\\\
+\frac{1}{2}\sum_{N,N^{\prime}}\;\sum_{k\neq
0}f_{k}(\alpha)\Bigg{[}\mathcal{J}_{N^{\prime}-N}\bigg{(}\frac{N+N^{\prime}}{L}\bigg{)}\,\Big{(}\tilde{b}^{\dagger}_{k}\tilde{b}_{k}+\tilde{b}_{-k}\tilde{b}^{\dagger}_{-k}\Big{)}+\mathcal{K}_{N^{\prime}-N}\bigg{(}\frac{N+N^{\prime}}{L}\bigg{)}\tilde{b}_{-k}\tilde{b}_{k}+\mathcal{K}^{*}_{N^{\prime}-N}\bigg{(}\frac{N+N^{\prime}}{L}\bigg{)}\tilde{b}^{\dagger}_{k}\tilde{b}^{\dagger}_{-k}\Bigg{]}\,|N^{\prime}\rangle\langle
N|,$ (60)

where $\mathcal{J}_{N^{\prime}-N}\big{(}\frac{N+N^{\prime}}{L}\big{)}\equiv
2\big{(}\mathcal{J}^{+,-}_{0,N^{\prime};0,N;(0)}-\mathcal{J}^{z,-}_{0,N^{\prime};0,N;(1)}-\mathcal{J}^{-,-}_{0,N^{\prime};0,N;(2)}\big{)}$
and $\mathcal{K}_{N^{\prime}-N}\big{(}\frac{N+N^{\prime}}{L}\big{)}\equiv
2\mathcal{J}^{+,+}_{0,N^{\prime};2,N;(0)}$. This equation is our first main
result: _We have rewritten our finite-range Hamiltonian as an effective model
of a semiclassical particle (the collective spin), described by states
$\\{|N\rangle\\}$, coupled to an ensemble of bosonic modes associated with
magnon excitations with non-vanishing momentum ($k\neq 0$). _

While the Hamiltonian (60) is quadratic in the bosons, the coupling is highly
non-linear as the coefficients of the quadratic bosonic Hamiltonian depend on
the configuration of the collective spin. Since the semiclassical particle
spans macroscopically different energies, we can simplify the expression by
focusing on a given energy shell, and rewrite the effective Hamiltonian
therein. Setting $N=nL/2+\delta_{N}$, $N^{\prime}=nL/2+\delta_{N}+r$, and
consistently neglecting $\mathcal{O}(1/L)$ terms, we obtain the following
quantum rotor-magnon Hamiltonian:

$\hat{H}_{\alpha,\text{eff}}(n)=L\mathcal{E}(n)+\omega(n)\sum_{\delta_{N}}\delta_{N}|\delta_{N}\rangle\langle\delta_{N}|+\bar{\omega}(n)\sum_{k\neq
0}\tilde{b}_{k}^{\dagger}\tilde{b}_{k}\\\
+\frac{1}{2}\sum_{\delta_{N},r}\;\sum_{k\neq
0}f_{k}(\alpha)\Bigg{[}\mathcal{J}_{r}(n)\,\Big{(}\tilde{b}^{\dagger}_{k}\tilde{b}_{k}+\tilde{b}_{-k}\tilde{b}^{\dagger}_{-k}\Big{)}+\mathcal{K}_{r}(n)\,\tilde{b}_{-k}\tilde{b}_{k}+\mathcal{K}^{*}_{r}(n)\,\tilde{b}^{\dagger}_{k}\tilde{b}^{\dagger}_{-k}\Bigg{]}\;|\delta_{N}+r\rangle\langle\delta_{N}|\,.$
(61)

Alternatively, introducing the collective action displacement
$\hat{\delta}_{N}=\sum_{\delta_{N}}\delta_{N}|\delta_{N}\rangle\langle\delta_{N}|$
and angle $e^{\pm i\hat{\varphi}}=\sum_{\delta_{N}}|\delta_{N}\pm
1\rangle\langle\delta_{N}|$ operators, we can rewrite the rotor-magnon
Hamiltonian as

$\hat{H}_{\alpha,\text{eff}}(n)=L\mathcal{E}(n)+\omega(n)\,\hat{\delta}_{N}+\bar{\omega}(n)\sum_{k\neq
0}\tilde{b}_{k}^{\dagger}\tilde{b}_{k}\\\ +\frac{1}{2}\sum_{k\neq
0}f_{k}(\alpha)\Big{[}\mathcal{J}(n,\hat{\varphi})\Big{(}\tilde{b}^{\dagger}_{k}\tilde{b}_{k}+\tilde{b}_{-k}\tilde{b}^{\dagger}_{-k}\Big{)}+\mathcal{K}(n,\hat{\varphi})\;\tilde{b}_{-k}\tilde{b}_{k}+\mathcal{K}^{*}(n,\hat{\varphi})\;\tilde{b}^{\dagger}_{k}\tilde{b}^{\dagger}_{-k}\Big{]}\,.$
(62)

where $\mathcal{J}(n,\varphi)=\sum_{r}\mathcal{J}_{r}(n)\,e^{ir\varphi}$ and
$\mathcal{K}(n,\varphi)=\sum_{r}\mathcal{K}_{r}(n)\,e^{ir\varphi}$ are
$2\pi$-periodic analytic functions of $\varphi$, smoothly depending on the
parameter $n$. This equation accomplishes our goal. In the first line we
recognize $\hat{H}_{\alpha=0}$, cf. Eq. (27), where $\delta_{S}$ is now
resolved in terms of total occupation of magnon modes. The second line
expresses interactions between magnons and collective spin.

A neat physical picture emerges:

  * •

The non-linear precession of the collective spin, with frequency $\omega$,
constitutes a _periodic drive_ for magnon modes with mode-dependent driving
amplitude $f_{k}(\alpha)$;

  * •

This drive induces a population of magnons, which in turn _back-reacts_ on the
collective spin precession.

The nature of these interacting dynamics is expected to crucially depend on
the presence of _resonances_ between the unperturbed magnon frequency
$\bar{\omega}$ and the rotor frequency $\omega$. In the next Subsection we
will present an analytical solution of Eq. (62) which substantiates this
physical intuition.

###  V.6 Exact solution of the effective Hamiltonian

The structure of Eq. (62) is superficially reminiscent of the Bogolubov theory
of a weakly interacting Bose gas Leggett (2006). However, such an analogy is
only strictly valid at the isotropic point $\gamma=0$ of our spin Hamiltonian,
where the total magnetization $\hat{S}^{z}$ is conserved. This constrains
$\hat{\delta}_{N}$ to $\sum_{k\neq 0}\tilde{b}^{\dagger}_{k}\tilde{b}_{k}$,
making the former variable redundant, just like the conservation of the total
number of particles in an interacting Bose gas enslaves the condensate
population to the condensate depletion. In this case, the condensate
fluctuations (associated with the collective phase $\hat{\varphi}$) can be
formally eliminated, and one is left with an effective quadratic theory for
$k\neq 0$ excitations (parametrically depending on the self-consistent
condensate fraction) Leggett (2006). Crucially, in this case, the stability of
the Bose-condensed long-range ordered ground state is _a priori guaranteed_ by
the $U(1)$ symmetry. In our spin system with $\gamma=0$, the eigenstates with
smallest $\delta_{S}$ are indeed the _ground states_ of the Hamiltonian in the
various symmetry sectors.

The main goal of this paper is, however, to understand the fate of the QMBS
when _no symmetry protects them_ — e.g. in the Ising model ($\gamma=1$), see
our numerical study in Sec. IV. There, the effective Hamiltonian (62) is _not_
equivalent to a standard interacting Bose gas and it incurs no further
simplification: the dependency of the second line on $\hat{\varphi}$ cannot be
trivially eliminated. The delocalization of large-spin eigenstates in Hilbert
space for $\alpha\neq 0$ is energetically possible and entropically favorable,
making their stability far from obvious.

The crucial realization is that _we can successfully diagonalize the effective
Hamiltonian_ (62). This will allow us to compute the properties of the
presumed eigenstates at finite $\alpha$ and to a posteriori determine the
regime of self-consistency of the approximation leading to Eq. (62). In the
text below we report a sketch of the diagonalization procedure, emphasizing
its crucial points; we relegate a detailed derivation to App. E. We remark
that, to the best of our knowledge, our exact solution of an interacting
Hamiltonian of the form (62) is novel and it may thus be of independent
interest.

The diagonalization procedure involves a sequence of two canonical
transformations. First, we introduce an ansatz $e^{i\hat{S}}$ with a generator
of the form

$\hat{S}=\frac{1}{2}\sum_{k\neq
0}\bigg{[}F^{(+)}_{k}(\hat{\varphi})\,\Big{(}\tilde{b}_{k}\tilde{b}_{-k}+\tilde{b}^{\dagger}_{-k}\tilde{b}^{\dagger}_{k}\Big{)}\\\
+iF^{(-)}_{k}(\hat{\varphi})\,\Big{(}\tilde{b}_{k}\tilde{b}_{-k}-\tilde{b}^{\dagger}_{-k}\tilde{b}^{\dagger}_{k}\Big{)}\bigg{]}\
.$ (63)

The ansatz (63) represents a kind of generalized Bogolubov transformation
where the role of the unknown “angles” is taken by unknown _functions_ of the
collective operator $\hat{\varphi}$. Our goal is to determine functions
${F}^{(\pm)}_{k}$ in such a way that the transformed coefficients of
$\tilde{b}_{k}\tilde{b}_{-k}$,
$\tilde{b}^{\dagger}_{-k}\tilde{b}^{\dagger}_{k}$ in the Hamiltonian get
cancelled. Parametrizing the hyperbolic rotation of modes $(k,-k)$ by a
rapidity
$\eta_{k}=\sqrt{\big{(}F^{(+)}_{k}\big{)}^{2}+\big{(}F^{(-)}_{k}\big{)}^{2}}$
and an angle $\xi_{k}$ with $\tan\xi_{k}=-F^{(+)}_{k}/F^{(-)}_{k}$,
diagonalization imposes the following pair of _ordinary differential
equations_ for $\eta_{k}(\varphi)$, $\xi_{k}(\varphi)$ (see App. E):

$\omega\partial_{\varphi}\big{(}\eta_{k}\sin\xi_{k}\big{)}=+\bar{\omega}\,\sinh\eta_{k}\cos\xi_{k}\\\
+f_{k}(\alpha)\Big{[}\mathcal{J}^{(0)}(\varphi)\sinh\eta_{k}\cos\xi_{k}+\mathcal{J}^{(+)}(\varphi)\big{(}\cosh\eta_{k}\cos^{2}\xi_{k}+\sin^{2}\xi_{k}\big{)}+2\mathcal{J}^{(-)}(\varphi)\sinh^{2}(\eta_{k}/2)\cos\xi_{k}\sin\xi_{k}\Big{]}\,,$
(64)
$\omega\partial_{\varphi}\big{(}\eta_{k}\cos\xi_{k}\big{)}=-\bar{\omega}\,\sinh\eta_{k}\sin\xi_{k}\\\
-f_{k}(\alpha)\Big{[}\mathcal{J}^{(0)}(\varphi)\sinh\eta_{k}\sin\xi_{k}+2\mathcal{J}^{(+)}(\varphi)\sinh^{2}(\eta_{k}/2)\cos\xi_{k}\sin\xi_{k}+\mathcal{J}^{(-)}(\varphi)\big{(}\cosh\eta_{k}\sin^{2}\xi_{k}+\cos^{2}\xi_{k}\big{)}\Big{]}\;,$
(65)

where $\mathcal{J}^{(0)}\equiv\mathcal{J}$,
$\mathcal{J}^{(+)}\equiv\mathrm{Re}\,\mathcal{K}$,
$\mathcal{J}^{(-)}\equiv\mathrm{Im}\,\mathcal{K}$. Note that we obtain
differential equations rather than standard equations because of the non-
linear rotor-magnon interactions.

To understand the existence conditions of analytic $2\pi$-periodic solutions
to Eqs. (64), (65) we first observe that for $\alpha\to 0$ the second lines
vanish, and hence $\eta_{k}\equiv 0$; linearization in $\eta_{k}\sim
f_{k}(\alpha)$ yields the simplified equations

$\begin{split}\omega\partial_{\varphi}\big{(}\eta_{k}\sin\xi_{k}\big{)}&=+\bar{\omega}\,\eta_{k}\cos\xi_{k}+f_{k}(\alpha)\mathcal{J}^{(+)}(\varphi)\,,\\\
\omega\partial_{\varphi}\big{(}\eta_{k}\cos\xi_{k}\big{)}&=-\bar{\omega}\,\eta_{k}\sin\xi_{k}-f_{k}(\alpha)\mathcal{J}^{(-)}(\varphi)\;,\end{split}$
(66)

which are equivalent (upon renaming $P\equiv\eta_{k}\sin\xi_{k}$,
$Q\equiv\eta_{k}\cos\xi_{k}$, and rescaling the “time” variable
$\varphi\equiv\omega t$) to a fictitious harmonic oscillator with natural
frequency $\bar{\omega}$ externally driven at frequency $\omega$. As is well
known, there always exists a unique periodic trajectory of this system with
the same frequency $\omega$ of the drive (i.e., a $2\pi$-periodic solution in
$\varphi$), provided $\bar{\omega}\neq r\omega$, for all integers $r$
appearing in the Fourier representation of $\mathcal{J}^{(\pm)}$. In
correspondence of such resonances $\\{\bar{\omega}=r\omega\\}$, all solutions
are unbounded. Remarkably, the existence and uniqueness of a solution with the
same periodicity as the drive persists under a general weak non-linearity
Gallavotti (2013). The main qualitative effect of the non-linearity is to
_thicken_ the resonances from discrete points to finite intervals
$\\{|\bar{\omega}-r\omega|\leq\delta_{r}\\}$ of width depending on the driving
amplitude, thus upper bounded as
$\delta_{r}\leq\delta_{0}\,e^{-\sigma|r|}\,\mathrm{Max}_{k\neq
0}|f_{k}(\alpha)|$ for some constant $\sigma>0$. This condition provides a
theoretical guarantee of the existence of a finite range of parameter values
for which the ansatz (63) successfully diagonalizes the magnon part of the
effective Hamiltonian. From a practical point of view, the non-linear
equations (64), (65) are known explicitly for a given model, and one can
numerically determine its range of solvability throughout the parameter space.

Diagonalization of the magnon part of the Hamiltonian brings us within a
stone’s throw from the full solution of our problem. The transformed
Hamiltonian reads

$e^{i\hat{S}}\,\hat{H}_{\alpha,\text{eff}}\,e^{-i\hat{S}}=L{\mathcal{E}}+\omega\hat{\delta}_{N}^{\prime}+\sum_{k\neq
0}\bigg{[}G_{k}(\hat{\varphi})\bigg{(}\tilde{\beta}^{\dagger}_{k}\tilde{\beta}_{k}+\frac{1}{2}\bigg{)}-\frac{\bar{\omega}}{2}\bigg{]},$
(67)

where $\hat{\delta}_{N}^{\prime}$ and
$\tilde{\beta}_{k},\tilde{\beta}^{\dagger}_{k}$ denote the dressed (i.e.
transformed) rotor and magnon operators, respectively, and

$G_{k}(\hat{\varphi})\equiv\cosh\eta_{k}(\hat{\varphi})\big{[}\bar{\omega}+f_{k}(\alpha)\mathcal{J}^{(0)}(\hat{\varphi})\big{]}\\\
+\sinh\eta_{k}(\hat{\varphi})\cos\xi_{k}(\hat{\varphi})\,f_{k}(\alpha)\mathcal{J}^{(+)}(\hat{\varphi})\\\
+\sinh\eta_{k}(\hat{\varphi})\sin\xi_{k}(\hat{\varphi})\,f_{k}(\alpha)\mathcal{J}^{(-)}(\hat{\varphi})\,.$
(68)

As the magnon part is now diagonal, we can diagonalize the rotor part
separately in each sector defined by a Fock occupation $\\{n_{k}\\}_{k\neq 0}$
of dressed magnon modes. In fact, it can be checked that a magnon-diagonal
canonical transformation $e^{i\hat{F}}$ ansatz, with

$\hat{F}=\sum_{k\neq
0}F^{(0)}_{k}(\hat{\varphi})\bigg{(}\tilde{\beta}^{\dagger}_{k}\tilde{\beta}_{k}+\frac{1}{2}\bigg{)}\,,$
(69)

successfully decouples the rotor from the magnons if the generator satisfies
$\omega\partial_{\varphi}F^{(0)}_{k}(\varphi)\equiv{G_{k}(\varphi)}$ for all
$k$ (see App. E for more details). In this case, we find that the sequential
application of $e^{i\hat{S}}$ and $e^{i\hat{F}}$ fully diagonalizes our
effective rotor-magnon Hamiltonian (62):

$e^{i\hat{F}}e^{i\hat{S}}\,\hat{H}_{\alpha,\text{eff}}\,e^{-i\hat{S}}e^{-i\hat{F}}=L\mathcal{E}(\alpha)+\omega\,\hat{\Delta}_{N}+\sum_{k\neq
0}\bar{\omega}_{k}(\alpha)\tilde{\beta}^{\dagger}_{k}\tilde{\beta}_{k}\,,$
(70)

with a reference energy
$\mathcal{E}(\alpha)=\mathcal{E}+\frac{1}{2}\sum_{k\neq
0}(\bar{\omega}_{k}(\alpha)-\bar{\omega})$ dressed by zero-point fluctuations
and a nontrivial “dispersion relation” for dressed magnons,

$\bar{\omega}_{k}(\alpha)\equiv\overline{G_{k}(\varphi)},$ (71)

where
$\overline{\bm{\cdot}(\varphi)}=\int_{0}^{2\pi}\frac{d\varphi}{2\pi}\bm{\cdot}(\varphi)$.
Here, $\bar{\omega}_{k}(\alpha)$ replaces the flat band $\bar{\omega}$ of the
limit $\alpha=0$, cf. Eq. (27); it may be interpreted as a spectrum of
quasiparticle excitations on top of a _highly excited_ eigenstate.

The many-body spectrum thus takes the form
$E_{\\{n_{k}\\},\delta_{N}}(\alpha)=L\mathcal{E}(\alpha)+\omega\,\delta_{N}+\sum_{k\neq
0}\bar{\omega}_{k}(\alpha)n_{k}$: the many-body eigenstates are labelled by
the integer quantum numbers $\delta_{N}$ and $\\{n_{k}\\}_{k\neq 0}$
continuously from $\alpha=0$ to $\alpha$ (away from resonances!). Importantly,
however, in the original unperturbed basis the eigenstates feature non-trivial
quantum entanglement between the rotor and the magnons, encoded in the inverse
canonical transformation $e^{-i\hat{S}}e^{-i\hat{F}}$, as is manifest in the
expressions (63) and (69).

###  V.7 Self-consistency of localized eigenstates for $0<\alpha<1$

In the previous Subsection we determined the exact eigenstates of the
effective rotor-magnon Hamiltonian (62) for arbitrary $\alpha$. An important
property of these eigenstates is the extent of delocalization — i.e. the
amount of bare rotor and bare magnon excitations they contain — as a function
of $\alpha$. These quantities can be evaluated explicitly and play a crucial
role to establish the regime of self-consistency of the theory, i.e.
consistency with the assumptions leading to the effective Hamiltonian (62): An
eigenstate of the effective Hamiltonian can represent a legitimate eigenstate
of the original long-range XY spin chain only if its wavefunction remains
localized in the subspace $\delta_{S}\ll L/2$, $|\delta_{N}|\ll L$ (associated
with the “vertical” and “horizontal” delocalizations in Fig. 11,
respectively).

Both these quantities can be computed within the exact solution above.
Considering for simplicity an eigenstate in the magnon vacuum tower
$\\{n_{k}\\}_{k\neq 0}=\emptyset$ (calculation is analogous for other
eigenstates), $\Big{\langle}\hat{\delta}_{S}\Big{\rangle}$ is readily
evaluated as (see App. E)

$\Big{\langle}\hat{\delta}_{S}\Big{\rangle}=\Bigg{\langle}\sum_{k\neq
0}\tilde{b}^{\dagger}_{k}\tilde{b}_{k}\Bigg{\rangle}=\sum_{k\neq
0}\overline{\sinh^{2}\eta_{k}(\varphi)}$ (72)

On the other hand, the amount of rotor fluctuations in an eigenstate can be
estimated by evaluating, e.g.,
$\Big{\langle}\hat{\delta}_{N}^{2}\Big{\rangle}$. Calculation yields an
explicit expression, reported in App. E.

The case $0<\alpha<1$ is particularly interesting, as the values
$f_{k_{\ell}}(\alpha)\equiv f_{\ell}(\alpha)$ are discrete and decay as
$\ell^{-(1-\alpha)}$ (see App. A). Hence, for sufficiently large $\ell$ we are
always in the linear regime where $\eta_{\ell}\propto f_{\ell}(\alpha)$. This
directly leads to the estimates:

$\Big{\langle}\hat{\delta}_{S}\Big{\rangle}\sim\sum_{\ell=-L/2}^{L/2}|f_{\ell}(\alpha)|^{2}\sim\left\\{\begin{split}\mathcal{O}(1)&\quad\text{for
}0<\alpha<1/2\,,\\\ \log L&\quad\text{for }\alpha=1/2\,,\\\
L^{2\alpha-1}&\quad\text{for }1/2<\alpha<1\,.\\\ \end{split}\right.$ (73)

In this regime, one finds the same qualitative behavior for
$\sqrt{\Big{\langle}\hat{\delta}_{N}^{2}\Big{\rangle}}$ (see App. E). This
prediction for $\langle\hat{\delta}_{S}\rangle$ has been shown in Fig.6 above
and Fig.15 below, along with the numerical data. The plot shows complete
compatibility even for finite system size. The joint occurrence

$\Big{\langle}\hat{\delta}_{S}\Big{\rangle}\ll
L/2\,,\quad\sqrt{\Big{\langle}\hat{\delta}_{N}^{2}\Big{\rangle}}\ll
L\,,\quad\text{for }0<\alpha<1$ (74)

graphically illustrated in the right panel of Fig. 11, constitutes a _solid
analytical argument for our description of QMBS being self-consistent for
$0<\alpha<1$._

On the other hand, for $\alpha\gtrsim 1$, the delocalization of eigenstates
extends to sectors with extensively depleted collective spin
$\delta_{S}=\mathcal{O}(L)$, as
$\Big{\langle}\hat{\delta}_{S}\Big{\rangle}\sim
L\int_{-\pi}^{\pi}\frac{dk}{2\pi}|f_{k}(\alpha)|^{2}$. The small but finite
density of magnons breaks the eigenstate self-consistency, and full
hybridization of the QMBS with the thermal bulk may be expected for large $L$.
The dynamical counterpart of this expectation is that non-linear magnon-magnon
interactions come into play and trigger the eventual thermalization of initial
states with large collective spin. As such processes are suppressed with the
magnon density, thermalization is (at least) parametrically slow as $\alpha$
is decreased toward $1$, in agreement with previous numerical observations
Žunkovič  _et al._ (2018); Halimeh  _et al._ (2017); Liu  _et al._ (2019). We
stress, however, that we are not in a position to rule out other possible
subtler effects which may further delay (or even suppress) thermalization for
$\alpha\gtrsim 1$.

###  V.8 Generality of the theory

Here we comment on the generalization of the results of this Section, covering
the full range of parameters of model (1).

####  V.8.1 Higher dimension

Our stability result becomes stronger for higher dimensional spin lattices.
The only change we need is from $f_{k}(\alpha)$ to the analogous function
$f_{\mathbf{k}}(\alpha)$ defined in the $d$-dimensional Brillouin zone.
Throughout our derivation, we just replace $\sum_{k}f_{k}(\alpha)\dots$ with
$\sum_{\mathbf{k}}f_{\mathbf{k}}(\alpha)\dots$ everywhere. As shown in App. A,
the properties of $f_{k}(\alpha)$ for $\alpha<1$ directly generalize to
$f_{\mathbf{k}}(\alpha)$ for $\alpha<d$; in particular, $f$ squeezes toward
$\mathbf{k}=\mathbf{0}$ as $L\to\infty$, and becomes a function of the
discrete wavenumbers $\mathbf{k}L=2\pi\bm{\ell}$, decaying algebraically as
$|\bm{\ell}|^{-(d-\alpha)}$. By analogy we denote
$f_{\mathbf{k}_{\bm{\ell}}}\equiv f_{\bm{\ell}}$.

It follows that the collective spin depletion $\langle\hat{\delta}_{S}\rangle$
of solvable eigenstates is

$\Big{\langle}\hat{\delta}_{S}\Big{\rangle}\sim\sum_{\bm{\ell}}|f_{\bm{\ell}}(\alpha)|^{2}\sim\left\\{\begin{split}\mathcal{O}(1)&\quad\text{for
}0<\alpha<d/2\,,\\\ \log L&\quad\text{for }\alpha=d/2\,,\\\
L^{2\alpha-d}&\quad\text{for }d/2<\alpha<d\\\ \end{split}\right.$ (75)

Similar estimates apply to
$\sqrt{\Big{\langle}\hat{\delta}_{N}^{2}\Big{\rangle}}$. This shows self-
consistency of QMBS for $0<\alpha<d$.

####  V.8.2 General spin-spin interactions

Our stability result applies equally well to arbitrary type of (long-range)
spin interactions. For example, concerning the XYZ spin interactions with
$\Delta\neq 0$ in Eq. (1), our derivation can be straightforwardly
generalized. The only addition we need is for terms arising from
$\tilde{S}^{z}_{k}\tilde{S}^{z}_{-k}$ in $\hat{V}_{k}$: The additional rotated
operators read

$\tilde{S}^{z}_{\pm k}(\bm{\Omega})=R_{z,+}(\bm{\Omega})\tilde{S}^{+}_{\pm
k}+R_{z,-}(\bm{\Omega})\tilde{S}^{-}_{\pm
k}+R_{z,z}(\bm{\Omega})\tilde{S}_{\pm k}^{z}$ (76)

with coefficients computed analogously:

$\begin{split}R_{z,+}(\bm{\Omega})=&R_{z,-}^{*}(\bm{\Omega})=-\sin(\theta)e^{i\phi},\\\
R_{z,z}(\bm{\Omega})=&+\cos(\theta).\\\ \end{split}$ (77)

These extra terms just result in the addition $+2\Delta R_{z,\mu}R_{z,\nu}$
inside the square brackets in Eq. (52). The rest of the discussion is
unchanged.

####  V.8.3 Higher spins

Our stability result becomes stronger when individual spins are larger. In
fact, an equivalent analysis can be performed when spins-$1/2$ are replaced by
spins-$s$ with arbitrary $s$. This replacement gives more room to the accuracy
of the bosonic description of spin states, which requires $\delta_{S}\ll Ls$.

The only important caveat here is the absence of self-interactions
$J_{\mathbf{r},\mathbf{r}}=0$ (which model e.g. single-ion anisotropy in
magnetic materials with higher local spins). We will comment on their effect
in Sec. VI below.

To summarize, our analysis straightforwardly carries over to arbitrary
lattices, multi-spin interactions, and spin sizes, resulting in the same
qualitative conclusion in all cases.

##  VI Predictions: eigenstate delocalization by semiclassical chaos

The theory of Sec. V predicts the stability of QMBS with i) sufficiently long
range of interactions and ii) absence of spectral resonances for $\alpha=0$.
This theory rationalizes the numerical findings of Sec. IV for the long-range
quantum Ising chain with $\alpha<1$ and $h\geq J_{0}$, which generically
satisfies both conditions.

What if conditions i) or ii) fail? In Sec. IV we already presented evidence
that violation of the condition i) (i.e. $\alpha>1$) results in asymptotic
hybridization of the seemingly anomalous finite-size eigenstates. Concerning
condition ii), the representation of the Hamiltonian in Sec. V.5 in terms of a
collective spin coupled to magnon modes gives us a clear and predictive
understanding of QMBS stability. In this Section we explore situations where,
according to our theory, stability is challenged by resonances in the mean-
field energy spectrum.

We begin by discussing proximity to rare accidental resonances in Sec.VI.1,
which generically exist for large system size. Afterwards, we demonstrate how
exponentially diverging classical trajectories of the mean-field Hamiltonian
result in instability of QMBS. This happens in correspondence of isolated
phase-space separatrices (often associated with ESQPT) discussed in Sec. VI.2,
as well as of fully developed classical chaos, e.g. generated by strong
periodic drives, which we analyze in Sec. VI.3. In all these cases numerical
data indicate instability of QMBS in agreement with our theory predictions.

Figure 12: QMBS hybridization due to an accidental spectral resonance. Data
are shown for the quantum Ising chain with $L=14$, $\alpha=0.7$, $J_{0}=1$ and
$h=1.01$, $1.02$ and $1.03$. Scatter plot of the observables
$\langle\delta_{S}\rangle_{E_{n}}$ (row a) and entanglement entropy
$S_{L/2}(E_{n})$ (row b) of each eigenstate $|E_{n}\rangle$ versus energy
density $E_{n}/L$.

###  VI.1 Accidental resonances

The derivation of a self-consistent exact solution for QMBS in Sec. V.6
required the absence of resonances $\bar{\omega}\neq r\omega$ between the two
fundamental rotor and magnon frequencies $\omega$ and $\bar{\omega}$. This is
a generic feature of energy shells of Hamiltonians with collective spin-spin
interactions. However, accidental resonances may occur for specific values of
model parameters and energy density. In this case the QMBS are expected to
hybridize with the rest of the spectrum. Observing such rare (zero-measure)
instances in finite systems requires a careful search in parameter space.

An example of a fine-tuned resonance is shown in the middle panel of Fig. 12,
where the spin depletion and the entanglement entropy of each eigenstate are
plotted for the quantum Ising chain Hamiltonian with $L=14$ and $\alpha=0.7$.
The accidental resonance takes place at $h=1.02J_{0}$ (central panel), causing
significant hybridization of one of the QMBS (large orange dots). The spikes
in the spin depletion $\braket{\delta_{S}}\simeq 1.2$ and entanglement entropy
$S_{L/2}$ for the central scar eigenstate are well visible. This effect
disappears upon an slight change in any parameter. This is illustrated in the
two lateral panels of the same Figure, where we report analogous plots with
$\approx 1\%$ variations of the transverse field, i.e., $h/J_{0}=1.01$ and
$1.03$.

###  VI.2 Excited-state quantum phase transitions

In this Section, we discuss another effect leading to instability of QMBS:
mean-field criticality at finite energy density. This happens whenever the
classical mean-field Hamiltonian $\mathcal{H}$ has saddle points, which are
accompanied by isolated unstable trajectories terminating on them (phase-space
separatrices) with diverging period. These systems are thus characterized by
singularities of the density of states at some finite energy density
$E_{\text{cr}}/L$. Such mean-field criticality can (but does not necessarily)
result from the spontaneous breaking of a discrete symmetry. In this case the
critical energy separates ordered eigenstates (below) from disordered ones
(above).

Such excited-state quantum phase transitions (ESQPTs) Cejnar  _et al._ (2021)
have been discussed at mean-field level, see also Refs. Stránskỳ  _et al._
(2015); Santos  _et al._ (2016b); Kloc  _et al._ (2017). The accumulation of
energy eigenvalues at a given energy density is necessarily associated with
resonances – in the language of Sec. V.5 one has $\omega(n)=0$. This leads to
the failure of eigenstate localization theory. We are thus led to conjecture
that quantum many-body chaos develops for $\alpha>0$ near the location of a
mean-field ESQPT.

We numerically assess the impact of finite-range interactions on an ESQPT. In
Fig. 13 we consider the quantum Ising chain Hamiltonian in Eq. (4) with
$\gamma=1$ and in the ordered phase $h<J_{0}$. In this range of parameters the
level spacing distribution agrees with Wigner-Dyson statistics for $\alpha>0$,
in agreement with the findings of Refs. Fratus and Srednicki (2016);
Russomanno  _et al._ (2021). In Fig. 13 we take $h=0.4J_{0}$ and $\alpha=0.7$
and check the stability of QMBS for $L=14$, $16$, $18$. In the energy region
corresponding to the transition, QMBS are unstable and hybridize with the rest
of the spectrum, while states well above or below the transition remain
stable. We find indications of scars hybridization even for small $\alpha>0$.

Figure 13: Fate of the ESQPT for finite-range interactions $\alpha>0$. Scatter
plot of the observables $\langle\delta_{S}\rangle_{E_{n}}$ (first row) and
entanglement entropy $S_{L/2}(E_{n})$ (second row) of each eigenstate
$|E_{n}\rangle$ versus energy density $E_{n}/L$, for $\alpha=0.7$, $h=0.4$,
with $L=14$, $16$, $18$ from left to right. The dashed vertical line indicates
the location $E_{\text{cr}}/L$ of the ESQPT in the mean-field limit
$\alpha=0$.

Interestingly, we observe that the energy window of instability of QMBS shows
a marked asymmetry, lying entirely _above_ the mean-field energy density value
$E_{\text{cr}}/L$ associated with the ESQPT (dashed vertical lines in Fig.
13). This global displacement towards higher energies can be attributed to
finite-size effects, as the window slowly drifts downwards in energy (i.e.
leftwards in the Figure) upon increasing $L$. However, in the weak coupling
regime (small $\alpha$), one can give an intuitive physical argument for the
persistence of the asymmetry for large $L$ based on energy conservation as
follows. As the collective spin is weakly coupled to empty bosonic modes, the
latter _subtract_ energy from the former, pushing it downward in the mean-
field energy landscape. Thus, for energy slightly below $E_{\text{cr}}/L$,
interactions push the collective spin _away_ from the resonance, leading to
stability; instead, for energy above $E_{\text{cr}}/L$, they push it _towards_
the resonance, leading to instability Lerose  _et al._ (2018, 2019a).

Figure 14: Fate of QMBS in the kicked variable-range quantum Ising model. Here
$L=18$. Scatter plot of the observables $\langle\delta_{S}\rangle_{E_{n}}$
(row a), entanglement entropy $S_{L/2}(E_{n})$ (row b), and overlap with the
polarized state $|\langle\uparrow|E_{n}\rangle|^{2}$ (row c), for each
eigenstate $|E_{n}\rangle$, versus energy density $E_{n}/L$. Column 1 has
$J_{0}=0.5$, corresponding to the weak-integrability breaking (KAM) regime,
and $\alpha=0.8$. Columns 2-5 have $J_{0}=8$, corresponding to the fully
chaotic regime of the mean-field dynamics, and increasing $\alpha=0$, $0.2$,
$0.5$, $0.8$, respectively.

###  VI.3 Semiclassical chaos: Periodic kicking

A cornerstone of our stability theory in Sec. V is classical integrability of
the mean-field Hamiltonian, which implies generic absence of resonances in the
unperturbed spectrum. We here discuss the effects on QMBS for $\alpha>0$
caused by the breaking of classical integrability at mean-field level
$\alpha=0$. The tight relationship between classical trajectories and
semiclassical eigenstates is well established Percival (1973). When the mean-
field integrability-breaking perturbation is small, Kolmogorov-Arnold-Moser
(KAM) theory Giorgilli (2022) guarantees that the most classical trajectories
are deformed tori. In this case the semiclassical spectrum is regular. Our
theory thus predicts generic stability of QMBS for finite $0<\alpha<d$.
Conversely, for stronger integrability breaking the classical phase space
gradually crosses over to fully developed chaos. In this case the
semiclassical spectrum for $\alpha=0$ is well-known to display statistical
properties compatible with random-matrix universality. As energy eigenvalues
are randomly arranged according to Wigner-Dyson statistics and uncorrelated
between adjacent towers in $\delta_{S}$, low-order resonances become generic,
and we expect full eigenstate delocalization for small $\alpha>0$.

We numerically verify the above prediction by studying the kicked version of
the quantum Ising model, with a time-periodic Hamiltonian defined by the
following two-step protocol,

$\hat{H}_{\alpha}(t)=\begin{dcases}-\frac{J_{0}}{\mathcal{N}_{\alpha,L}}\sum_{j=1}^{L}\sum_{r=1}^{L/2}\frac{\hat{\sigma}^{x}_{j}\hat{\sigma}^{x}_{j+r}}{r^{\alpha}}\quad
t\in\left[-\frac{T}{4},\frac{T}{4}\right)\mod T\\\
-h\sum_{j=1}^{L}\hat{\sigma}^{z}_{j}\qquad\qquad
t\in\left[\frac{T}{4},\frac{3}{4}T\right)\mod T\ ,\end{dcases}$ (78)

where $T$ is the period of the drive. The evolution over one period reads

$\displaystyle\begin{split}\hat{U}&=T_{t}\exp\left(-i\int_{0}^{T}H(t)dt\right)=e^{i\hat{W}_{\alpha}\frac{T}{4}}e^{ih\hat{S}^{z}\frac{T}{2}}e^{i\hat{W}_{\alpha}\frac{T}{4}}\
,\end{split}$ (79)

where
$\hat{W}_{\alpha}=\frac{J_{0}}{\mathcal{N}_{\alpha,L}}\sum_{j=1}^{L}\sum_{r=1}^{L/2}\frac{\hat{\sigma}^{x}_{j}\hat{\sigma}^{x}_{j+r}}{r^{\alpha}}$.
In our numerical simulations, we fix $T=1$ and $h=\pi/2$. For $\alpha=0$ this
is a textbook model of quantum chaos known as _quantum kicked top_. Depending
on the value of the interaction strength $J_{0}$, the model exhibits a
transition between a regular phase-space – described by KAM theory – and a
chaotic one Haake  _et al._ (1987); Haake (1991). Its dynamical properties and
their relations to quantum information dynamics have been intensively
investigated Zarum and Sarkar (1998); Miller and Sarkar (1999); Chaudhury  _et
al._ (2009); Piga  _et al._ (2019); Wang  _et al._ (2004); Trail  _et al._
(2008); Ghose  _et al._ (2008a); Lombardi and Matzkin (2011); Stamatiou and
Ghikas (2007); Ghose  _et al._ (2008b); Madhok  _et al._ (2014); Pappalardi
_et al._ (2018); Sieberer  _et al._ (2019); Pilatowsky-Cameo  _et al._ (2020);
Pappalardi and Kurchan (2023).

We perform exact diagonalization of the unitary operator in Eq. (79) and study
the Floquet spectrum

$\hat{U}\ket{\phi_{n}}=e^{i\phi_{n}}\ket{\phi_{n}}\ ,$ (80)

with $\phi_{n}\in(-\pi,\pi]$. We examine the behaviour of Floquet eigenstates
as a function of the chaoticity parameter $J_{0}$ and of $\alpha$. This is
illustrated in Fig. 14, where – paralleling the analysis of the time-
independent Hamiltonian version in Fig. 5 – we report a scatter plot the
collective spin depletion, the entanglement entropy, and the overlap with a
fully polarized state, for each Floquet eigenstate, versus the Floquet quasi-
energy $\phi_{n}/\pi$, with fixed $L=18$.

Results are remarkably sharp even for the small system size we can simulate.

For small $J_{0}=0.5$ the mean-field classical phase-space is regular.
Correspondingly, the spectral properties in this regime (column 1 of Fig. 14)
show clear signatures of stability of QMBS for all $0<\alpha\lesssim 1$. The
conditions identified by our theory in Sec. V thus appear necessary and
sufficient for robust QMBS to exist. By contrast, for $J_{0}=8$ and
$\alpha=0$, the classical phase-space is fully chaotic. The quantized single-
spin spectra within sectors labelled by $\delta_{S}$ thus exhibit random-
matrix features (while retaining exact degeneracy among equal-weight $SU(2)$
irreducible representation spaces). This is shown in column 2 of Fig. 14. As
soon as permutational symmetry is broken by finite $\alpha>0$ (columns 3-5 of
Fig. 14) we observe that the eigenstates quickly hybridize, in a way
compatible with the onset of strong ETH. In particular, _no signatures of QMBS
phenomena appear_. We conclude that classical chaos at the mean-field level
completely destroys the mechanism of eigenstate localization and suppresses
all signatures of QMBS, for arbitrarily slow decay of the interactions.

We further analyze the behavior upon increasing the system size $L$. In Fig.
15 we plot the collective spin depletion
$\braket{\hat{\delta}_{S}}=L/2-\langle|\hat{\mathbb{S}}|\rangle$ averaged over
the $L/4+1$ eigenstates with minimum $\braket{\hat{\delta}_{S}}$ (i.e., the
candidate QMBS). As shown in the left panel, in the strongly chaotic regime
$\braket{\hat{\delta}_{S}}$ increases rapidly with $L$, while the growth is
strongly suppressed in the regular regime. In this case the scaling of
$\braket{\hat{\delta}_{S}}$ with $L$ looks roughly compatible with the
behavior predicted by the theory of eigenstate localization, Eq. (73),
highlighted by the dashed lines in the right panel.

As an important byproduct of our analysis, we have demonstrated that _in
periodically driven quantum many-body systems with slowly decaying
interactions and non-chaotic mean-field dynamics, energy absorption from the
drive and entropy growth (“heating”) is suppressed for arbitrarily long
times_. Conversely, in presence of chaotic mean-field dynamics, heating is
fast even for small $\alpha>0$. This demonstration underlies the findings of
Ref. Lerose  _et al._ (2019c). Suppression of heating opens the door to
hosting long-range ordered non-equilibrium phases of matter such as discrete
time crystals Pizzi  _et al._ (2021); Giachetti  _et al._ (2022) protected by
QMBS stability, similarly to Refs. Maskara  _et al._ (2021); Bluvstein  _et
al._ (2021). See also Ref. Daniel  _et al._ (2023).

Figure 15: Collective spin depletion averaged over the $L/4+1$ most polarized
Floquet eigenstates of the kicked quantum Ising chain vs system size. Left:
Comparison between the behavior in the KAM ($J=0.5$, dashed) and chaotic
($J=8$, solid) regimes, for three values of $\alpha$ (see legend). In the
right panel we zoom on the former, and compare them with the theory prediction
in Eq. (73) (dashed lines).

###  VI.4 Higher spin models

Higher spin Hamiltonians with $s>1/2$ can exhibit chaotic mean-field dynamics
even without an external drive. This happens when the Hamiltonian features
additional “self-interaction” terms with $\mathbf{i}=\mathbf{j}$ in Eq. (1),
for example $+D\sum_{{\mathbf{j}}}(\hat{\sigma}_{\mathbf{j}}^{z})^{2}$. Such
terms break the conservation of the collective spin magnitude
$\hat{\mathbb{S}}^{2}$. In this case the $\alpha=0$ Hamiltonian involves
$n=2s>1$ collective degrees of freedom Sciolla and Biroli (2011), and its
classical limit may display chaotic behavior depending on the Hamiltonian
parameters.

Our theory thus suggests that spin-$s$ Hamiltonians with slowly decaying
interactions $0<\alpha<d$ _and_ self-interactions exhibit a similar
phenomenology as the periodically kicked model discussed in Sec. VI.3: QMBS
arising from $\alpha=0$ KAM regimes and strong ETH arising from $\alpha=0$
fully chaotic regimes. Note that mean-field integrability in this enlarged
setup would lead to a generalization of our effective Hamiltonian in Sec. V.5,
consisting of multiple rotors coupled to magnon modes. We envisage that this
generalized effective Hamiltonian can be approached with the same technique
and leads to self-consistent eigenstate localization under the assumption of
strong incommensurability of all the unperturbed fundamental frequencies.

##  VII Conclusions and perspectives

In this paper we investigated the existence of _robust_ quantum many-body
scars (QMBS) in quantum spin lattices. Our theory identifies a stability
mechanism of QMBS based on two conditions: sufficiently slow decay of
isotropic interactions $0<\alpha<d$ and classical integrability of the
underlying few-body mean-field Hamiltonian. We provided extensive numerical
evidence supporting our analytical predictions.

The results of this work have broader consequences:

  * •

Dynamics in mean-field models ($\alpha=0$) are well understood in terms of the
classical limit of few collective degrees of freedom Sciolla and Biroli
(2011). Mean-field integrability-breaking gives rise to the established
Kolmogorov-Arnold-Moser (KAM) scenario Giorgilli (2022), whereby non-ergodic
motion is largely stable for small perturbations, and chaos gradually invades
phase space as the perturbation gets stronger. However, no general KAM-like
scenario has been established for systems with many interacting degrees of
freedom. On the contrary, it is largely believed that generic infinitesimal
perturbations completely destroy integrability. A crucial question remains
open on whether specific classes of interacting (classical or quantum) many-
body systems allow for a KAM-like scenario. Long-range interacting systems are
natural candidates, as the parameter $\alpha$ allows to interpolate between
few-body classical dynamics ($\alpha=0$) and conventional many-body quantum
dynamics ($\alpha=\infty$). The findings of this paper imply _a sense of
structural stability of non-ergodic dynamics arising from mean-field
integrability, reminiscent of the KAM scenario_. To the best of our knowledge
this might be the first controlled extension of KAM-like phenomena to many-
body physics.

  * •

As a byproduct, we have shown when and how _heating can be suppressed in
periodically driven long-range interacting spin systems_. We found that
regularity of the mean-field phase space for $\alpha=0$ determines stability
of large-spin eigenstates to finite-range interactions $0<\alpha<d$. Existence
of such robust QMBS implies bounded energy absorption from the drive and
bounded entropy growth.

  * •

Excited state quantum phase transitions (ESQPTs) are usually discussed only in
the context of mean-field models ($\alpha=0$). Our results elucidate the _fate
of ESQPTs upon decreasing the interaction range_ : For finite $0<\alpha<d$ the
unavoidable spectral resonances at the ESQPT cause hybridization of the large-
spin eigenstates in a finite spectral window right above the critical energy.
This work thus provides a step towards a possible many-body theory of ESQPTs.

  * •

The class of QMBS wavefunctions discovered here extends the scope of
collectively entangled states arising from all-to-all interactions, such as
Dicke states (which are actually akin to the known exact QMBS wavefunctions
Moudgalya  _et al._ (2022)). For instance, theoretical bounds on multipartite
entanglement of these states will inform extensions of the quantum metrology
and sensing toolbox Braunstein and Caves (1994); Giovannetti  _et al._ (2011);
Tóth (2012); Hyllus  _et al._ (2012); Pezzè  _et al._ (2018) as already
pointed out for other instances of QMBS Dooley (2021); Desaules  _et al._
(2022b); Dooley  _et al._ (2023). Experimental applications can be envisaged
in the atomic-molecular-optical platforms featuring slowly decaying
interactions Blatt and Roos (2012); Britton  _et al._ (2012); Richerme  _et
al._ (2014). These include dipolar lattices in two but not in three
dimensions, as the strong anisotropy of interactions has been theoretically
and experimentally shown to cause rapid thermalization Lepoutre  _et al._
(2019); Alaoui  _et al._ (2022).

This work furthermore opens future avenues of investigation:

  * -

Our approach — more specifically Eqs. (53) and (60) — may form the basis of _a
new class of computational methods_ to study non-equilibrium quantum spin