# Out-of-equilibrium dynamics of quantum many-body systems with long-range
interactions

Nicolò Defenu  Institute for Theoretical Physics, ETH Zürich, Wolfgang-Pauli-
Str. 27, 8093 Zürich, Switzerland<EMAIL_ADDRESS>Alessio Lerose
Department of Theoretical Physics, University of Geneva, Quai Ernest-Ansermet
30, 1205 Geneva, Switzerland<EMAIL_ADDRESS>Silvia Pappalardi
Institut für Theoretische Physik, Universität zu Köln, Zülpicher Straße 77,
50937 Köln, Germany<EMAIL_ADDRESS>
###### Abstract

Experimental progress in atomic, molecular, and optical platforms in the last
decade has stimulated strong and broad interest in the quantum coherent
dynamics of many _long-range interacting_ particles. The prominent collective
character of these systems enables novel non-equilibrium phenomena with no
counterpart in conventional quantum systems with local interactions. Much of
the theory work in this area either focussed on the impact of variable-range
interaction tails on the physics of local interactions or relied on mean-
field-like descriptions based on the opposite limit of all-to-all infinite-
range interactions. In this Report, we present a systematic and organic review
of recent advances in the field. Working with prototypical interacting quantum
spin lattices without disorder, our presentation hinges upon a versatile
theoretical formalism that interpolates between the few-body mean-field
physics and the many-body physics of quasi-local interactions. Such a
formalism allows us to connect these two regimes, providing both a formal
quantitative tool and basic physical intuition. We leverage this unifying
framework to review several findings of the last decade, including the
peculiar non-ballistic spreading of quantum correlations, counter-intuitive
slowdown of entanglement dynamics, suppression of thermalization and
equilibration, anomalous scaling of defects upon traversing criticality,
dynamical phase transitions, and genuinely non-equilibrium phases stabilized
by periodic driving. The style of this Report is on the pedagogical side,
which makes it accessible to readers without previous experience in the
subject matter.

###### Contents

  1. 1 Introduction
  2. 2 Equilibrium properties of long-range interacting quantum spin systems
    1. 2.1 Variable-range quantum XY model
    2. 2.2 Equilibrium phase diagram in a nutshell
    3. 2.3 Low-energy spectrum with infinite-range interactions ($\alpha=0$)
      1. 2.3.1 Mean-field theory as an exact classical limit
      2. 2.3.2 Collective quantum fluctuations and excitations
      3. 2.3.3 “Spin-wave” excitations
    4. 2.4 Finite-range interactions ($\alpha>0$)
      1. 2.4.1 Perturbation to mean-field
      2. 2.4.2 Quantum paramagnetic phase
      3. 2.4.3 Quantum ferromagnetic phase
    5. 2.5 Structure of the spectrum beyond linear spin-wave theory
  3. 3 Low-energy dynamics
    1. 3.1 Spreading of correlations
      1. 3.1.1 Weak long-range regime ($\alpha>d$)
      2. 3.1.2 Strong long-range regime ($0<\alpha<d$)
      3. 3.1.3 Other directions
    2. 3.2 Metastability
      1. 3.2.1 State of the art
      2. 3.2.2 Quasi-stationary states and spectral properties
      3. 3.2.3 Equilibration in presence of disorder
    3. 3.3 Kibble-Zurek mechanism
      1. 3.3.1 State of the art
      2. 3.3.2 Quasi-static dynamics for $\alpha=0$
      3. 3.3.3 Adiabaticity breaking
      4. 3.3.4 Full counting statistics of defects
    4. 3.4 Dynamical quantum phase transitions - Loschmidt echo
      1. 3.4.1 Spherical spin-wave theory
      2. 3.4.2 Step approximation
      3. 3.4.3 Strong long-range regime
  4. 4 Dynamics in highly excited states
    1. 4.1 Quench dynamics of fully-connected spin systems ($\alpha=0$)
      1. 4.1.1 Mean-field classical limit
      2. 4.1.2 Dynamical phase transitions - Order parameter
      3. 4.1.3 Semiclassical dynamics of quantum fluctuations
      4. 4.1.4 Scrambling dynamics
      5. 4.1.5 Entanglement dynamics
    2. 4.2 Quench dynamics of long-range interacting spin systems ($\alpha>0$)
      1. 4.2.1 Dynamics of quantum fluctuations with finite interaction range
      2. 4.2.2 Prethermal freezing of spin-wave excitations
      3. 4.2.3 Impact of finite-range interactions on dynamical phase transitions
      4. 4.2.4 Scrambling dynamics with variable-range interactions
      5. 4.2.5 Entanglement entropy dynamics: Spin-squeezing vs Quasiparticle picture
  5. 5 Dynamical phases induced by periodic driving
    1. 5.1 Kapitza phases
      1. 5.1.1 Fully-connected limit $\alpha=0$: Non-equilibrium phases by driving
      2. 5.1.2 Quantum many-body Kapitza phases for $\alpha>0$
      3. 5.1.3 Prethermalization and heating
    2. 5.2 Discrete time crystals
      1. 5.2.1 Mean-field DTC
      2. 5.2.2 Finite-size and finite-range effects
      3. 5.2.3 Order parameter
  6. 6 Conclusions and perspectives
  7. A Semiclassical spectrum
  8. B Semiclassical description of fully-connected systems
    1. B.1 Semiclassical approach
    2. B.2 Classical limit of permutationally invariant systems
    3. B.3 Beyond global permutational symmetry
  9. C Asymptotic estimates for ${f}_{\mathbf{k}}(\alpha)$
    1. C.1 Strong long-range regime ($0<\alpha<d$)
    2. C.2 Weak long-range regime ($\alpha>d$)
  10. D Exact solution of quasi-static drive for a single mode
    1. D.1 Fidelity and defect density
    2. D.2 Slow quench to the critical point
    3. D.3 Full ramp
  11. E Entanglement dynamics in other long-range interacting models
    1. E.1 Quantum Ising chain in a tilted field
    2. E.2 Short-range perturbations to collective spin models
  12. F Floquet Hamiltonian and high-frequency expansion

##  1 Introduction

An increasing interest in quantum many-body physics with long-range
interactions is being driven by growing experimental capabilities in
controlling and manipulating atomic, molecular, and optical systems (AMO).
Currently, various platforms such as Rydberg atoms, dipolar quantum gases,
polar molecules, quantum gases in optical cavities, and trapped ions, have
native two-body long-range interactions which can be modeled as algebraically
decaying $J/(\Delta r)^{\alpha}$ with the distance $\Delta r$ [1, 2, 3, 4, 5,
6]. The exponent $\alpha$ can in some cases be experimentally tuned — e.g.
through off-resonant coupling of internal levels of trapped ions to motional
degrees of freedom, or trapping neutral atoms in photonic modes of a cavity.
Additionally, the effective interaction range can be efficiently tuned in
systems of Rydberg atoms in one- and two-dimensional arrays [7, 8] or by
Rydberg dressing [9, 10].

The versatility of the aforementioned AMO platforms spurred intense
theoretical and experimental explorations. These studies established that
long-range interactions provide clear routes to circumventing the constraints
imposed by either conventional thermalization [11] or conventional bounds on
information spreading [12]. Accordingly, the prominent collective character of
systems with long-range interactions can lead to a kaleidoscope of novel
phenomena which cannot be observed in systems with local interactions. Major
examples include: the observation of “super-luminal” correlation and
entanglement spreading [13, 14, 3] (to be contrasted with the conventional
light-cone behavior in presence of local interactions [15]); dynamical phase
transitions [16, 17, 18, 19, 20, 21, 22]; exotic defect scaling [23, 24];
self-organized criticality [25]; time-translation symmetry breaking [26, 27,
28]; quantum many-body chaos [29, 30, 31]. As such, control of long-range
interacting assemblies stands out as a promising ingredient for future
quantum-technological applications, including quantum metrology and quantum
computation.

While this great diversity of platforms and research directions largely
contributes to generate widespread excitement about long-range interactions,
it has at the same time certain drawbacks. The backgrounds and interests of
the numerous research groups active in this area span a very wide spectrum. On
one hand, experimental interpretations are often based on a few-body, mean-
field-like way of thinking [32, 33, 34, 35]. Albeit remarkably simple and
powerful, this perspective may fail to fully capture the complexity of non-
equilibrium phenomena with long-range interactions. On the other hand,
theoretical investigations have often prioritized mathematically rigorous
efforts aimed at characterizing the departure from known properties of
locally-interacting systems [36, 37, 38, 39, 40, 41, 42, 43, 44]. Albeit
sometimes in synergy with experiments [45], this perspective may obscure the
construction of an intuitive physical picture applicable to the broad range of
out-of-equilibrium phenomena mentioned above. Despite recent attempts to
recompose the corresponding mosaic in equilibrium [46], this complementarity
of perspectives on similar phenomena still struggle to come together and
cement a unified research field and community. As a consequence, the current
understanding of the out-of-equilibrium dynamics of long-range interacting
quantum many-body systems still seems to lack a systematic organization
comparable to that of quantum locally interacting [11] or classical long-range
interacting [47] systems.

The purpose of this Report is to provide a systematic and intuitive
theoretical approach to non-equilibrium phenomena arising from non-random
long-range interactions in quantum many-body systems. Our effort aims at
bridging the various complementary views in this wide research area and
creating a unifying framework. We will review a selection of significant
findings in the field, emphasizing how they can be encompassed within a common
basic theoretical language and formalism. The approach reviewed in this Report
is suited to bridge the simple mean-field description — which applies to
infinite-range interactions, i.e. $\alpha=0$ — to the description of systems
with quasi-local interactions, i.e. $\alpha\gg d$, which allow a well-defined
notion of locality in spite of non-local interaction tails. The _strong long-
range regime_ in between will be the focus of this Report; we will frequently
emphasize the leitmotiv that _the physics in this regime interpolates between
conventional few-body and conventional many-body physics_. The reach of this
unifying framework will be illustrated using prototypical models of
interacting quantum spin lattices. This choice does not only serve the purpose
of directly relating our results to paradigmatic locally-interacting systems
[48, 49], but it is also allows us to make direct connections with the major
AMO experimental platforms recalled above.

This Report is organized as follows:

  * •

Our journey will start in Sec. 2 with a review of _equilibrium properties_ of
ferromagnetic quantum spin systems exemplified by a variable-range quantum XY
model 2.1, including a discussion of the equilibrium phase diagram upon
varying parameters and interaction range (via $\alpha$) 2.2 and a critical
examination of the mean-field limit 2.3. Hence, in Sec. 2.4 we will review the
low-energy description in terms of bosonic excitations (spin waves) across the
phase diagram, with emphasis on spectral properties arising from a long
interaction range. Finally in Sec. 2.5 we will discuss spectral properties
beyond linear spin-wave theory.

  * •

The low-energy description reviewed in Sec. 2 can be used to investigate near-
equilibrium dynamics. This setup allows to study the peculiar properties of
spatial propagation of quantum correlation in presence of long-range
interactions [50, 51] as well as their unusual equilibration dynamics [52],
reviewed in Sec. 3.1 and 3.2 respectively. Both these phenomena can be studied
in quantum quenches lying within the supercritical phase and, therefore, only
represent a small departure from equilibrium. More surprisingly, we are going
to show that the low-energy description is also capable of addressing
dynamical scaling phenomena arising after quenches across the critical point.
This is the case of the universal defect formation following a quasi-static
sweep across the quantum critical point [53, 54], see Sec. 3.3 and the rise of
dynamical quantum phase transitions [55], which we will treat in Sec. 3.4 (the
latter, however, will require to modify the simple low-energy description
employed before).

  * •

Section 4 is devoted to the study of _dynamics far away from equilibrium_ ,
induced by _quantum quenches_. We will first consider in Sec. 4.1 the fully-
connected model with all-to-all uniform interactions, and examine the simplest
instances of dynamical phenomena in this limit, which can be understood in
terms of few-body semiclassical dynamics. For finite-range interactions,
however, the motion of semiclassical collective degrees of freedom is coupled
to many quantum-fluctuation modes with various wavelengths. In Sec. 4.2 we
review a systematic approach to the resulting complex many-body problem,
originally developed in Refs. [56, 57]. As a first implication stemming from
this approach we will review lower bounds on thermalization time scales
associated with long-range interactions, establishing the genuinely non-
equilibrium nature of dynamical phenomena in these systems. Hence, we will
examine the impact of many-body quantum fluctuations on dynamical criticality
and quantum information spreading far away from equilibrium upon tuning the
interaction range. Throughout, we will highlight the role of long-range
interactions in generating novel phenomena.

  * •

In Sec. 5 we will employ the methodology of Sec. 4 to describe coherent
dynamics subject to _periodic driving_. Here we will review how long-range
interactions allow to stabilize genuinely non-equilibrium phases, without an
equilibrium counterpart, in low-dimensional quantum systems of the kind
routinely realized in AMO experiments. This will include phases that may be
viewed as quantum many-body realizations of the celebrated Kapitza pendulum
(Sec. 5.1) as well as discrete-time crystals, which spontaneously break time-
translation symmetry (Sec. 5.2).

  * •

Finally, in the conclusive Section, we spell out the topics which are _not_
covered in this Report: from effects of inhomogeneities, to frustrated,
random, or noisy interactions, to dissipative and monitored dynamics.

Throughout the presentation, our goal is to provide both physical intuition
and systematic theoretical understanding of experimentally relevant phenomena.
We kept the style of the Report on the pedagogical side, as we hope this work
will also be useful to inexperienced readers who only possess a ground
knowledge of quantum many-body theory and are interested in taking their first
dive into the realm of quantum dynamics in presence of long-range
interactions.

##  2 Equilibrium properties of long-range interacting quantum spin systems

In this Section we summarize and discuss basic _equilibrium_ properties of
quantum spin lattices with variable interaction range, which will come in
useful in the rest of this Report. For definiteness we will focus on a class
of ferromagnetic XY quantum spin models, introduced in Sec. 2.1, and review
its equilibrium phase diagram in Sec. 2.2. We will then work out its low-
energy description in terms of bosonic excitations (“spin waves”) both in the
fully-connected limit (Sec. 2.3) and with finite-range interactions (Sec.
2.4), with emphasis on the peculiar features of the quasiparticle spectrum
such as discreteness [58, 52], divergent group velocity [59, 60, 50, 61], and
dressing effects. Finally, in Sec. 2.5 we will discuss finer low-energy
properties beyond spin-wave theory, including domain-wall (de)confinement [62,
63].

In this Section we will keep the model parameters fully general. In the rest
of the Report we will frequently restrict the model for simplicity, but all
the results can always be straightforwardly extended, drawing on the general
setup introduced here and in Sec. 4.2.1 below.

###  2.1 Variable-range quantum XY model

Throughout this work we will consider a prototypical model implemented in AMO
platforms, a quantum XY spin model with tunable interaction range. We take a
$d$-dimensional square lattice of $N=L^{d}$ quantum spins-$s$ described by a
Hamiltonian of the form

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

In this equation
$\hat{\sigma}^{\mu}_{\mathbf{r}}=\hat{s}^{\mu}_{\mathbf{r}}/s$ are operators
corresponding to the normalized spin components in the $\mu=x,y,z$ direction,
acting on site $\mathbf{r}=(r_{1},\dots,r_{d})$ of the lattice,
$r_{1},\dots,r_{d}=1,\dots,L$. This represents a generalization of the
standard spin-$1/2$ case, where $\hat{\sigma}^{\alpha}_{\mathbf{r}}$ reduce to
the standard Pauli matrices. Such a normalization allows us to keep track of
the role of quantum fluctuations, which are suppressed in the classical limit
$s\to\infty$. The quantity $\gamma$ parametrizes the XY anisotropy. In this
Report we will consider anisotropic spin systems, i.e. $\gamma\neq 0$; for
definiteness we assume $\gamma>0$ (negative values equivalent upon rotating
the spins around the $z$-axis by $\pi/2$), and we will frequently set
$\gamma=1$ (quantum Ising model). We will occasionally comment on the
isotropic limit $\gamma\to 0$ when relevant. The quantity $h$ represents the
transverse magnetic field strength, which we assume $h\geq 0$ (negative values
are equivalent upon rotating the spins around the $x$-axis by $\pi$).
Throughout this report we will always use units such that Planck’s constant is
$\hbar\equiv 1$.

The ferromagnetic couplings $J_{\mathbf{r},\mathbf{r^{\prime}}}(\alpha)\equiv
J_{\Delta r}(\alpha)$ depend on the distance $\Delta
r=\lvert\lvert\mathbf{r}-\mathbf{r^{\prime}}\rvert\rvert$ between the two
involved sites, and we will be interested in tuning their spatial range
through the parameter $\alpha$. Specifically, we consider interactions
algebraically decaying with the distance,111Note that for spins-$1/2$ the
terms $\mathbf{r}=\mathbf{r^{\prime}}$ produce an inconsequential additive
constant $E=\sum_{\mathbf{r}}J_{\mathbf{r},\mathbf{r}}/2$, as Pauli matrices
square to $1$. Diagonal terms may be important for higher-spin Hamiltonians.
As a rule we will set $J_{\mathbf{r},\mathbf{r}}=0$. We will occasionally
comment on interesting phenomena associated with spin self-interactions
further below.

$J_{\mathbf{r},\mathbf{r^{\prime}}}(\alpha)=\frac{J}{\lvert\lvert\mathbf{r}-\mathbf{r^{\prime}}\rvert\rvert^{\alpha}}.$
(2)

To impose periodic boundary conditions, various equivalent choices of distance
function are possible; we take
$\lvert\lvert\mathbf{r}-\mathbf{r^{\prime}}\rvert\rvert\equiv\sqrt{\sum_{\mu=1}^{d}[\min(|r_{\mu}-r^{\prime}_{\mu}|,L-|r_{\mu}-r^{\prime}_{\mu}|)]^{2}}$.
The overall constant $J>0$ is chosen in such a way to fairly compare the
models with different $\alpha$, i.e. to make the mean-field interaction
strength $J_{0}$ independent of $\alpha$:

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

This prescription — known as Kac normalization [64] — is necessary to make the
thermodynamic limit well defined for $\alpha\leq d$, where the divergence
$\mathcal{N}_{\alpha,L}\sim L^{d-\alpha}$ with the system size ensures that
energy scales extensively. For $\alpha>d$ the Kac rescaling factor saturates
to a finite value in the thermodynamic limit
$\mathcal{N}_{\alpha,L}\to\mathcal{N}_{\alpha}$.

Summary: We consider a $d$-dimensional quantum XY spin model with non-random
ferromagnetic interactions that decay algebraically with exponent $\alpha$.
The interaction strength is rescaled to make the mean field independent of
$\alpha$.

###  2.2 Equilibrium phase diagram in a nutshell

Quite generally ††margin:  $T=0$ phase diagram for ferromagnetic interactions
$J_{\mathbf{r},\mathbf{r}^{\prime}}\geq 0$ the system has an equilibrium zero-
temperature phase transition for small enough $\lvert h\rvert$, associated
with the spontaneous breaking of its $\mathbb{Z}_{2}$ spin-inversion symmetry
of the $x$-component. The longitudinal magnetization
$\braket{\hat{\sigma}_{\mathbb{r}}^{x}}$ undergoes an abrupt change from
$\braket{\hat{\sigma}_{\mathbb{r}}^{x}}=0$ in the unique paramagnetic ground
state for $\lvert h\rvert>h_{\text{cr}}$ to
$\braket{\hat{\sigma}_{\mathbb{r}}^{x}}_{\pm}=\pm m(h)\neq 0$ in the two
degenerate ferromagnetic ground states for $\lvert h\rvert<h_{\text{cr}}$.
When interactions are local, the universality class of this quantum phase
transition is the same as that of the $(d+1)$-dimensional classical Ising
model [48].

The possible emergence ††margin:  $T>0$ phase diagram of an ordered phase at
_finite temperature_ $T>0$, i.e. of long-range ordered excited states at
finite energy density, depends on the dimensionality $d$ and on the decay
exponent $\alpha$ of the interactions. For strictly finite range, finite-
temperature order is stable only for $d\geq 2$; in this case, the universal
properties of the thermal phase transition are the same as that of the
corresponding $d$-dimensional classical Ising model [48, 65]. Increasing the
interaction range, however, enhances the effective lattice connectivity,
somewhat similarly to the effect of increasing the lattice dimensionality [66,
67, 68]. While frustration prevents antiferromagnetic interactions from
creating collective ordering, ferromagnetic interactions do cooperate to
suppress the effect of spatial fluctuations, generally resulting in a
qualitative enhancement of the system’s ability to order as the interaction
range increases [69].

(a) Ising Model

(b) Finite temperature

Figure 1: Phase diagram of long-range interacting quantum Ising model. (a)
Mean-field theory properly describes the universal scaling behavior in the
cyan shaded region. The white background with the red $\rm LR$ label
represents the region in the phase diagram where the universal behavior is
correlated and influenced by the presence of relevant long-range couplings.
The boundary of the LR region is simply $\alpha=d+2$ in mean-field theory
(vertical shaded line), but gets displaced to $\alpha_{*}=d+2-\eta_{\rm sr}$
by two-loop corrections (red line). Finally, the white area on the right of
the red line signals the region of irrelevant long-range couplings where the
universal behavior is controlled by the local part of the interactions. Figure
reproduced from Ref. [70]. (b) Finite temperature phase-diagram of the one-
dimensional Ising model, in the transverse field $h$, interaction exponent
$\alpha$, and temperature $T$ space, in the units of $J$ (without Kac
normalization). For $\alpha<1$, the system is in the mean-field regime
(striped region), for $\alpha<5/3$ the mean-field universality is exact, while
for $\alpha>3$ the model is in the same universality class of the short-range
Ising model (dark grey). Figure adapted from Ref.[71].

The analogy ††margin:  Raise of effective dimension by $\alpha$ between
integer-dimension long-range systems and local systems in lower fractional
dimensions has been quantitatively tested in multiple studies in recent years,
both in classical [72, 68, 73, 74] and quantum [75, 70] long-range systems.
Leading-order perturbation theory results support the exact correspondence
between the universal behavior of long-range interacting systems with
dimension $d$ and decay exponent $\alpha$ and locally-interacting system with
dimension $d_{\rm eff}=2(d+z)/(\alpha-d)$, where $z$ is the dynamical critical
scaling exponent [75, 76, 77]. Advanced renormalization group studies
highlighted deviations from this correspondence, which only occur beyond the
leading order and hence remain small [70]. Therefore, it is possible to employ
the effective-dimension relation above to get the qualitative shape of the
phase diagram: For $\alpha<\frac{5}{3}d$ the universal scaling behavior is
captured by mean-field theory, while for $\frac{5}{3}d<\alpha<\alpha_{*}$ the
system displays correlated critical behavior influenced by the presence of
long-range interactions. Finally, for $\alpha>\alpha_{*}$ the interaction
tails become irrelevant and the critical exponents coincide with the ones of
the model with local interactions [70].

The location of $\alpha_{*}$ was subject to multiple controversies, but the
result $\alpha_{*}=d+2-\eta_{\rm sr}$ [72, 68], with $\eta_{\rm sr}$ the
anomalous dimension of the model with local interactions, appears now to be
established [70] in agreement with extensive numerical simulations in
classical models [73]. Due to the dependence of $\alpha_{*}$ on the universal
equilibrium properties of the local model, this boundary _does_ depend on the
particular model. For the model considered in this Report, the phase diagram
of the equilibrium critical problem is displayed in Fig. 1: in Fig. 1a we
report the universality properties for different dimension [70], while in Fig.
1b we show the phase diagram for finite temperature of the one-dimensional
model [71]. Decreasing the decay exponent below the dimension of the system,
i.e. $\alpha<d$, does not produce any major implications in the equilibrium
critical scaling, but it does modify the thermodynamic properties. In the
regime $\alpha<d$, the system becomes non-additive and the boundary
contribution to thermodynamic quantities cannot be neglected, leading to the
violation of several established equilibrium properties including the
equivalence of thermodynamic ensembles [78].

Given this scenario, long-range interacting systems can be classified in the
following way [46]: for $\alpha<d$ they are in the so-called _strong long-
range regime_ , for $d<\alpha<\alpha^{*}$ they are in the _weak long-range
regime_ , while for $\alpha>\alpha^{*}$ one retrieves short-range
properties.222We warn the readers that the nomenclature we adopt here is far
from being universally established in the vast literature on long-range
interactions.  In this Report we will mainly focus on quantum spin systems
around the strong long-range regime, i.e. $0\leq\alpha\lesssim d$, and we aim
at providing a cohesive picture for their distinctive dynamics.

In the light of the above discussion, the role of the spatial dimension is
diminished in systems with variable-range interactions, as the relevant
parameter in equilibrium is the effective dimension $d_{\rm eff}$. Yet, the
one-dimensional case is particularly interesting: In presence of local
interactions, $d=1$ systems cannot exhibit ordering at finite temperature,
because isolated topological defects of a ferromagnetically ordered pattern
(domain-wall-like excitations) cost a finite energy [69]; a longer range of
ferromagnetic interactions induces _binding_ between domain walls, and hence a
tendency to stabilize ferromagnetic order. The effect of long-range
interactions is thus most dramatic for $d=1$: The algebraically decaying
interactions in Eq. (1) allow to stabilize ferromagnetic order in the
thermodynamic limit upon decreasing $\alpha$ below $2$. This happens as the
interaction potential between two domain walls becomes _confining_ at large
distances, such that free isolated domain walls cost an infinite energy.

Summary: In equilibrium, the universal critical properties with $J/(\Delta
r)^{\alpha}$-interactions are close to those of the locally interacting
version of the system ($\alpha=\infty$) in a higher effective dimension
$d_{\rm eff}=2(d+z)/(\alpha-d)$. For $d=1$ and $\gamma\neq 0$, finite-
temperature ordering becomes possible for $\alpha\leq 2$.

###  2.3 Low-energy spectrum with infinite-range interactions ($\alpha=0$)

Let us start ††margin:  $\alpha=0$ by discussing the exactly solvable
infinite-range limit $\alpha\to 0$. This will be the starting point to analyze
the behavior for $\alpha>0$.

####  2.3.1 Mean-field theory as an exact classical limit

Increasing ††margin:  Collective spin the range of interactions $\alpha\to 0$
weakens spatial fluctuations, leading the system toward its mean-field limit —
similarly to the effect of increasing the system dimensionality $d\to\infty$.
This can be seen explicitly by rewriting the Hamiltonian (1) in terms of the
collective spin components

${\hat{S}}^{\mu}=\sum_{i=1}^{N}{\hat{s}}^{\mu}_{i}\,,\quad\mu=x,y,z\,,$ (4)

which gives the expression

${\hat{H}}_{\alpha=0}=-\frac{J_{0}}{N\,s^{2}}\left(\frac{1+\gamma}{2}(\hat{S}^{x})^{2}+\frac{1-\gamma}{2}(\hat{S}^{y})^{2}\right)-\frac{h}{s}\,\hat{S}^{z}\,\
,$ (5)

where we used $J=J_{0}/(N-1)\approx J_{0}/N$. This expression highlights that
the $\alpha=0$ Hamiltonian is a function of a single degree of freedom: the
collective spin. All other non-collective spin modes are frozen and do not
participate in dynamics. The collective spin magnitude
${{\hat{S}}}^{2}=(\hat{S}^{x})^{2}+(\hat{S}^{y})^{2}+(\hat{S}^{z})^{2}={S(S+1)}$
with $S=Ns,Ns-1,Ns-2,\dots,0$ or $1/2$ is conserved,

$\big{[}{\hat{S}}^{2},\hat{H}_{\alpha=0}\Big{]}=0\ .$ (6)

The Hilbert space sector $\mathcal{H}_{\hat{S}^{2}=S(S+1)}$ associated with
the quantum number $S$ contains $g_{N,S}$ copies of a spin-$S$ representation
of SU(2), where
$g_{N,S}=\rm{dim}\,\mathcal{H}_{\hat{S}^{z}=S}-\rm{dim}\,\mathcal{H}_{\hat{S}^{z}=S+1}$.
This combinatorial number depends implicitly on $s$; in the simplest case
$s=1/2$ we have

$g_{N,S}=\binom{N}{N/2-S}-\binom{N}{N/2-S-1}=\frac{2S+1}{N+1}\binom{N+1}{N/2-S}.$
(7)

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

For all states ††margin:  Classical limit with large $S$ growing with $N$, the
thermodynamic limit $N\to\infty$ is equivalent to a semiclassical limit for
the collective spin: The rescaled spin satisfies commutation relations of the
form

$\left[\frac{\hat{S}^{\mu}}{S},\frac{\hat{S}^{\nu}}{S}\right]=\frac{i}{S}\;\epsilon_{\mu\nu\rho}\;\frac{\hat{S}^{\rho}}{S}\
;$ (8)

and the Hamiltonian can we rewritten in terms of the rescaled spin as

${\hat{H}}_{\alpha=0}=(S/s)\bigg{\\{}-J_{0}\rho\left[\frac{1+\gamma}{2}\bigg{(}\frac{\hat{S}^{x}}{S}\bigg{)}^{2}+\frac{1-\gamma}{2}\bigg{(}\frac{\hat{S}^{y}}{S}\bigg{)}^{2}\right]-h\,\frac{\hat{S}^{z}}{S}\bigg{\\}}\,,$
(9)

where $\rho\equiv S/(Ns)$ is a constant depending on the collective spin
sector, with $0\leq\rho\leq 1$. Thus, the system manifestly has an effective
Planck constant $\hbar_{\text{eff}}\equiv 1/S$. Keeping in mind that a
meaningful thermodynamic limit requires to take $J_{0}\rho$ as a constant
independent of $N$, we conclude that the limit $N\to\infty$ realizes a
classical limit with a continuous spin

$\frac{\langle{\hat{\vec{S}}}\rangle}{S}\leadsto\vec{\mathcal{S}}$ (10)

of (conserved) length $1$ governed by the classical Hamiltonian
$\hat{H}_{\alpha=0}/(S/s)\leadsto\mathcal{H}_{\text{cl}}$,

$\mathcal{H}_{\text{cl}}(\vec{\mathcal{S}})=-{\rho
J_{0}}\left(\frac{1+\gamma}{2}(\mathcal{S}^{x})^{2}+\frac{1-\gamma}{2}(\mathcal{S}^{y})^{2}\right)-h\mathcal{S}^{z}\,,$
(11)

Canonical variables can be taken as, e.g., $\mathcal{S}^{z}=\cos\theta$ and
$\arctan_{2}(\mathcal{S}^{x},\mathcal{S}^{y})=\phi$.

The absolute ground state ††margin:  Ground state minimizes energy across all
sectors; for ferromagnetic interactions the ground state is realized for
maximal collective spin polarization, $S=Ns$, i.e. for $\rho=1$. A rigorous
implication of the classical limit [79] is that, as ${N\to\infty}$, the ground
state expectation values $\braket{\vec{\hat{S}}}_{\text{GS}}/S$ of the
collective spin components converge to the minimum point
$\vec{\mathcal{S}}^{*}$ of the classical Hamiltonian $\mathcal{H}_{\text{cl}}$
on the unit sphere. For later purpose it is convenient to define a rotated
reference frame $(\mathbf{X},\mathbf{Y},\mathbf{Z})$ adapted to the ground
state polarization, i.e., such that $\mathbf{Z}\equiv\vec{\mathcal{S}}^{*}$.
Using spherical coordinates we can parametrize

$\mathbb{X}\equiv\left(\begin{matrix}\cos\theta\cos\phi\\\
\cos\theta\sin\phi\\\
-\sin\theta\end{matrix}\right),\quad\mathbb{Y}\equiv\left(\begin{matrix}-\sin\phi\\\
\cos\phi\\\
0\end{matrix}\right),\quad\mathbb{Z}\equiv\left(\begin{matrix}\sin\theta\cos\phi\\\
\sin\theta\sin\phi\\\ \cos\theta\end{matrix}\right).$ (12)

Crucially, the quantum uncertainty
$\sqrt{\big{\langle}\big{(}\hat{S}^{X}\big{)}^{2}\big{\rangle}\big{\langle}\big{(}\hat{S}^{Y}\big{)}^{2}\big{\rangle}}$
associated with spin fluctuations in the transverse directions $\mathbf{X}$
and $\mathbf{Y}$ spans a phase-space area of order $h_{\rm eff}=2\pi/S$, which
is vanishingly small as $N\to\infty$.333This point will be further discussed
at length in Sec. 4.1.3.

The discussion above is valid for generic infinite-range Hamiltonians. For our
model in Eq. (5), minimization of $\mathcal{H}_{\rm cl}$ on the unit sphere
gives the ground-state polarization
$\vec{\mathcal{S}}^{*}=(\pm\sin\theta^{*},0,\cos\theta^{*})$, with

$\theta^{*}=\left\\{\begin{split}0&\quad\text{ for }h>h_{\rm cr}\equiv\rho
J_{0}(1+\gamma)\\\ \arccos\left(\frac{h}{\rho
J_{0}(1+\gamma)}\right)&\quad\text{ for }0\leq h\leq h_{\rm cr}\equiv\rho
J_{0}(1+\gamma)\end{split}\right.$ (13)

The ferromagnetic phase transition at $h=h_{\rm cr}$ is associated with the
bifurcation of the minimum. The ferromagnetic energy
$\mathcal{E}\equiv\mathcal{H}_{\rm cl}(\vec{\mathcal{S}}^{*})=-\rho
J_{0}(1+\gamma)/2-h^{2}/[\rho J_{0}(1+\gamma)]$ is extremal for $\rho=1$, in
agreement with the claim anticipated above and with intuition. See plots in
Figs. 2a and 2b.

Summary: The fully-connected Hamiltonian with $\alpha=0$ is a function of
collective spin variables only. The thermodynamic limit realizes a
semiclassical limit for the collective spin with an effective Planck constant
$\hbar_{\text{eff}}\propto 1/N$.

####  2.3.2 Collective quantum fluctuations and excitations

It is important ††margin:  Zero-point collective fluctuations to stress that,
in spite of the exact classical limit, the ground-state wavefunction is _not_
a product state of $N$ spins pointing in the direction
$\vec{\mathcal{S}}^{*}$: Collective interactions generate global
(_multipartite_) quantum entanglement among all spins. Such quantum
correlations stem from quantum fluctuations of the collective spin around the
average direction $\vec{\mathcal{S}}^{*}$. Such effects can be understood via
semiclassical analysis to leading order in $\hbar_{\rm eff}$.

Let us first compute the low-energy spectrum of the infinite-range Hamiltonian
(5) thought as a single-spin Hamiltonian, with $S$ growing with $N$. The
collective spin moves in an energy landscape whose depth grows with $N$. The
ground state wavefunction is localized around its global minimum (or minima).
Expansion of $\mathcal{H}_{\rm cl}$ around the minimum gives access to the
ground-state fluctuations and low-lying harmonic excitations. This can be
conveniently done via a Holstein-Primakoff transformation [80]: Recalling
$\hat{S}^{\pm}=\hat{S}^{x}\pm i\hat{S}^{y}$,

$\left\\{\begin{split}\hat{S}^{-}&=\hat{b}^{\dagger}\sqrt{2S-\hat{b}^{\dagger}\hat{b}}\,,\\\
\hat{S}^{+}&=\sqrt{2S-\hat{b}^{\dagger}\hat{b}}\;\;\hat{b}\,,\\\
\hat{S}^{z}&=S-\hat{b}^{\dagger}\hat{b}\,.\end{split}\right.$ (14)

These equations represent an exact embedding of a quantum spin into a bosonic
mode.

This procedure is simplest in the paramagnetic phase.††margin:  Single-spin
spectrum in the symmetric phase For large $h\gg h_{\rm cr}$ the ground state
approaches the uncorrelated state fully polarized along $z$, and the
elementary excitations approach the tower of spin lowering excitations. For
finite $h>h_{\rm cr}$ the collective spin fluctuates along the transverse
directions — more prominently along the “soft” direction $x$ and more weakly
along the “stiff” direction $y$.444Recall that we assumed $\gamma>0$ for
definiteness. Such fluctuations can be described by mapping $\hat{S}^{x}$ and
$\hat{S}^{y}$ to canonical bosonic operators via Eq. (14),

$\left\\{\begin{split}\hat{S}^{x}&\approx\sqrt{S}\;\hat{q}\,,\\\
\hat{S}^{y}&\approx\sqrt{S}\;\hat{p}\,,\\\
\hat{S}^{z}&=S-\hat{n}_{0}=S-\frac{\hat{q}^{2}+\hat{p}^{2}-1}{2}\,.\end{split}\right.$
(15)

Using $[\hat{q},\hat{p}]=i$ one can check that for large $S$ the spin
commutation relations are satisfied by the right-hand sides of Eqs. (15) to
leading order. In a classical phase-space description, the approximation given
by the above truncated Holstein-Primakoff transformation corresponds to
replacing the surface of the sphere by its tangent plane at the North pole.

Using Eq. (15), the Hamiltonian (5) can be approximated by neglecting terms of
order $1/S$, and hence easily diagonalized. We find:

$\displaystyle\hat{H}_{\alpha=0}$ $\displaystyle\approx-N\rho
h+\frac{h}{s}\frac{\hat{q}^{2}+\hat{p}^{2}-1}{2}-\frac{\rho
J_{0}}{s}\bigg{(}\frac{1+\gamma}{2}\hat{q}^{2}+\frac{1-\gamma}{2}\hat{p}^{2}\bigg{)}$
(16a) $\displaystyle=-N\rho
h+\frac{1}{s}\bigg{(}\frac{\omega_{>}-\omega_{>}^{(0)}}{2}\bigg{)}+\frac{1}{s}\omega_{>}\;\hat{n}\
,$ (16b)

where

$\omega_{>}=\sqrt{[h-\rho J_{0}(1-\gamma)][h-\rho
J_{0}(1+\gamma)]},\qquad\omega_{>}^{(0)}=h\ .$ (17)

The first term in the last line of Eq. (16b) represents the classical energy,
and the second one is the variation of the zero-point energy due to quantum
fluctuations around the classical configuration. In the last term, $\hat{n}$
is the harmonic excitation quanta of energy $\omega_{>}$ (not to be confused
with the “bare” spin-lowering excitation quanta $\hat{n}_{0}$).

For ††margin:  Phase transition $h>h_{\text{cr}}$, the number
$\braket{\hat{n}_{0}}=\braket{\hat{q}^{2}+\hat{p}^{2}-1}/2$ of bare collective
spin excitations in the ground state is finite, and it diverges as $h\searrow
h_{\text{cr}}$, signaling a critical phenomenon (see Fig. 2d). Indeed the
energy gap $\omega_{>}/s$ closes at $h=h_{\text{cr}}$, with a mean-field
critical exponent $1/2$ (see Fig. 2c). For $h<h_{\text{cr}}$ the frequency
$\omega_{>}$ becomes imaginary, which signals instability of the paramagnetic
state.

In ††margin:  Single-spin spectrum in the broken-symmetry phase order to
determine the ground state and the elementary excitations in the broken-
symmetry phase, let us start from some general considerations. For
$h<h_{\text{cr}}$ the classical landscape presents two symmetric minima, as
discussed above. Below the energy
$E_{\text{dyn}}\equiv\mathcal{H}_{\text{cl}}(\theta=0)$ of the classical
phase-space separatrix, two symmetric families of classical trajectories fill
the two energy wells. In the thermodynamic limit, this corresponds to two
towers of pairwise degenerate energy levels, associated with wavefunctions
localized in the two wells. At finite size $N$, however, the energy
eigenstates below the critical energy are nondegenerate and alternately even
and odd with respect to the $\mathbb{Z}_{2}$ symmetry of the Hamiltonian. For
large $N$, they approach even and odd superpositions of the localized
wavefunctions. The energy splitting between each pair of quasidegenerate
eigenstates is proportional to the quantum tunneling amplitude across the
energy barrier, which is exponentially small in the height of the barrier
[81], and hence exponentially small in $N$. Accordingly, tunneling between the
two broken-symmetry sectors is practically suppressed even for moderate system
sizes, and it is extremely fragile to tiny symmetry-breaking perturbations.
For these reasons it makes sense to consider the two towers of symmetry-
breaking states independently of each other.

(a) Order parameter

(b) Ground-state energy density

(c) Excitation gap

(d) Polarization depletion

(e) Entanglement entropy

Figure 2: Equilibrium properties of the fully-connected quantum Ising model
across the phase diagram. All quantities are computed for the Hamiltonian (5)
with $\gamma=1$, and display singularities at the phase transition $h=2J_{0}$
with mean-field critical exponents. (The behavior of the corresponding
quantities for a general anisotropic XY model $\gamma\neq 0$ is analogous.)
Panel (a): The order parameter is determined by Eq. (13). Panel (b): The
ground-state energy density is the minimum of the classical energy (11) on the
unit sphere. Panel (c): The dark-blue curve corresponds to $\omega_{<}$ in Eq.
(22) and $\omega_{>}$ in Eq. (17), respectively below and above the phase
transition; the light-blue curve corresponds to $\omega_{\rm sw}$, cf. Eq.
(24). Panel (d): The number of collective spin excitations is computed by
diagonalizing Eqs. (16a) and (20), respectively below and above the phase
transition. Panel (e): The bipartite (half-system) entanglement entropy in
fully-connected spin models is a function of $\langle\hat{n}_{0}\rangle$ only,
see Sec. 4.1.5 below; in the quantum ferromagnetic phase, the finite-size
ground state approaches a symmetric superposition of two symmetry-breaking
ground states for large $N$, which yields an extra bit of entropy; the
divergence at criticality is logarithmic in system size, see Refs. [82, 83].

To compute the spectrum explicitly, it is convenient to introduce a procedure
which will lend itself to powerful generalizations in the rest of this Report.
We rewrite the components of the collective spin in a rotated frame
$(\mathbb{X},\mathbb{Y},\mathbb{Z})$, cf. Eq. (12), by angle $\theta$ in the
$xz$-plane, i.e.,

$\hat{S}^{x}=\cos\theta\,\hat{S}^{X}+\sin\theta\,\hat{S}^{Z},\qquad\hat{S}^{y}=\hat{S}^{Y},\qquad\hat{S}^{z}=-\sin\theta\,\hat{S}^{X}+\cos\theta\,\hat{S}^{Z}\,.$
(18)

Performing a Holstein-Primakoff transformation with rotated quantization axis
$\mathbb{Z}$ and neglecting terms of order $1/S$,

$\left\\{\begin{split}\hat{S}^{X}&\approx\sqrt{S}\;\hat{q},\\\
\hat{S}^{Y}&\approx\sqrt{S}\;\hat{p},\\\
\hat{S}^{Z}&=S-\hat{n}_{0}=S-\frac{\hat{q}^{2}+\hat{p}^{2}-1}{2},\end{split}\right.$
(19)

we get

$\begin{split}\hat{H}_{\alpha=0}\approx&-N\left(h\rho\cos\theta+J_{0}\frac{1+\gamma}{2}\rho^{2}\sin^{2}\theta\right)\\\
&-\sqrt{\frac{\rho N}{s}}\sin\theta\Big{(}h-\rho
J_{0}(1+\gamma)\cos\theta\Big{)}\hat{q}\\\ &+\frac{1}{s}\bigg{[}\left(\rho
J_{0}(1+\gamma)\sin^{2}\theta+h\cos\theta\right)\frac{\hat{q}^{2}+\hat{p}^{2}-1}{2}-\rho
J_{0}(1+\gamma)\cos^{2}\theta\,\frac{\hat{q}^{2}}{2}-\rho
J_{0}(1-\gamma)\,\frac{\hat{p}^{2}}{2}\,\bigg{]}.\end{split}$ (20)

In order for the bosonic variables to describe quantum fluctuations it is
necessary to align the frame with the classical configuration, in such a way
that linear terms in the second line vanish. This condition leads to
$\theta^{*}$ as in Eq. (13). The resulting quadratic Hamiltonian can then be
readily diagonalized:

$\hat{H}_{\alpha=0}\approx-\frac{N}{2}\bigg{(}\frac{h^{2}}{J_{0}(1+\gamma)}+\rho^{2}J_{0}(1+\gamma)\bigg{)}+\frac{1}{s}\bigg{(}\frac{\omega_{<}-\omega_{<}^{(0)}}{2}\bigg{)}+\frac{1}{s}\omega_{<}\;\hat{n},$
(21)

where

$\omega_{<}=\sqrt{\left[\rho^{2}J_{0}^{2}(1+\gamma)^{2}-h^{2}\right]\frac{2\gamma}{1+\gamma}},\qquad\omega_{<}^{(0)}=\rho
J_{0}(1+\gamma)\,.$ (22)

Analogously to Eq. (16), the first term in the last line of Eq. (21)
represents the classical energy, the second one expresses the shift in the
zero-point energy due to quantum fluctuations around the classical minimum
configuration, while the last one [arising from diagonalization of Eq. (20)]
is the energy of the harmonic excitations, with $n=0,1,2,\dots$.

In Fig. 2 ††margin:  Validity of semiclassical analysis we plotted the exact
ground state energy density $\mathcal{E}$ [panel (b)], the energy gap of
collective spin excitations $\omega_{>,<}$ (dark-blue curve) [panel (c)], and
the number of “bare” collective spin excitations $\langle\hat{n}_{0}\rangle$
[panel (d)], of the infinite-range quantum Ising model ($\gamma=1$) in the
thermodynamic limit $N\to\infty$, as a function of the ratio $h/J_{0}$. The
results (16) and (21) are asymptotically exact for $n/N\to 0$, and _fully
nonperturbative_ in the Hamiltonian parameters $h,J_{0},\gamma$. Systematic
improvements in powers of $n/N$ can be worked out with a more refined analysis
[84]. This is particularly relevant to understand the finite-size scaling
$\omega\sim N^{-1/3}$ of the energy gap at criticality $h=h_{\text{cr}}$: see
A for an elementary semiclassical derivation.

(For completeness, Fig. 2e also reports the ground-state bipartite
entanglement entropy across the phase diagram. This quantity can be computed
numerically for large $N$ [82] and compared with analytical calculations in
the large-$N$ limit based on semiclassical fluctuations [83]. This analytical
procedure can be deduced as particular case of the more general discussion on
entanglement dynamics in Sec. 4.1.5 below; for this reason, we do not discuss
this here.)

Summary: The collective spin low-energy spectrum is described by bosonic
excitations, obtained by a Holstein-Primakoff expansion around the classical
ground state.

####  2.3.3 “Spin-wave” excitations

The analysis ††margin:  Complete many-body spectrum above concerns collective
spin quantum fluctuations and excitations within a fixed sector with
collective spin length $S$ — and we are ultimately interested in the ground
state sector with maximal $S=Ns$. Different families of spin excitations lower
the collective spin length to $S=Ns-n_{\rm sw}$, with $n_{\rm
sw}=0,1,2,\dots$. (For reasons that will become clear below, we will refer to
the quantum number $n_{\rm sw}$ as the _total occupation of spin-wave modes
with non-vanishing momenta_.) Their spectrum can also be straightforwardly
obtained from semiclassical arguments: Recalling the definition $\rho=S/(Ns)$
above, we have

$\rho=1-\frac{n_{\text{sw}}}{Ns}.$ (23)

Substituting into Eqs. (16) and (21) and consistently neglecting terms of
higher order in $1/N$, we obtain the complete spectrum of low-lying
excitations above the ground state to leading order in $n/N$ and
$n_{\text{sw}}/N$:

$\begin{split}H_{\alpha=0}&\approx-
Nh+\frac{\omega_{>}-\omega_{>}^{(0)}}{2s}+\frac{1}{s}\big{(}\omega_{>}\,\hat{n}+h\,\hat{n}_{\text{sw}}\big{)},\\\
H_{\alpha=0}&\approx-\frac{N}{2}\bigg{(}\frac{h^{2}}{(1+\gamma)J_{0}}+J_{0}(1+\gamma)\bigg{)}+\frac{\omega_{<}-\omega_{<}^{(0)}}{2s}+\frac{1}{s}\big{(}\omega_{<}\,\hat{n}+(1+\gamma)J_{0}\,\hat{n}_{\text{sw}}\big{)},\end{split}$
(24)

valid for $h>(1+\gamma)J_{0}$ and $h<(1+\gamma)J_{0}$, respectively. (Here
$\omega$’s are taken at $\rho=1$.) In Fig. 2c we additionally reported the
“spin-wave” excitation gap $\omega_{\rm sw}=h$ or $J_{0}(1+\gamma)$ in the two
phases. Note that the Hilbert space sector dimension grows exponentially with
$n_{\rm sw}$ [cf. the exact expression in Eq. (7)]; however, because of
permutational invariance, these energy levels are exactly degenerate.

As discussed so far, the properties of infinite-range spin Hamiltonians can be
efficiently computed either analytically (via a large-$N$ asymptotic
expansion) or numerically (via exact diagonalization of the single-spin
problem for $S\leq N/2\approx 10^{5}$). In closing this Subsection it is worth
to briefly mention that the Hamiltonian (5) is equivalent to the Lipkin-
Meshkov-Glick model of nuclear physics [85, 86, 87], which is actually Bethe-
ansatz solvable [88]; however, this solution is not practically useful for
large $N$, and semiclassical or numerical techniques give much easier access
to the relevant information.

Summary: “Spin-wave” excitations — lowering the collective spin length —
remain gapped and dispersionless throughout the phase diagram for $\alpha=0$.

###  2.4 Finite-range interactions ($\alpha>0$)

The tendency ††margin:  Parameter $\alpha$ as a perturbation of long-range
interactions to form collective spin alignment and to preserve it even in
excited states becomes increasingly prominent as $\alpha$ is decreased. To
quantify this aspect, it is convenient to view a long-range interacting system
as a “perturbation” of the infinite-range interacting system with all-to-all
interactions ($\alpha=0$).

####  2.4.1 Perturbation to mean-field

This viewpoint can be made explicit by rewriting the Hamiltonian in momentum
space. To this aim, we Fourier transform the spin operators
$\hat{s}^{\mu}_{\mathbf{r}}$ for $\mu=x,y,z$:

$\tilde{S}^{\mu}_{\mathbf{k}}=\sum_{\mathbf{r}}e^{i\mathbf{k}\cdot\mathbf{r}}\hat{s}^{\mu}_{\mathbf{r}}\,,$
(25)

with

$\mathbf{k}\equiv\mathbf{k}_{\boldsymbol{\ell}}=2\pi\boldsymbol{\ell}/L,\quad\boldsymbol{\ell}=(\ell_{1},\dots,\ell_{d}),\quad\ell_{a}=0,\pm
1,\pm 2,\dots,\pm\lfloor L/2\rfloor$ (26)

(for $L$ even $\ell_{a}={\pm L/2}$ coincide). We also define
$\tilde{S}^{\pm}_{\mathbf{k}}=\tilde{S}^{x}_{\mathbf{k}}\pm
i\tilde{S}^{y}_{\mathbf{k}}$. Note that

$\tilde{\vec{S}}_{\mathbf{k}=\mathbf{0}}\equiv\hat{\vec{S}}=\sum_{\mathbf{r}}\hat{\vec{s}}_{\mathbf{r}}$
(27)

is the system’s collective spin. It is straightforward to separate the
variable-range quantum XY Hamiltonian (1) into the $\alpha$-independent
collective part — given by the $\mathbf{k}=\mathbf{0}$ terms — and the
“perturbation” controlled by $\alpha$:

$\hat{H}_{\alpha}=\hat{H}_{\alpha=0}+\hat{V}_{\alpha}$ (28)

with555Note that in this expression the various ${\mathbf{k}}$-modes are _not_
dynamically decoupled, since
$\big{[}\tilde{S}^{\mu}_{\mathbf{k}},\tilde{S}^{\nu}_{\mathbf{q}}\big{]}=i\epsilon^{\mu\nu\lambda}\tilde{S}^{\lambda}_{{\mathbf{k}+\mathbf{q}}}$.

$\hat{V}_{\alpha}=-\frac{J_{0}}{4s^{2}N}\sum_{\mathbf{k}\neq\mathbf{0}}f_{\mathbf{k}}(\alpha)\Big{[}\big{(}\tilde{S}^{+}_{\mathbf{k}}\tilde{S}^{-}_{-\mathbf{k}}+\tilde{S}^{-}_{\mathbf{k}}\tilde{S}^{+}_{-\mathbf{k}}\big{)}+\gamma\big{(}\tilde{S}^{+}_{\mathbf{k}}\tilde{S}^{+}_{-\mathbf{k}}+\tilde{S}^{-}_{\mathbf{k}}\tilde{S}^{-}_{-\mathbf{k}}\big{)}\Big{]}\,.$
(29)

In Eq. (29) we defined the function

$f_{\mathbf{k}}(\alpha)=\sum_{\mathbf{r}\neq\mathbf{0}}\frac{\cos(\mathbf{k}\cdot\mathbf{r})}{\lvert\lvert\mathbf{r}\rvert\rvert^{\alpha}}\Bigg{/}\sum_{\mathbf{r}\neq\mathbf{0}}\frac{1}{\lvert\lvert\mathbf{r}\rvert\rvert^{\alpha}}$
(30)

which depends implicitly on the dimensionality $d$ of the lattice. By
construction, $f_{\mathbf{k}=\mathbf{0}}(\alpha)=1$.

Figure 3: Plots of the function ${f}_{\alpha,k}$ (30) for $d=1$. (Left panel):
${f}_{k}(\alpha)$ is shown for several values of $\alpha$, for $N=L=500$. The
function squeezes towards $k=0$ for $0\leq\alpha\leq 1$. For $1<\alpha<2$,
${f}_{k}(\alpha)$ becomes a finite function with a cusp behavior for small
$k$, while for $\alpha\gg 2$ it is a cosine-like function. (Central panel):
${f}_{k}(\alpha)$ is shown for $\alpha=0.7$ and increasing values of $N$.
Qualitatively similar behavior occurs for $0\leq\alpha\leq 1$. Squeezing
towards a delta function as $N\to\infty$ occurs with a speed $N^{-(1-\alpha)}$
for $\alpha<1$ and $1/\ln N$ for $\alpha=1$. (Right panel): a “zoom” of the
plot in the bottom left panel is shown, for larger values of $N$. The rescaled
function in the vicinity of $k=0$ converges to a finite limiting curve as
$N\to\infty$. This _discrete_ structure approaches a continuum as
$\alpha\nearrow 1$.

When $\alpha\to 0$ ††margin:  Properties of $f_{\mathbf{k}}(\alpha)$ the
couplings $f_{\mathbf{k}\neq\mathbf{0}}(\alpha)$ turn off [Eq. (33)], and
$\hat{H}_{\alpha}$ reduces to a Hamiltonian describing a single collective
degree of freedom. The effect of spatially modulated interactions $\alpha\neq
0$ is then to couple the collective spin to all finite-wavelength modes
describing spatially non-trivial spin fluctuations, resulting in complex
interacting many-body dynamics. The form of this coupling is dictated by the
function $f_{\mathbf{k}}(\alpha)$, which plays a crucial role in the physics
of long-range interacting systems. In C we derive the following asymptotic
estimates:

$\displaystyle{f}_{\mathbf{k}_{\boldsymbol{\ell}}\neq\mathbf{0}}(\alpha)\equiv{f}_{\boldsymbol{\ell}\neq\mathbf{0}}(\alpha)$
$\displaystyle\underset{|\boldsymbol{\ell}|\to\infty}{\thicksim}$
$\displaystyle\frac{A(\alpha)}{|\boldsymbol{\ell}|^{d-\alpha}}+\frac{B(\alpha)}{|\boldsymbol{\ell}|^{(d+1)/2}}\,\qquad\qquad\text{for
}0<\alpha<d\,;$ (31) $\displaystyle{f}_{\mathbf{k}\neq\mathbf{0}}(\alpha)$
$\displaystyle\underset{\mathbf{k}\to\mathbf{0}}{\thicksim}$ $\displaystyle
1-\bar{A}(\alpha)|\mathbf{k}|^{\alpha-d}-\bar{B}(\alpha)|\mathbf{k}|^{2}\quad\text{for
}\alpha>d\,.$ (32)

The sharp changes in behavior are summarized in Fig. 3, where we plot
${f}_{k}(\alpha)$ for a range of values of $\alpha$ and $d=1$. Its shape
shrinks from $f_{k}(\alpha\to\infty)=\cos k$ to

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

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

$\displaystyle f_{k_{\ell}}(\alpha)\equiv f_{\ell}(\alpha)$ $\displaystyle\sim
c(\alpha)|\ell|^{-(1-\alpha)}$ $\displaystyle\quad\text{for }\alpha<1;$ (34a)
$\displaystyle f_{k}(\alpha)$ $\displaystyle\sim 1-c(\alpha)|k|^{\alpha-1}$
$\displaystyle\quad\text{for }1<\alpha<3;$ (34b) $\displaystyle f_{k}(\alpha)$
$\displaystyle\sim 1-c(\alpha)k^{2}$ $\displaystyle\quad\text{for }\alpha>3.$
(34c)

For long-range interactions $0<\alpha<1$, the values of $f_{k}(\alpha)$
progressively squeeze onto the vertical axis as $L\to\infty$; upon zooming
near $k=0$ one finds a sequence of _discrete_ finite values, see the right
panel of Fig. 3 [58, 52]. This phenomenon can be physically interpreted as
follows: interactions decay so slowly with the spatial distance that the
system behaves as a permutationally invariant system over finite length
scales, hence observables are unable to resolve finite wavelengths. Only modes
with extensive wavelengths $k_{\ell}\propto 1/L$ may impact the physical
properties. As $\alpha$ is increased to values larger than $d$, all modes
$k\neq 0$ get eventually activated.

Despite its simplicity, ††margin:  Discrete spectrum for $0<\alpha<d$ the
result in Eq. (31) has significant physical implications: As we will show
below, the low-energy spectrum of a quantum system with long-range
interactions remains discrete in the thermodynamic limit. In this precise
sense we may say that _long-range interacting systems with $0<\alpha<d$
interpolate between few-body and many-body physics_.

At the same time, ††margin:  Divergent velocity for $\alpha<d+1$ Eq. (32)
showcases another fundamental properties of long-range interacting systems:
the singularity at small momenta gives rise to a divergent velocity of
propagation of quantum information across the system for $\alpha<d+1$,
violating the famous Lieb-Robinson light-cone bound of short-range interacting
systems [12]. This property is actually completely general, as it does not
rely on any low-energy description. We will further discuss its consequences
in Sec. 3.1.1.

Summary: Finite-range interactions can be seen as a perturbation to the fully-
connected Hamiltonian. Physically, this perturbation couples the collective
spin to spin fluctuation modes at non-vanishing momentum. Long-range
interactions preferentially generate coupling to long-wavelength modes only.

####  2.4.2 Quantum paramagnetic phase

We are ††margin:  Bosonization of paramagnetic spin fluctuations now ready to
compute the low-energy spectrum and properties of the variable-range quantum
XY model (1). For large $h$ and arbitrary $\alpha$ the paramagnetic ground
state is the spin-coherent eigenstate of $S^{z}$,

$|{\rm
GS}_{h=\infty}\rangle=|\boldsymbol{\Uparrow}\,\rangle\equiv\bigotimes_{\mathbf{r}}|\uparrow_{\mathbf{r}}\rangle.$
(35)

For finite $h/J_{0}$ the ground state is a ($\alpha$-dependent) distortion of
this state dressed by spin-lowering excitations. A convenient approach to
describe spin fluctuations in $|{\rm GS}_{h<\infty}\rangle$ is by mapping them
to bosonic modes using the Holstein-Primakoff transformation [80].

Unlike our discussion in Sec. 2.3.2, we now have to keep track of spatially-
resolved fluctuating spin modes, which we can conveniently do by working at
the level of individual microscopic spins. Recalling the standard definitions
$\hat{s}^{\pm}_{\mathbf{r}}=\hat{s}^{x}_{\mathbf{r}}\pm
i\hat{s}^{y}_{\mathbf{r}}$, we can bosonize individual spin fluctuations
around the positive $\mathbf{z}$ axis by setting

$\hat{s}^{-}_{\mathbf{r}}\mapsto\hat{b}^{\dagger}_{\mathbf{r}}\sqrt{2s-\hat{b}^{\dagger}_{\mathbf{r}}\hat{b}_{\mathbf{r}}},\quad\hat{s}^{+}_{\mathbf{r}}\mapsto\sqrt{2s-\hat{b}^{\dagger}_{\mathbf{r}}\hat{b}_{\mathbf{r}}}\;\hat{b}_{\mathbf{r}},\quad\hat{s}^{z}_{\mathbf{r}}\mapsto
s-\hat{b}^{\dagger}_{\mathbf{r}}\hat{b}_{\mathbf{r}},$ (36)

where $\hat{b}_{\mathbf{r}},\hat{b}^{\dagger}_{\mathbf{r}}$ are canonical
bosonic annihilation and creation operators. Performing the substitutions (36)
in Eq. (1) we obtain an exact representation of the variable-range XY quantum
spin model as a non-linear bosonic Hamiltonian, where the state
$|\boldsymbol{\Uparrow}\,\rangle$ in Eq. (35) corresponds to the Fock space
vacuum $|\emptyset\rangle$.

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. The states $|\uparrow\,\rangle$ and $|\downarrow\,\rangle$ are
mapped onto $|0\rangle\equiv|\emptyset\rangle$ and $|1\rangle\equiv
b^{\dagger}|\emptyset\rangle$. The operators on the right-hand sides of Eqs.
(36) act non-trivially on the full bosonic space; however, they are block-
diagonal, as their matrix elements between the physical spin subspace and its
orthogonal complement are vanishing; their action on the physical spin
subspace coincides with the operators on the left-hand sides.

It is convenient ††margin:  Bosonic representation of XY quantum spin model to
write the bosonic Hamiltonian directly in momentum space. To this aim, we
define the Fourier-transformed bosonic modes666Note that we take a unitary
Fourier transformation on the bosonic modes, while the convention for spins in
Eq. (25) was such that
$\tilde{S}^{x,y,z}_{\mathbf{k}=\mathbf{0}}=\hat{S}^{x,y,z}$ (collective spin
projections).

$\tilde{b}^{\dagger}_{\mathbf{k}}=\frac{1}{\sqrt{N}}\sum_{\mathbf{r}}e^{i\mathbf{k}\cdot\mathbf{r}}\hat{b}^{\dagger}_{\mathbf{r}}\,.$
(37)

We now formally expand the Holstein-Primakoff mapping (36) in $1/s$ and
Fourier-transform term by term:

$\left\\{\begin{split}\tilde{S}^{-}_{\mathbf{k}}&\approx(2Ns)^{1/2}\,\tilde{b}^{\dagger}_{\mathbf{k}}-\frac{1}{2(2Ns)^{1/2}}\sum_{\mathbf{q}_{1},\mathbf{q}_{2}}\tilde{b}^{\dagger}_{\mathbf{q}_{1}}\tilde{b}^{\dagger}_{\mathbf{q}_{2}}\tilde{b}_{\mathbf{q}_{1}+\mathbf{q}_{2}-\mathbf{k}},\\\
\tilde{S}^{+}_{\mathbf{k}}&\approx(2Ns)^{1/2}\,\tilde{b}_{-\mathbf{k}}-\frac{1}{2(2Ns)^{1/2}}\sum_{\mathbf{q}_{1},\mathbf{q}_{2}}\tilde{b}^{\dagger}_{\mathbf{q}_{1}+\mathbf{q}_{2}+\mathbf{k}}\tilde{b}_{\mathbf{q}_{1}}\tilde{b}_{\mathbf{q}_{2}},\\\
\tilde{S}^{z}_{\mathbf{k}}&=Ns\;\delta_{\mathbf{k},0}-\sum_{\mathbf{q}}\tilde{b}^{\dagger}_{\mathbf{q}+\mathbf{k}}\tilde{b}_{\mathbf{q}}.\end{split}\right.$
(38)

It is worth to stress here the connection with the previously introduced
expansion. First of all, we immediately recognize that the bosonic mode with
$\mathbf{k}=\mathbf{0}$ coincides with the previously introduced collective
bosonic mode in Eq. (19). Furthermore, by expanding $\hat{S}^{2}$ using Eqs.
(38), one can check that $\hat{n}_{\mathbf{k}=\mathbf{0}}$ cancels to leading
order [80]:

$\hat{n}_{\rm sw}\equiv
Ns-S=\sum_{\mathbf{k}\neq\mathbf{0}}\tilde{b}^{\dagger}_{\mathbf{k}}\tilde{b}_{\mathbf{k}}.$
(39)

This substantiates the naming of this quantity introduced before Eq. (23).

Making the substitutions (38) into Eq. (28) we obtain an expression of the
form

$\hat{H}_{\alpha}=\frac{1}{s}\bigg{[}(Ns)^{1}\mathcal{E}_{0}+(Ns)^{0}\hat{H}_{2}+(Ns)^{-1}\hat{H}_{4}+\dots\,\bigg{]},$
(40)

where:

$\mathcal{E}_{0}=\mathcal{H}(\mathbf{z})=-h$ (41)

is the classical (mean-field) energy density of the paramagnetic state;

$\hat{H}_{2}=\sum_{\mathbf{k}}h\,\tilde{b}^{\dagger}_{\mathbf{k}}\tilde{b}_{\mathbf{k}}-\sum_{\mathbf{k}}{J_{0}}\,f_{\mathbf{k}}(\alpha)\left(\frac{\tilde{b}_{\mathbf{k}}\tilde{b}^{\dagger}_{\mathbf{k}}+\tilde{b}^{\dagger}_{-\mathbf{k}}\tilde{b}_{-\mathbf{k}}}{2}+\gamma\frac{\tilde{b}_{\mathbf{k}}\tilde{b}_{-\mathbf{k}}+\tilde{b}^{\dagger}_{-\mathbf{k}}\tilde{b}^{\dagger}_{\mathbf{k}}}{2}\right)$
(42)

describes semiclassical (Gaussian) spin fluctuations;

$\hat{H}_{4}=\frac{J_{0}}{2}\sum_{\mathbf{k},\mathbf{q}_{1},\mathbf{q}_{2}}f_{\mathbf{k}}(\alpha)\;\times\\\
\Big{[}\Big{(}\tilde{b}^{\dagger}_{-\mathbf{k}+\mathbf{q}_{1}+\mathbf{q}_{2}}\tilde{b}_{\mathbf{q}_{1}}\tilde{b}_{\mathbf{q}_{2}}\tilde{b}^{\dagger}_{\mathbf{k}}+\tilde{b}^{\dagger}_{\mathbf{q}_{1}}\tilde{b}^{\dagger}_{\mathbf{q}_{2}}\tilde{b}_{\mathbf{k}+\mathbf{q}_{1}+\mathbf{q}_{2}}\tilde{b}_{-\mathbf{k}}+\tilde{b}_{\mathbf{k}}\tilde{b}^{\dagger}_{\mathbf{q}_{1}}\tilde{b}^{\dagger}_{\mathbf{q}_{2}}\tilde{b}_{-\mathbf{k}+\mathbf{q}_{1}+\mathbf{q}_{2}}+\tilde{b}^{\dagger}_{-\mathbf{k}}\tilde{b}^{\dagger}_{\mathbf{k}+\mathbf{q}_{1}+\mathbf{q}_{2}}\tilde{b}_{\mathbf{q}_{1}}\tilde{b}_{\mathbf{q}_{2}}\Big{)}\\\
+\gamma\Big{(}\tilde{b}^{\dagger}_{-\mathbf{k}+\mathbf{q}_{1}+\mathbf{q}_{2}}\tilde{b}_{\mathbf{q}_{1}}\tilde{b}_{\mathbf{q}_{2}}\tilde{b}_{-\mathbf{k}}+\tilde{b}_{\mathbf{k}}\tilde{b}^{\dagger}_{\mathbf{k}+\mathbf{q}_{1}+\mathbf{q}_{2}}\tilde{b}_{\mathbf{q}_{1}}\tilde{b}_{\mathbf{q}_{2}}+\tilde{b}^{\dagger}_{-\mathbf{k}}\tilde{b}^{\dagger}_{\mathbf{q}_{1}}\tilde{b}^{\dagger}_{\mathbf{q}_{2}}\tilde{b}_{-\mathbf{k}+\mathbf{q}_{1}+\mathbf{q}_{2}}+\tilde{b}^{\dagger}_{\mathbf{q}_{1}}\tilde{b}^{\dagger}_{\mathbf{q}_{2}}\tilde{b}_{\mathbf{k}+\mathbf{q}_{1}+\mathbf{q}_{2}}\tilde{b}^{\dagger}_{\mathbf{k}}\Big{)}\Big{]}$
(43)

represents the 2-body non-linear interactions between spin fluctuations. One
can similarly derive $(Ns)^{-2}\hat{H}_{6}$ etc.

While ††margin:  Linear spin-wave theory the full exact bosonic representation
is cumbersome, its usefulness rests on the approximability of highly polarized
spin states with bosonic states. To this aim we introduce the number of bosons

$\hat{n}_{\rm tot}\equiv\hat{n}_{0}+\hat{n}_{\rm sw}=Ns-\hat{S}^{z}$ (44)

and we approximate well-polarized states with $n_{\rm tot}\ll Ns$ with dilute
Fock states with $n_{\rm tot}$ boson. In such corner of the Hilbert space, the
bosonic modes turn out to provide an accurate description of spin states and
operators. Intuitively, by inspecting Eqs. (38), one recognizes that the
action of the non-linear terms (second on the right-hand side) on a dilute
Fock state is suppressed by a density factor $n_{\rm tot}/N$ compared to the
action of the leading terms. Thus, up to an error of order
$\mathcal{O}[(n_{\rm tot}/N)^{2}]$, we may identify
$\tilde{S}^{-}_{\mathbf{k}}\propto\tilde{b}^{\dagger}_{\mathbf{k}}$,
$\tilde{S}^{+}_{\mathbf{k}}\propto\tilde{b}_{-\mathbf{k}}$. We now show that
the ground state of long-range interacting spin models lives exactly in this
corner of the spin space.

An approximate solution of the bosonic Hamiltonian (40) can be found by
neglecting the terms with $\hat{H}_{4}$ and higher order — an approximation
usually termed linear spin-wave (LSW) theory. The quality of the result
heavily depends on the parameters and in particular on $\alpha$. Our purpose
is to show that the LSW description of low-energy properties becomes _exact_
for $\alpha<d$, and quantify its accuracy for $\alpha>d$.

The quadratic spin-wave Hamiltonian can be diagonalized via a standard
Bogolubov transformation,
$\tilde{b}_{\mathbf{k}}=\cosh\theta_{\mathbf{k}}\beta_{\mathbf{k}}+\sinh\theta_{\mathbf{k}}\beta^{\dagger}_{-\mathbf{k}}$,
with

$\tanh(2\theta_{\mathbf{k}})\equiv\frac{\gamma
J_{0}f_{\mathbf{k}}(\alpha)}{h-J_{0}f_{\mathbf{k}}(\alpha)}\,.$ (45)

The result is

$Ns\,\mathcal{E}_{0}+\hat{H}_{2}=Ns\,\mathcal{E}_{2}+\sum_{\mathbf{k}}\omega_{\mathbf{k},>}(\alpha)\,\hat{\beta}^{\dagger}_{\mathbf{k}}\hat{\beta}_{\mathbf{k}}\
,$ (46)

where we identify the excitation spectrum

$\omega_{\mathbf{k},>}(\alpha)=\sqrt{\big{[}h-J_{0}f_{\mathbf{k}}(\alpha)\big{]}^{2}-\gamma^{2}\big{[}J_{0}f_{\mathbf{k}}(\alpha)\big{]}^{2}}=\sqrt{\big{[}h-J_{0}f_{\mathbf{k}}(\alpha)(1-\gamma)\big{]}\big{[}h-J_{0}f_{\mathbf{k}}(\alpha)(1+\gamma)\big{]}},$
(47)

and the ground-state energy

$Ns\,\mathcal{E}_{2}=Ns\,\mathcal{E}_{0}+\frac{1}{2}\sum_{\mathbf{k}}\big{[}\omega_{\mathbf{k},>}(\alpha)-\omega_{>}^{(0)}\big{]}$
(48)

where $\omega_{>}^{(0)}=h$ [cf. Eq. (17)].

Within LSW theory, the ground-state wavefunction is given by

$\begin{split}|{\rm
GS}_{2}\rangle&=\prod_{\mathbf{k}}\exp\bigg{[}\frac{\theta_{\mathbf{k}}}{2}\big{(}\tilde{b}_{\mathbf{k}}\tilde{b}_{-\mathbf{k}}-\tilde{b}^{\dagger}_{-\mathbf{k}}\tilde{b}^{\dagger}_{\mathbf{k}}\big{)}\bigg{]}\,|\emptyset\rangle\propto\prod_{\mathbf{k}}\exp\bigg{(}-\frac{\epsilon_{\mathbf{k}}}{4\gamma
J_{0}f_{\mathbf{k}}(\alpha)}\tilde{b}^{\dagger}_{-\mathbf{k}}\tilde{b}^{\dagger}_{\mathbf{k}}\bigg{)}\,|\emptyset\rangle\,\end{split}$
(49)

where

$\epsilon_{\mathbf{k}}\equiv
2h-2J_{0}f_{\mathbf{k}}(\alpha)-\omega_{\mathbf{k},>}(\alpha)\geq 0\,.$ (50)

Meaningfulness of the LSW solution is determined by $\omega_{\mathbf{k},>}$
being real. This requires $h\geq h_{\rm cr}$, where

$h_{\rm cr}\equiv J_{0}(1+\gamma).$ (51)

The minimum of $\omega_{\mathbf{k},>}$ is attained as
$\mathbf{k}\to\mathbf{0}$. To expand around this limit we write
$f_{\mathbf{k}}(\alpha)\equiv 1-\sigma_{\mathbf{k}}(\alpha)$. Calculation
gives

$\omega_{\mathbf{k},>}(\alpha)\underset{\mathbf{k}\to\mathbf{0}}{\thicksim}2h\sqrt{a+b\sigma_{\mathbf{k}}(\alpha)},$
(52)

with dimensionless coefficients $a$ and $b$.777Explicitly,
$a=1-2J_{0}/h+(1-\gamma^{2})(J_{0}/h)^{2}$ and
$b=2(J_{0}/h)[1-(1-\gamma^{2})(J_{0}/h)]$. For $h>h_{\rm cr}$ one has $a>0$
and thus
$\omega_{\mathbf{k},>}\underset{\mathbf{k}\to\mathbf{0}}{\thicksim}2h\sqrt{a}+h(b/\sqrt{a})\sigma_{\mathbf{k}}$.
For short-range interactions $\alpha\geq 3$ the gapped dispersion relation is
parabolic, $\sigma_{\mathbf{k}}\sim|\mathbf{k}|^{2}$; for longer range
$d<\alpha<d+2$ it behaves as $\sigma_{\mathbf{k}}\sim|\mathbf{k}|^{\alpha-d}$;
for $\alpha<d$ the spectrum becomes discrete,
$\sigma_{\mathbf{k}_{\boldsymbol{\ell}}}\equiv\sigma_{\boldsymbol{\ell}}$.888It
is interesting to note that the LSW description of the ground state is
meaningful even for large $\alpha$, where LSW theory completely fails to
capture the possible topological nature of excitations. At the critical point
$h=h_{\rm cr}$, one has $a=0$ and hence
$\omega_{\mathbf{k},>}\underset{\mathbf{k}\to\mathbf{0}}{\thicksim}2h\sqrt{b\sigma_{\mathbf{k}}}$,
signaling closure of the spectral gap at $\mathbf{k}=\mathbf{0}$. However, for
$\alpha<d$, the spectrum of spin-wave excitations with
$\mathbf{k}\neq\mathbf{0}$ is discrete.

To assess ††margin:  Self-consistency of LSW theory the accuracy of LSW theory
we evaluate the depletion of spin polarization, i.e.

$\langle\hat{n}_{\rm
tot}\rangle=Ns-\langle\hat{S}^{z}\rangle=\sum_{\mathbf{k}}\langle{\rm
GS}|\tilde{b}^{\dagger}_{\mathbf{k}}\tilde{b}_{\mathbf{k}}|{\rm GS}\rangle.$
(53)

Approximating the ground-state $|{\rm GS}\rangle$ by the LSW theory ground-
state $|{\rm GS}_{2}\rangle$ in Eq. (49) we obtain the explicit expression

$\langle\hat{n}_{\rm
tot}\rangle=\frac{1}{2}\sum_{\mathbf{k}}\frac{\epsilon_{\mathbf{k}}}{\omega_{\mathbf{k},>}}\,.$
(54)

This quantity depends on $h$, $\alpha$ and $\gamma$; in particular, it is
suppressed as $h\to\infty$ or $\alpha\to 0$ or $\gamma\to 0$.

Figure 4:  Ground-state spin depletion density $\langle\hat{n}_{\rm
tot}\rangle/N$ in the quantum paramagnetic phase, cf. Eq. (54), for $\gamma=1$
and $d=1$ (variable-range quantum Ising chain).

In Fig. 4 we plot the depletion per spin $\langle\hat{n}_{\rm tot}\rangle/N$
given by Eq. (54) at fixed $\gamma=1$ (quantum Ising model) and for $d=1$. As
is evident, the effect of spin fluctuations is enhanced as the interactions
become short-ranged, i.e. $\alpha\to\infty$, or as the critical point $h_{\rm
cr}=2J_{0}$ is approached. All the qualitative aspects of this plot can be
understood analytically. Specifically, for $h>h_{\rm cr}$, we have

$\frac{\langle\hat{n}_{\rm tot}\rangle}{Ns}=\frac{1}{16}\bigg{(}\frac{h_{\rm
cr}}{h}\bigg{)}^{2}\frac{1}{Ns}\sum_{\mathbf{k}}f^{2}_{\mathbf{k}}(\alpha)+\mathcal{O}\Bigg{(}\bigg{(}\frac{h_{\rm
cr}}{h}\bigg{)}^{3}\frac{1}{Ns}\sum_{\mathbf{k}}f^{3}_{\mathbf{k}}(\alpha)\Bigg{)}.$
(55)

The behavior of the right-hand side as $N\to\infty$ depends qualitatively on
$\alpha$: For $\alpha>d$ the limit is a finite number,

$\frac{1}{N}\sum_{\mathbf{k}}f^{2}_{\mathbf{k}}(\alpha)\sim\int_{-\pi}^{\pi}\dots\int_{-\pi}^{\pi}\frac{dk_{1}\dots
dk_{d}}{(2\pi)^{d}}f^{2}_{\mathbf{k}}(\alpha)$ (56)

(cf. left panel of Fig. 3). As $\alpha\searrow d$ the function
$f_{\mathbf{k}}(\alpha)$ squeezes on the vertical axis (cf. Fig. 3),
suppressing the value of the integral. This means that the spin depletion
becomes subextensive for $\alpha<d$: using
$f_{\boldsymbol{\ell}}(\alpha)\sim|\boldsymbol{\ell}|^{-(1-\alpha)}$ [cf. Eq.
(34)], one finds

$\langle\hat{n}_{\rm
tot}\rangle\sim\sum_{|\boldsymbol{\ell}|<L/2}|f_{\boldsymbol{\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.$ (57)

On the other hand, for $h=h_{\rm cr}$, we have

$\frac{\langle\hat{n}_{\rm
tot}\rangle}{Ns}=\frac{1}{4Ns}\sum_{\mathbf{k}}\Bigg{(}\frac{1+\sigma_{\mathbf{k}}(\alpha)}{\sqrt{\sigma_{\mathbf{k}}(\alpha)}}-1\Bigg{)}.$
(58)

Here the behavior of the right-hand side as $L\to\infty$ depends even more
strongly on $\alpha$, in particular for one-dimensional systems $d=1$: For
$\alpha>3$ the sum is divergent,
$\frac{1}{N}\sum_{k}\frac{1}{\sqrt{\sigma_{k}(\alpha)}}\sim\int_{-\pi}^{\pi}\frac{dk}{2\pi}\frac{1}{|k|}=\infty$.
This divergence witnesses the inadequacy of LSW theory to describe critical
behavior of one-dimensional systems with short-range interactions. Contrarily,
for $1<\alpha<3$ one has $\sigma_{k}(\alpha)\sim|k|^{\alpha-1}$, and the
integral is convergent. As $\alpha\searrow 1$ the depletion per spin is
suppressed, making LSW theory increasingly accurate. Finally, for $\alpha<1$,
one finds the same subextensive scaling as in Eq. (57). Note, however, that
the collective spin mode $k=0$ yields an additional divergent (but still
subextensive) contribution $\langle\hat{n}_{0}\rangle\sim N^{1/3}$ to
$\langle\hat{n}_{\rm tot}\rangle$ at the critical point, which can be shown by
semiclassical analysis (see A); such a contribution is thus dominant for
$0<\alpha<\frac{2}{3}d$ and subleading for $\alpha>\frac{2}{3}d$.

The bottom line of this Section is that the Holstein-Primakoff description of
spin fluctuations is _exact_ in the thermodynamic limit for $0<\alpha<d$, and
otherwise an increasingly accurate approximation as $\alpha$ is decreased
towards $d$. Importantly, this result is true _regardless_ of the value of
$s$, down to $s=1/2$. Accuracy for low $s$ may be surprising at first sight,
considering Eq. (36). Its origin can be traced back to the observation that
the truncated Holstein-Primakoff mapping gives exact matrix elements within
the subspace with at most one boson on each site, for arbitrary $s$. Thus,
what really controls the quality of the approximation is the ground-state
spin-wave _density_ : For $0<\alpha<d$ the probability of finding two or more
bosons in a given site in $|{\rm GS}_{2}\rangle$ is vanishingly small in the
thermodynamic limit, and it is finite but parametrically small for
$\alpha\gtrsim d$.

Summary: The paramagnetic ground state can be determined via linear spin-wave
theory, which becomes exact as the strong long-range regime is approached.

####  2.4.3 Quantum ferromagnetic phase

To derive ††margin:  Bosonic representation of ferromagnetic spin fluctuations
the low-energy spectrum in the quantum ferromagnetic phase for $\alpha>0$, we
promote the frame rotation in Eq. (18) from the level of the collective spin
to the level of individual spins:

$\hat{s}^{x}_{\mathbf{r}}=\cos\theta\,\hat{s}^{X}_{\mathbf{r}}+\sin\theta\,\hat{s}^{Z}_{\mathbf{r}},\qquad\hat{s}^{y}_{\mathbf{r}}=\hat{s}^{Y}_{\mathbf{r}},\qquad\hat{s}^{z}_{\mathbf{r}}=-\sin\theta\,\hat{s}^{X}_{\mathbf{r}}+\cos\theta\,\hat{s}^{Z}_{\mathbf{r}}\,.$
(59)

Hence we perform a Holstein-Primakoff expansion of individual spins with
quantization axis $\mathbf{Z}$ and Fourier-transform,

$\left\\{\begin{split}\tilde{S}^{X}_{\mathbf{k}}&\approx(Ns)^{1/2}\,\frac{\tilde{b}^{\dagger}_{\mathbf{k}}+\tilde{b}_{-\mathbf{k}}}{\sqrt{2}}\,,\\\
\tilde{S}^{Y}_{\mathbf{k}}&\approx(Ns)^{1/2}\,\frac{\tilde{b}^{\dagger}_{\mathbf{k}}-\tilde{b}_{-\mathbf{k}}}{\sqrt{2}i}\,,\\\
\tilde{S}^{Z}_{\mathbf{k}}&=Ns\;\delta_{\mathbf{k},0}-\sum_{\mathbf{q}}\tilde{b}^{\dagger}_{\mathbf{q}+\mathbf{k}}\tilde{b}_{\mathbf{q}}\,,\end{split}\right.$
(60)

and substitute into the Hamiltonian (28). As in Eq. (40) we obtain a formal
series in inverse powers of $Ns$, including the classical energy
$(Ns)^{1}\mathcal{H}(\mathbf{Z})$, the quadratic bosonic Hamiltonian

$\hat{H}_{2}=\left(h\cos\theta+J_{0}(1+\gamma)\sin^{2}\theta\right)\sum_{\mathbf{k}}\tilde{b}^{\dagger}_{\mathbf{k}}\tilde{b}_{\mathbf{k}}\\\
-J_{0}\sum_{\mathbf{k}}f_{\mathbf{k}}(\alpha)\bigg{[}\left(\frac{1+\gamma}{2}\cos^{2}\theta+\frac{1-\gamma}{2}\right)\frac{\tilde{b}_{\mathbf{k}}\tilde{b}^{\dagger}_{\mathbf{k}}+\tilde{b}^{\dagger}_{-\mathbf{k}}\tilde{b}_{-\mathbf{k}}}{2}\\\
+\left(\frac{1+\gamma}{2}\cos^{2}\theta-\frac{1-\gamma}{2}\right)\frac{\tilde{b}_{\mathbf{k}}\tilde{b}_{-\mathbf{k}}+\tilde{b}^{\dagger}_{-\mathbf{k}}\tilde{b}^{\dagger}_{\mathbf{k}}}{2}\bigg{]}\,,$
(61)

as well as quartic and higher-order interactions involving an even number of
bosons. However, unlike in Eq. (40), we also get an additional term
$(Ns)^{1/2}\hat{H}_{1}$ linear in
$\hat{q}_{\mathbf{k}=\mathbf{0}}=(\hat{b}_{\mathbf{k}=\mathbf{0}}+\hat{b}^{\dagger}_{\mathbf{k}=\mathbf{0}})/\sqrt{2}$
— cf. the first line of Eq. (20) — as well as other odd terms
$(Ns)^{-1/2}\hat{H}_{3}$ and so on.

The rotation angle $\theta^{*}$ must be determined by imposing that the
expectation value of $\hat{q}_{\mathbf{k}=\mathbf{0}}$ vanishes. To lowest
order this gives the mean-field solution in Eq. (13). Paralleling the
derivation in the previous Section we can then solve the LSW Hamiltonian
$\hat{H}_{2}(\theta=\theta^{*})$, which yields the spectrum

$\omega_{\mathbf{k},<}(\alpha)=\sqrt{\left[J_{0}^{2}(1+\gamma)^{2}-h^{2}f_{\mathbf{k}}(\alpha)\right]\left[1-f_{\mathbf{k}}(\alpha)\frac{1-\gamma}{1+\gamma}\right]}$
(62)

as well as the zero-point energy shift
$\frac{1}{2}\sum_{\mathbf{k}}(\omega_{\mathbf{k},<}(\alpha)-\omega_{<}^{(0)})$,
where $\omega_{<}^{(0)}=J_{0}(1+\gamma)$ [cf. Eq. (22)].

The analysis ††margin:  Self-consistency of LSW theory and beyond of spectral
properties and of the spin depletion in the quantum paramagnetic phase can be
repeated for the quantum ferromagnetic phase, with qualitatively similar
conclusions. The mean-field description of local observables is _exact_ for
$0<\alpha<d$ in the thermodynamic limit. For $\alpha>d$ finite corrections to
the mean-field results arise. Such corrections can be evaluated within the
bosonic formalism [57, 89]. In particular, the downward shift of the quantum
critical point $h_{\rm cr,\alpha}=h_{\rm cr,\alpha=0}-\delta h_{\rm
cr,\alpha}$ due to quantum fluctuations amounts to [57, 89]

$\frac{\delta h_{\rm cr,\alpha}}{h_{\rm
cr,\alpha=0}}\;=\;\frac{\gamma}{s}\,\frac{2\
+3\gamma}{4(1+\gamma)^{2}}\,\int_{-\pi}^{\pi}\dots\int_{-\pi}^{\pi}\frac{dk_{1}\dots
dk_{d}}{(2\pi)^{d}}f_{\mathbf{k}}^{2}(\alpha)\,.$ (63)

The right-hand side is in fact vanishing for $0<\alpha<d$ and grows finite for
$\alpha>d$. Note that effects of quantum fluctuations are suppressed as
$s\to\infty$.

For completeness, let us mention that for $\alpha<d$, the ground state shares
the same basic properties of the fully-connected limit [90]. Long-range
interactions however can induce unexpected entanglement properties. For
instance, for $d<\alpha<d+1$, the ground state entanglement entropy the long-
range Dyson Hierarchical model obeys an area law at criticality [91], due to
its special Tree Tensor Network structure [92]. On the other hand, numerical
studies for antiferromagnetic long-range systems have shown violations of
area-law scaling also in the gapped phase [93, 94, 95, 96].

Summary: The ferromagnetic ground state can be determined via linear spin-wave
theory in a rotated frame. This approach is exact in the strong long-range
regime and it determines corrections to the location of quantum critical
point.

###  2.5 Structure of the spectrum beyond linear spin-wave theory

In this final Subsection we comment on the structure of the many-body low-
energy spectrum beyond LSW theory. To grasp such effects we will make use of
degenerate perturbation theory — i.e., of the Schrieffer-Wolff transformation
— around points in the two phases where $\hat{H}_{\alpha}$ becomes diagonal.

Quite generally,††margin:  Beyond LSW in paramagnetic phase spin waves provide
a rather complete description of low energy properties in the quantum
paramagnetic phase, even beyond LSW theory. This is best understood in the
regime of large external field $h$, where the ground state is $|{\rm
GS}\rangle=|\boldsymbol{\Uparrow}\rangle$ and excited states can be described
as a set of individual spin lowering excitations. The degeneracy of blocks
with multiple spin excitations is split by the interactions. The effective
block Hamiltonian for large $h$ is obtained by projecting out interaction
terms that do not conserve $S^{z}$, i.e. the part of $\hat{H}_{\alpha}$
proportional to $\gamma$.999E.g., it can be checked that $\gamma$ drops out
from the LSW spectrum $\omega_{\mathbf{k},>}(\alpha)$ in Eq. (47) to lowest
order in $1/h$.  The Holstein-Primakoff transformation maps this effective XX
quantum spin model to a model of lattice bosons with variable-range hopping
$\propto 1/(\Delta r)^{\alpha}$. Non-linearity of the mapping is associated
with the suppression of matrix elements for hopping processes from/to multiply
occupied sites.101010In particular, transition amplitudes to states with more
than $2s$ bosons at any site are strictly vanishing. Such unconventional
multi-boson interactions [cf. the second line of Eq. (43)] produce a
scattering phase shift for quantized spin-wave excitations, which can be
determined by solving the few-body problem. A qualitatively similar scenario
is expected for finite but large enough $h/J_{0}$, upon rotating the bare
Holstein-Primakoff bosons $b_{\mathbf{k}},b^{\dagger}_{\mathbf{k}}$ to the
dressed spin-wave basis $\beta_{\mathbf{k}},\beta^{\dagger}_{\mathbf{k}}$ in
$\hat{H}_{\alpha}$.

On the other hand,††margin:  Beyond LSW in ferromagnetic phase the
phenomenology is drastically different in the quantum ferromagnetic phase, due
to the strong binding tendency of local spin excitations. This is most easily
understood starting from the “classical” Ising model, i.e. Eq. (1) with
$\gamma=1$, $h=0$:

$\hat{H}_{\alpha}=-\sum_{\mathbf{r},\mathbf{r^{\prime}}}J_{\lvert\lvert\mathbf{r}-\mathbf{r^{\prime}}\rvert\rvert}(\alpha)\,\hat{\sigma}^{x}_{\mathbf{r}}\hat{\sigma}^{x}_{\mathbf{r^{\prime}}}\,.$
(64)

This Hamiltonian is diagonal in the $x$-basis. It has two degenerate ground
states $|{\rm GS}_{+}\rangle=|\boldsymbol{\Leftarrow}\rangle$ and $|{\rm
GS}_{-}\rangle=|\boldsymbol{\Rightarrow}\rangle$ with energy $E_{\rm
GS}=-J_{0}N$, and excited states can be described as a set of spin-changing
excitations with respect to either ground state. Here the unperturbed energy
levels depend on the details of the spin configuration: The lowest elementary
excitations are individual isolated spin excitations,

$E\Big{(}s^{x}_{\mathbf{r}}=s-1,\,s^{x}_{\mathbf{R}\neq\mathbf{r}}=s\Big{)}=E_{\rm
GS}+\frac{J_{0}}{s}\,;$ (65)

the subspace with two spin excitations has a configuration-dependent energy,

$E\Big{(}s^{x}_{\mathbf{r}}=s-1,\,s^{x}_{\mathbf{r}^{\prime}}=s-1,\,s^{x}_{\mathbf{R}\neq\mathbf{r},\mathbf{r}^{\prime}}=s\Big{)}=E_{\rm
GS}+2\frac{J_{0}}{s}-\frac{J_{\mathbf{r},\mathbf{r^{\prime}}}(\alpha)}{s^{2}}\,.$
(66)
$E\Big{(}s^{x}_{\mathbf{r}}=s-2,\,s^{x}_{\mathbf{R}\neq\mathbf{r}}=s\Big{)}=E_{\rm
GS}+2\frac{J_{0}}{s}\qquad\text{(for $s>1/2$ only)}\,.$ (67)

This implies an attractive potential between spin excitations,111111Note that
spin excitations on the same site (relevant for $s>1/2$ only) do not feel any
attraction or repulsion, unless the Hamiltonian features self-interaction
terms.

$V_{\mathbf{r},\mathbf{r^{\prime}}}(\alpha)=-\frac{J_{\mathbf{r},\mathbf{r^{\prime}}}(\alpha)}{s^{2}}=-\frac{J_{0}}{2s^{2}\mathcal{N}_{\alpha,L}}\,\frac{1}{\lvert\lvert\mathbf{r}-\mathbf{r^{\prime}}\rvert\rvert^{\alpha}}\,.$
(68)

Similarly one can compute the unperturbed energy of more complex spin
configurations with three or more spin excitations.

Such excited states acquire a non-trivial dispersion relation upon turning on
$h\neq 0$ or $\gamma\neq 1$. Using lowest-order degenerate perturbation
theory, it is straightforward to check that processes generated by
$1-\gamma\neq 0$ induce a variable-range hopping of individual spin
excitations, whereas processes generated by $h\neq 0$ do not induce any
resonant transitions to lowest order. We thus retrieve a dispersion relation
$\sim\frac{J_{0}}{s}[1-f_{\mathbf{k}}(\alpha)(1-\gamma)/2]$ for individual
spin excitations, in agreement with Eq. (62) from LSW theory. The number
$N_{b}$ of stable spin-wave bound states depends on the relative magnitude of
“classical potential well depths”, controlled by $\alpha$, and “quantum
hopping amplitudes”, controlled by $1-\gamma$ and $h$. This number grows
unbounded upon reducing the quantum fluctuations. Estimating $N_{b}$ as well
as the lifetime of unstable bound states depending on the interaction range
and in one or higher dimensions is in general a challenging problem which, to
the best of our knowledge, has not been discussed extensively; see however
Refs. [62, 97].

All the observations above ††margin:  Short-range limit carry over to short-
range interacting systems $\alpha=\infty$, _provided_ the system
dimensionality is large enough ($d\geq 2$).

In one dimension ††margin:  1d: Domain walls and confinement LSW is still a
meaningful description of the paramagnetic spectrum (asymptotically exact for
large external field). In the quantum ferromagnetic phase, however, LSW theory
completely misses the relevant degrees of freedom, i.e. topological domain-
wall-like excitations. In the simplest case $s=1/2$ these fractionalized spin
excitations can be described as fermions, as the exact solution of the XY
quantum spin chain [98] makes manifest. The qualitative effect of longer-range
interactions is then to create an effective attractive potential $v_{\Delta
r}(\alpha)$ between domain walls at a distance $\Delta r$ (NB not to be
confused with the attractive potential $V_{\Delta r}(\alpha)$ between
individual spin excitations introduced above!). Taking for simplicity the
classical Ising limit $\gamma=1$, $h=0$ as a reference, one can
straightforwardly compute the excess energy of a spin configuration with two
domain walls separated by $\Delta r$ sites: Assuming $\alpha>1$,121212For
$0<\alpha<1$ we have $v_{\Delta r}(\alpha)=2J_{0}\Delta r$ in the
thermodynamic limit, in agreement with naive LSW theory. In this case,
however, the spatial configuration of the $\Delta r$ flipped spins becomes
immaterial. Thus, it is not meaningful to speak about “domain-wall
confinement”.

$v_{\Delta r}(\alpha)=\frac{2J_{0}}{\zeta(\alpha)}\bigg{(}\sum_{r=1}^{\Delta
r}r\frac{1}{r^{\alpha}}+r\sum_{r={\Delta
r}+1}^{\infty}\frac{1}{r^{\alpha}}\bigg{)}\,.$ (69)

This potential grows from $v_{\Delta r=1}(\alpha)=2J_{0}$ to $v_{\Delta
r=\infty}(\alpha)=2J_{0}\zeta(\alpha-1)/\zeta(\alpha)$ for $\alpha>2$, or to
$\infty$ for $\alpha\leq 2$. The asymptotic behavior of $V_{\Delta r}(\alpha)$
at large distance $\Delta r$ is

$v_{\Delta r}(\alpha)\;\underset{\Delta
r\to\infty}{\thicksim}\;\frac{2J_{0}}{\zeta(\alpha)}\times\left\\{\begin{array}[]{ll}\displaystyle\frac{\Delta
r^{2-\alpha}}{(2-\alpha)(\alpha-1)}&\text{for }1<\alpha<2,\vspace*{1mm}\\\
\log\Delta r&\text{for }\alpha=2,\\\ \displaystyle\zeta(\alpha-1)-\frac{\Delta
r^{-(\alpha-2)}}{(\alpha-2)(\alpha-1)}&\text{for }\alpha>2.\end{array}\right.$
(70)

For $\alpha>2$ finitely many bound states coexist with unbound deconfined
domain walls. As anticipated above, the cost of having a deconfined domain
wall blows up as $\alpha\searrow 2$, which witnesses the stabilization of
long-range order by long-range interactions in 1d. Upon decreasing $\alpha$
the lowest excitation in the spectrum — the tightest bound state between two
domain walls — is increasingly well described by LSW theory.

We finally note that while long-range interacting quantum spin chains do not
naturally map to local lattice gauge theories [99], except in special cases
[100], the spatial confinement, the spatial confinement of domain walls bears
qualitative resemblance with quark confinement in high energy physics [101].
This bridge led to insights on the anomalous non-equilibrium dynamics of these
systems [62, 63, 102]. Furthermore, although domain-wall deconfinement
prevents finite-temperature ordering for $\alpha>2$, it has been shown that
the existence of low-lying bound states is associated with a severe
suppression of the dynamical melting rate of the order parameter after shallow
quantum quenches [103].

Summary: The low-energy spectrum in the quantum ferromagnetic phase hosts a
rich structure of spin-wave bound states for $d\geq 2$, or for $d=1$ provided
$\alpha$ is low enough. In $d=1$ deconfined domain-wall-like excitations
appear for $\alpha>2$ along with confined spin-wave-like bound states.

##  3 Low-energy dynamics

The previous Section shows how the main impact of long-range interactions on
low-energy equilibrium properties of the Hamiltonian in Eq. (1) can be traced
back to the quasi-particle spectrum, in turn determined by the function
$f_{\mathbf{k}}(\alpha)$, see Fig. 3. Interestingly, the same is true for
several types of non-equilibrium phenomena at low energies. In the present
Section, we focus on dynamics following weak perturbations of the ground
state, which is captured by the quadratic quasi-particle Hamiltonian (such as
Eq. (42)). This allows us to capture the dynamics of quantum correlations at
low energies [50, 51, 104, 41], see Sec. 3.1, the appearance of long-live
metastable state [42, 52], i.e. the quasi-stationary states (QSSs), see Sec.
3.2, the universal defect formation upon slowly traversing criticality [105,
106, 107, 53, 54], see Sec. 3.3, and the appearance of dynamical quantum phase
transitions in the Loschmidt echo (DQPTs) [108, 109], see Sec. 3.4. Since Sec.
3.1 and 3.2 focus on super-critical quenches, the dynamics occurs in the near
equilibrium regime, where the spin-wave expansion around the equilibrium state
remains applicable. On the other hand, universal defect scaling and DQPTs are
observed for quenches across the critical point, making the applicability of
the low-energy theory a priori questionable. Nevertheless, we are going to
show how the salient features of those critical quenches actually arise from a
low density of excitations above the ground state.

###  3.1 Spreading of correlations

††margin:  Lieb-Robinson bound

In systems governed by local Hamiltonians, out-of-equilibrium quantum
correlations are known to spread within a “light cone”: The propagation of
information in non-relativistic quantum lattice systems with bounded local
Hilbert space obeys a speed limit given by the _Lieb-Robinson theorem_ [110].
This states that the support of an operator $\hat{A}_{\mathbb{r}}$ initially
localized in a finite region around site $\mathbb{r}$ and evolving in the
Heisenberg representation with a local Hamiltonian $\hat{H}$, spreads in space
with a finite (model-dependent) velocity $v_{\text{LR}}$. Formally, for any
locally interacting lattice system there exist positive constants $\xi,\mu$
and $v_{\text{LR}}$ such that the correlator between two operators at distance
$\Delta r=|\mathbb{r}^{\prime}-\mathbb{r}|$

$||[\hat{A}_{\mathbb{r}}(t),\hat{B}_{\mathbb{r}^{\prime}}(0)]||\leq\xi||\hat{A}_{\mathbb{r}}||\;||\hat{B}_{\mathbb{r}^{\prime}}||\;e^{-\mu\,\max\left(0,\Delta
r-v_{\text{LR}}t\right)}\ ,$ (71)

where $||\cdot||$ is the operator norm. Namely, the weight of the time-evolved
operator outside the “light-cone” region $t\geq\Delta r\,/v_{\text{LR}}$ is
exponentially suppressed with $\Delta r-v_{\text{LR}}t\to\infty$. In other
words, it takes at least a time proportional to the distance $t\propto\Delta
r$ to send information at a distance $\Delta r$. Such light-cone propagation
of information is by now theoretically well understood in short-range
interacting systems, and it goes hand-in-hand with a linear dynamical increase
of bipartite entanglement out of equilibrium [111, 112, 113, 114].
Experimental observation of linear light-cone propagation [15, 115] has been
accompanied by abundant numerical confirmations [116, 117, 118, 119].

In presence of long-range interactions, the standard behavior of locally
interacting systems changes substantially: the bounds on the group velocity
may not hold anymore, and the spreading of correlations and information
scrambling may be drastically boosted.

Figure 5:  Spatial spreading of correlations in systems with power-law
interactions. (a) Connected correlation function in a long-range trapped ion
platform following a global quench with $\alpha\approx 0.64$. Image adapted
from Ref. [3]. (b) Violation of the Lieb-Robinson bound in Eq. (71) for long-
range interacting systems with power-law interactions for $\alpha>d$. Image
adapted from Ref. [120].

The study of the impact of algebraically decaying interactions on correlation
spreading has been addressed as a function of the different values of the
power-law exponent $\alpha$. Part of the current understanding is based on
assessing the behavior of the spatial spreading of connected correlations,
e.g.

$G_{\alpha\beta}(r,t)=\langle\hat{\sigma}^{\alpha}_{i+r}(t)\hat{\sigma}^{\beta}_{i}(0)\rangle-\langle\hat{\sigma}^{\alpha}_{i+r}(t)\rangle\langle\hat{\sigma}^{\beta}_{i}(0)\rangle\
,$ (72)

in paradigmatic quantum spin chains or tight-binding models; the expectation
value is taken over some initial state $\ket{\psi_{0}}$ and
$\alpha,\beta=x,y,z$.

Generalized bounds have been derived for long-range systems [121, 122], see
Ref.[123] for a recent comprehensive review. The related experiments and
numerical investigations have, however, led to conflicting pictures [59, 124,
14, 3, 60, 50, 125]. For instance, experiments on ion chains [3] and numerical
simulations within truncated Wigner approximation [126] for the one-
dimensional long-range XY model point towards bounded, super-ballistic,
propagation for all values of $\alpha$. In contrast, experiments on the long-
range transverse Ising model reported ballistic propagation of correlation
maxima with, however, observable leaks that increase when $\alpha$ decreases
[14]. Moreover, time-dependent density matrix renormalization group (t-DMRG)
and variational Monte-Carlo (t-VMC) numerical simulations indicate the
existence of three distinct regimes, namely instantaneous, sub-ballistic, and
ballistic, for increasing values of the exponent $\alpha$, see Ref. [127, 59,
124, 60, 50, 125].
In the following, we will see how these difficulties can be overcome in the
restricted setting of near-equilibrium dynamics, by studying correlation
spreading within linear spin-wave theory.

Summary: The Lieb-Robinson bound forbids super-ballistic spreading of quantum
correlations in locally interacting systems. Long-range interactions allow to
circumvent this constraint.

####  3.1.1 Weak long-range regime ($\alpha>d$)

Let us first consider the case of the Ising Hamiltonian, i.e. Eq. (1) with
$\gamma=1$, but restrict our study to the spin-wave representation in Eq.
(46). In this Section, we aim at characterizing the universal scaling of
correlations following Ref. [51]. Let us simplify the spin-wave dispersion
relation in Eq. (47) by considering its low-momentum asymptotic expression,

$\displaystyle\omega_{\mathbf{k}}\mapsto\omega^{\rm
low}_{\mathbf{k}}=\Delta+c{k}^{\zeta},$ (73)

where the gap $\Delta=\sqrt{h\left(h+2J_{0}f_{0}(\alpha)\right)}$ is finite,
$c=\sqrt{\frac{h}{h+2J_{0}f_{0}(\alpha)}}J_{0}\frac{\partial
f_{0}(\alpha)}{\partial{k}}$, and $\zeta=\alpha-d$. As long as $\alpha>d$ the
quasi-particle energy remains finite, while the group velocity diverges for
$d<\alpha<d+1$. ††margin:  Quenches within the paramagnetic phase The system
is prepared in its ground state and the coupling is suddenly quenched from
$J_{0}^{i}\to J_{0}^{f}\equiv J_{0}$ at the initial time $t=0$. When
considering longitudinal spin correlations, i.e. Eq. (72) with
$\alpha=\beta=x$ one can employ the formula

$\displaystyle
G_{xx}\left(r,t\right)=g(r)-\int_{\mathcal{B}}\frac{d^{d}{\mathbf{k}}}{\left(2\pi\right)^{d}}\mathcal{F}\left({k}\right)\frac{e^{i\left(k\cdot
r+2\omega_{{\mathbf{k}}}t\right)}+e^{i\left(k\cdot
r-2\omega_{{\mathbf{k}}}t\right)}}{2},$ (74)

where the integral spans the first Brillouin zone $\mathcal{B}$. In the
following we are going to ignore the time-independent function $g(r)$ and
focus on the time evolution of the correlations, which can be readily obtained

$\displaystyle\mathcal{F}(k)=\frac{2h\left(J_{0}^{\rm i}-J_{0}^{\rm
f}\right)f_{k}\left(\alpha\right)}{\left[h+2J_{0}^{\rm
f}f_{k}\left(\alpha\right)\right]\sqrt{h\left[h+2J_{0}^{\rm
i}f_{k}\left(\alpha\right)\right]}}.$ (75)

(a) Ising Model

(b) XY Model

Figure 6: Spreading of connected spin-spin correlation function. Panel (a)
displays $G_{xx}(r,t)$ for the quantum Ising chain with $\alpha=1.7$ for a
sudden quench for the quench in the paramagnetic phase from $h/J_{0}^{\rm
i}=50$ to $h/J_{0}=1$. The green line is the correlation front which scales
sub-ballistically (white dashed line represents ballistic propagation). The
dashed black line are the analytic scaling $r^{\zeta}$ obtained in Eq. (79).
Panel (b) reports $G_{zz}(r,t)$ for the quantum $XY$ model with $\alpha=3$ in
$d=2$ for a quench starting from the fully polarized state along $x$ and
evolved with the Hamiltonian in Eq. (1) with $\gamma=h=0$ (i.e. a quench from
$\gamma=1$ to $\gamma=0$). The dashed black lines show the scaling of the
maxima which is linear in the axis $r^{\zeta_{XY}}$. Panel (a) is adapted from
Ref. [51] and panel (b) from Ref. [104].

The amplitude of the quench is directly proportional to the difference between
the initial $J_{0}^{\rm i}$ and final $J_{0}^{\rm f}$ couplings. Both these
values are chosen to maintain the system within the paramagnetic phase
$h>h_{\rm cr}$. The time-dependent correlation function $G_{xx}(r,t)$ of the
long-range Ising model obtained by Eq. (74) is displayed in Fig. (6a) for
$\alpha=1.7$. The front of the correlation is highlighted by a green line. Its
scaling is not linear but algebraic as expected for long-range interactions
[121, 122]. Nevertheless, the front propagation does not saturate the
conventional super-ballistic bounds [128], rather it displays a sub-ballistic
scaling, i.e. $t\sim r^{\beta_{\rm front}}$ with $\beta_{\rm front}>1$, which
is represented as a solid green line. Inside the correlation front the scaling
changes and for the Ising model the correlation maxima (light yellow areas)
propagate ballistically with $t\sim r$.

It is interesting to use the stationary phase approximation in order to
evaluate Eq. (74) in the large size and long-time limit. Indeed, for
$t,r\to\infty$ the integral in Eq. (74) is dominated by the configurations
with

$\displaystyle\nabla_{k}\omega_{k}=r/t$ (76)

where the group velocity diverges as $k^{\zeta-1}$ in the $k\to 0$ limit.
Thus, quasi-particles with momentum $k_{\rm sp}=\left(2|c|\zeta
t/r\right)^{1/\left(1-\zeta\right)}$ fulfil Eq. (76) at any given point
$(t,r)$ in space-time. The leading contribution to the correlation front
propagation comes from the low-energy divergence of quasi-particle group
velocity. In order to evaluate the leading contribution to the correlation
function we assume that the amplitude function obeys $\lim_{k\to
0}\mathcal{F}(k)\sim k^{\eta}$, leading to

$G_{xx}(r,t)\propto\frac{t^{\gamma}}{r^{\chi}}\cos\left[A_{\zeta}\left(\frac{t}{r^{\zeta}}\right)^{\frac{1}{1-\zeta}}-2\Delta
t+\frac{\pi}{4}\right],$ (77)

with $\gamma=\frac{\eta+d/2}{1-\zeta}$,
$\chi=\frac{\eta+d(2-\zeta)/2}{1-\zeta}$, and
$A_{z}=2|c|(1-\zeta)(2|c|\zeta)^{\frac{\zeta}{1-\zeta}}$. It follows from Eq.
(77) that the correlation front is obeys the relation $t^{\gamma}\approx
r^{\chi}$ and

$t\propto r^{\beta_{\rm front}},\qquad\beta_{\rm front}=\chi/\gamma.$ (78)

Interestingly, the propagation of the wave-front deos not depend only on the
universal scaling exponent $\zeta$ but also on the specific correlation
function into consideration, since the exponent $\eta$, which describes the
low energy scaling of $\mathcal{F}\left({k}\right)$ in Eq. (77), enters in the
determination of the ration $\chi/\gamma$. As the local limit is approached,
the quasi-particle velocity ceases to diverge $\zeta\to 1$ and linear
spreading of the wave-front is recovered so that Eq. (77) reproduces the Lieb-
Robinson expectation [12]. The relation $\chi=\gamma+d/2$ yields $\beta_{\rm
front}>1$ and imposes sub-ballistic wave-front propagation.

The quench protocol into consideration stays within the paramagnetic phase,
which is characterized by a gapped dispersion relation, see Eq. (73). This, in
turns, leads to $\lim_{k\to 0}\mathcal{F}(k)\sim O(1)$ and the scaling
exponent of the correlation function vanishes, i.e. $\eta=0$. Thus, only the
exponent $\zeta$ determines the front propagation scaling $\beta_{\rm
front}=2+d-\alpha$. The theoretical prediction for $\alpha=1.7$ and $d=1$
produces $\beta_{\rm front}=1.3$, which is in perfect agreement with the one
observed in the numerical computation displayed in Fig. 6a. The formula
$\beta_{\rm front}=2+d-\alpha$ also matches the exact result obtained in Ref.
[60] for $d=1$ and $\alpha=3/2$, also confirmed by t-VMC calculations.

Within the causal region delimited by the wave-front, the local maxima are
determined by the maxima of the cosine function in Eq. (77). Thus, the
correlation maxima occur at the time $t_{\textrm{max}}$, whose value does not
depend on the shape of $\mathcal{F}(k)$, but only on the value of the $\zeta$
exponent, yielding

$t_{\textrm{max}}\propto r^{\zeta},$ (79)

at least for a gapless dispersion relation. According to this analysis, the
maxima of the correlations, located at the time $t_{\max}$ spread super-
ballistically for weak long-range interactions. This has to be contrasted with
the sub-ballistic scaling obtained in Eq. (78) for the front propagation.

The result in Eq. (79) is consistent with the experimental observation on
trapped ions [3] as well as with the truncated Wigner approximation analysis
[126, 129]. However, for the long-range Ising model the dynamical protocol
under consideration remains within the paramagnetic phase, leading to a finite
gap $\Delta\neq 0$. Therefore, the argument of the cosine function in Eq. (77)
is insensitive to the non-analytic scaling of the dispersion relation in the
low-energy limit, becoming constant in the large $t$ and $r$ limit with
$t/r\sim const$. Thus, Eq. (79) has to be substituted with $t_{\rm max}\propto
r$ for a gapped dispersion relations. It follows that the local maxima are
always ballistic, $\beta_{\rm max}=1$ for quenches within gapped phases, see
Fig. 6a. The ballistic motion of local maxima has also been observed with a
trapped ion quantum simulator [14].

Based on the above discussion††margin:  Quenches in the gapless phase , it can
be deduced that the correlations spreading reflect the low-energy properties
of the long-range model. It is evident that the scaling of correlations is
universal in long-range systems, in the sense that it reflects the low-energy
properties of the model. Then, a very different picture is obtained by
studying a quantum quench within a gapless phase. In order to investigate this
dependence we prepare the system in the state fully polarized along the
direction $x$ and evolve it with the Hamiltonian in Eq. (1) with $\gamma=h=0$,
i.e. the long-range XY Hamiltonian, which has $U(1)$ symmetry. Following Ref.
[104] we consider the case $d=2$. The dispersion relation can be obtained by
Eq. (62) by setting $\gamma=h=0$

$\displaystyle\omega_{k}=J_{0}\sqrt{1-f_{{k}}(\alpha)}.$ (80)

As expected, changing the symmetry of the final Hamiltonian modifies the low-
energy dispersion relation, which now scales as $\omega^{\rm low}_{k}\propto
k^{\zeta_{\rm XY}}$ with $\zeta_{\rm XY}=(\alpha-d)/2$, leading to a diverging
quasi-particle group velocity in the $k\to 0$ limit for $d<\alpha<d+2$. A
straightforward computation for the long-range XY model [104, 51] produces
$\eta=\zeta$ leading to $\beta_{\rm front}=1+d(2+d-\alpha)/(2\alpha)$. On the
other hand, Eq. (79) remains unchanged and it yields $\beta_{\rm
max}=\zeta_{\rm XY}$, as it is visible in Fig. 6b and verified by DMRG
calculations in Ref. [104].

Summary: In the weak long-range regime, the correlations front spreads non-
linearly, with exponents that depend on the details of the underlying low-
energy dispersion.

####  3.1.2 Strong long-range regime ($0<\alpha<d$)

We now consider the Hamiltonian (1) in the strong long-range regime
$0<\alpha<d$. Following Ref. [50], in this Subsection the interactions are
_not_ Kac-normalized, i.e. we set $J\equiv J_{0}$ in Eq. (3). This leads us to
discuss the effect of a divergent quasi-particle energy for $k\to 0$ onto the
correlation spreading. [Note that this is different from the rest of the
Report where we focus on the discrete spectrum in Eq. (34) at low $k$!] Within
this framework, we approximate the low-energy dispersion relation with the
expression ††margin:  Divergent dispersion

$\omega_{k}\approx\frac{e_{0}}{k^{\gamma}},$ (81)

where $e_{0}=\sqrt{2hJ_{0}}$ and $\gamma=\frac{d-\alpha}{2}$. Including the
modified dispersion relation into the Eq. (74) one gets

$G_{xx}(r,t)\sim\int_{\omega}d\Omega\int_{\pi/L}^{\pi}dk\,k^{d-1+\gamma}e^{\imath
kr\cos\left(\theta\right)}\left[1-\cos\left(2e_{0}tk^{-\gamma}\right)\right],$
(82)

where the factor $k^{\gamma}$ comes from the low energy limit of the amplitude
function in Eq. (75), i.e. $\mathcal{F}(k)\sim k^{\gamma}$. Due to the
divergent nature of the quasi-particle spectrum, one can introduce a low-
energy cutoff $\sim 1/L$ in the momentum integral in Eq. (82). After expanding
the exponential term Eq. (82) in powers of the distance $r$, the integration
is performed term by term [50]. Then, after discarding finite term in the
system size one finds

$G_{xx}(r,t)\sim\lim_{L\rightarrow\infty}\frac{\sin\left(L^{\gamma}\tau\right)}{\tau}\frac{\int
d\Omega e^{\imath\frac{R}{L}\cos\left(\theta\right)}}{L^{2\gamma+D}},$ (83)

where is the dimensionless time variable $\tau=2e_{0}t$.

Due to the algebraic divergence of the quasi-particles spectrum the time scale
for signal spreading in the system vanishes in the thermodynamic limit.
Accordingly, the vanishing of the signal spreading time displays the same
scaling exponent $\gamma$ as the divergence of the quasi-particle energy,
which for the Ising model reads $\gamma=\frac{d-\alpha}{2}$. This analytic
derivation has been also corroborated by numerical analysis of the spin-wave
dynamics in Ref. [50]. It is worth noting that the same scaling has been
derived within a generalized Lieb-Robinson bound in long-range fermionic
systems [41].

Summary: In the strong long range regime without Kac normalization, the
divergent quasi-particle energy leads to instantaneous correlation spreading
in the thermodynamic limit.

####  3.1.3 Other directions

The ††margin:  Numerical results on the Ising model present Subsection has
been devoted to the study of correlation spreading within linear spin-wave
theory [60, 50, 51, 104]. This approximation proved capable to capture the
salient features of several numerical simulations [59, 124, 130, 127, 126]. In
particular, in the case of the quantum Ising chain [Eq. (1) with $d=1$],
numerical matrix-product state calculations have shown that the emergence of a
short-range-like light-cone behavior for $\alpha>2$ [127] as confirmed by the
study in Sec. 3.1.1. On the other hand, for $1<\alpha<2$, the model displays
clear light cone but with an infinite propagation speed of almost all
excitations. This is linked to the divergence of the maximum group velocity,
which leads to a scenario of multispeed prethermalization [104], see again
Sec. 3.1.1. For $\alpha<1$, all studies report a clear nonlocal regime, with
instantaneous transmission of correlations between distant sites, in agreement
with the study reported in Sec. 3.1.2.

Despite ††margin:  Extensions the successes of linear spin-wave theory, it
would be interesting to reconsider these results using the time-dependent
framework that we will present in Sec. 4.2.1. One may compute $G(r,t)$ by
expressing spin operators in a time-dependent frame and expanding them using
Holstein-Primakoff bosons. This very analysis has been performed has been
performed in a related model in Ref. [57] in connection with dynamical phase
transitions; see also the analogous calculation for scrambling dynamics in Eq.
(143).

Finally, ††margin:  Rigorous results a very successful research direction
based on rigorous mathematical investigations was pursued to generalize the
Lieb-Robinson bound for power-law decaying interactions [123]. In a seminal
work of 2006, Hastings and Koma [121] showed that it takes a time
$t\gtrsim\log r$ to propagate information at distance $r$ for all $\alpha>d$.
However this bound is far from tight, since it does not recover the linear
light-cone in Eq.(71) for large $\alpha$. Several efforts in the past years
have led to a greatly improved picture [41, 131, 132, 133, 134, 135, 136, 122,
137, 36, 138, 38, 139]. Firstly, it was proved the existence of the linear
light-cone $t\gtrsim r$ for $\alpha>2d+1$ [40, 140]. Secondly, it was shown
that this result becomes $t\gtrsim r^{\min(\alpha-2d,1)}$ for $\alpha>2d$
[139, 120]. On the other hand, in the strong long-range regime $0<\alpha<d$,
correlations between distant degrees of freedom can propagate instantaneously,
since the bounds on the light-cone time-scale can vanish with the system size
[41, 39, 137]. So far, the best estimate for interacting systems is $t\gtrsim
N^{\alpha/d-1}\log N$, that can be made tighter $t\gtrsim N^{\alpha/d-1/2}$
for free models with $\alpha<d/2$ [39]. Notably, the violations of Lieb-
Robinson bound have been experimentally probed on trapped ions quantum
simulators for $0.6\lesssim\alpha\lesssim 1.2$ in Refs. [14, 3].

Summary: In addition to low-energy approximations, the spreading of
correlations has been tackled with various approaches ranging from numerical
simulations to mathematically rigorous bounds. The current established
scenario is rather complete and satisfactory.

###  3.2 Metastability

In the following, we are going to show that metastability in quantum strong
long-range systems may be traced back to their discrete quasi-particle
spectrum, which hinders the applicability of the kinematical chaos hypothesis
[141].

####  3.2.1 State of the art

Diverging equilibration times in the thermodynamic limit are a well-known
feature of long-range interacting systems. This phenomenon is widespread
within the classical long-range physics world [142, 143], but multiple
theoretical observation occurred also in the quantum realm [42, 144]. Long-
lived pre-thermalization is also expected to occur in cold atom clouds
confined into optical resonators [145], where semi-classical analysis of the
Fokker-Planck equation directly connects to the Hamiltonian mean-field model
[146, 147], the workhorse of classical long-range physics [148]. Recent
studies have directly linked the absence of equilibration in strong long-range
quantum systems to the discreteness of the quasi-particle spectrum, see Eq.
(34a). This results in a violation of Boltzmann’s H-theorem and leads to the
emergence of finite Poincaré recurrence times in the thermodynamic limit [52].
This section discusses the appearance of diverging equilibration times for
quantum long-range systems in the thermodynamic limit. This is consistent with
the properties discussed in Section 2.4, which are common to both large long-
range systems and finite local ones. Examples include the inability to
completely disregard boundary effects over bulk phenomena [149, 150], the
existence of concave entropy regions [151], and the presence of a macroscopic
energy gap between the ground state and the first excited state [152, 153]. It
is worth noting that our description mostly pertains isolated quantum systems,
while multiple theoretical and experimental observations in cavity systems
evidenced a substantial role of dissipation [154, 155].

The crucial aspect is that the excitation spectrum of non-interacting systems
does not become continuous in the thermodynamic limit. The eigenvalues of a
long-range coupling matrix have been shown to remain discrete even in the
infinite components limit, forming a pure point spectrum [156] similar to that
observed in strongly disordered systems [157, 158, 159, 160]. A discussion on
the spectral discreteness of long-range couplings in the thermodynamic limit
has been presented in Ref. [52] for a few quadratic models and used to explain
the observation of diverging equilibration times in a long-range Ising model,
quenched across its quantum critical point [42]. We refer the readers to Sec.
2.4 and C.

####  3.2.2 Quasi-stationary states and spectral properties

The ††margin:  Quasi-stationary states first evidence of quasi-stationary
states (QSS) in quantum systems was described in the prototypical example of
the long-range quantum Ising chain [see Eq. (1)]. QSS were shown to appear for
quenches starting well inside the paramagnetic phase in the $h\to\infty$ limit
and ending deep in the ferromagnetic phase at $h=0$. Here, the system is
prepared in the transversally polarized ground state and evolved according to
the classical ferromagnetic Hamiltonian in Eq. (1) in the absence of the
transverse field $h=0$. As a result, the expectation of the global operator
$m_{z}=\langle\sum_{i}\sigma^{z}_{i}\rangle/N$ evolves from the initial value
$\lim_{t\to 0}m_{z}=1$ to the equilibrium expectation
$\lim_{t\to\infty}m_{z}=0$, if the system actually equilibrates. These
observations have been extended to any choice of the initial and final
magnetic fields $h_{i}$, $h_{f}$ using the Jordan-Wigner representation of the
Ising model.

The appearance of the QSS has been frequently linked to the scaling of
equilibration times of critical observables, such as the magnetization [161,
143, 47]. However, persistent time fluctuations have also been found in
generic thermodynamic observables of classical systems, such as the evolution
of internal energy in systems of particles with attractive power-law pair
interactions [95]. Similar phenomena ††margin:  Low-energy quenches can be
observed in our system, by considering just the leading order low-energy
theory. In order to simplify the study we restrict our analysis to the
paramagnetic quantum Ising chain, whose quasi-particle dispersion is

$\displaystyle\omega_{k}=\sqrt{h(h-2J_{0}f_{k}(\alpha))}\,,$ (84)

cf. Eq. (47). It is worth noting that the present spin-wave approximation
corresponds to the time-dependent Hartree-Fock approximation of the Ising and
$O(N)$ rotor models. Accordingly, several phenomena occurring in the out-of-
equilibrium low-energy dynamics of the Ising Hamitlonian can be also observed
in the large-$N$ limit of $O(N)$ models [162, 163], including
prethermalization [164, 165], defect formation [166], dynamical phase
transitions [109]. In particular, the dynamics induced by a sudden quench
leads to universal relaxation properties [163, 167, 109].

In this regime, equilibration does not occur in the non-additive regime due to
the discrete quasi-particle spectrum $f_{k}(\alpha)$. In order to demonstrate
this fact let us consider a sudden magnetic field quench $h^{\rm i}\to h^{\rm
f}$ in the Hamiltonian (46). The quench occurs within the normal phase
$h>h_{\rm cr}$ so that no magnetization occurs. Nevertheless, a finite spin-
wave density will arise due to the sudden quench and will contribute to the
evolution of any internal observable of the system. In order to make a direct
parallel with the classical case described in Ref. [168] we consider the
evolution of the spin-wave Kinetic energy

$\displaystyle
K(t)=\sum_{k}\langle\hat{p}_{k}^{2}\rangle/2N=-\frac{\omega_{k}}{4}\left\langle\left(\beta^{\dagger}_{k}-\beta_{k}\right)^{2}\right\rangle\
,$ (85)

where $\beta_{k}$ and $\beta^{\dagger}_{k}$ diagonalize the quadratic
Hamiltonian in Eq. (42). The ††margin:  Solution of the quadratic dynamics
calculation is rather straightforward, since the system is assumed to lie in
the ground-state before the sudden quench. After the quench, each spin-wave
occupies a squeezed state, so that the system lies in the quantum state
$\Pi_{k}\hat{S}_{k}(\zeta)|0\rangle$, where $|0\rangle$ is the vacuum and the
squeezing operator $\hat{S}_{k}(\zeta)|0\rangle$ reads

$\displaystyle\hat{S}_{k}(\zeta)=\exp\frac{\left(\zeta^{*}(\hat{\beta}_{k})^{2}-\zeta(\hat{\beta}^{\dagger}_{k})^{2}\right)}{2}.$
(86)

The squeezing parameter $r$ is defined by rewriting $\zeta$ in polar
coordinates $\zeta=re^{i\phi}$. Then, it is rather straightforward to rewrite
the squeezing parameter in terms of the effective oscillator length $\xi(t)$

$\displaystyle\tanh{r_{k}}=\sqrt{\frac{\left(\frac{1}{2\xi_{k}(t)^{2}}-\omega_{k}\right)^{2}+\frac{\dot{\xi_{k}}(t)^{2}}{\xi_{k}(t)^{2}}}{\left(\frac{1}{2\xi_{k}(t)^{2}}+\omega_{k}\right)^{2}+\frac{\dot{\xi}_{k}(t)^{2}}{\xi_{k}(t)^{2}}}}.$
(87)

To ††margin:  Persistent oscillations obtain the spin-wave dynamics is then
sufficient to solve the Ermakov equation, which describes the evolution of the
effective length [169]

$\displaystyle\ddot{\xi}_{k}(t)+\omega_{k}^{2}(t)\xi_{k}(t)^{2}=\frac{1}{4\xi_{k}(t)^{3}}.$
(88)

(a) Kinetic Energy

(b) Dynamical fluctuations

Figure 7: Equilibration of long-range spherical model. Panel (a) displays the
dynamical evolution of the kinetic energy following the sudden quench
$h_{i}\to h_{f}$. After a steady decay during the initial dynamics $t\lesssim
10^{2}$, dynamical oscillations set to a finite value that remain steady in
the long-time regime. The amplitude of dynamical fluctuations for the kinetic
energy after a time $T$ is quantified by the quantity $Q_{K}(T)$ see Eq. (89)
displayed in panel (b).

The solution of the sudden quench dynamics is readily obtained by introducing
$\omega_{k}(t)=\theta(-t)\,\omega_{k,\rm i}+\theta(t)\,\omega_{k,\rm f}$ in
Eq. (88). The resulting dynamical evolution for the spin-wave kinetic energy
is displayed in Fig. 7 for $\alpha\in[0.15,0.35,0.65,0.95]$. In analogy with
the classical case the observable $K(t)$ displays persistent dynamical
oscillations, which do not wash out in the thermodynamic limit. The smaller
the $\alpha$ the wider the amplitude of those fluctuations.

A ††margin:  Origin due to the discrete spectrum simple explanation of this
phenomenon is found into the fully connected limit ($\alpha\to 0$), where the
function $f_{k}(\alpha)$ separates between two distinct energy levels in the
thermodynamic limit: a non-degenerate ground-state with energy $-J$ and a
$N-1$ degenerate excited states with zero energy.

In presence of any given set of boundary conditions, the degeneracy is split
and the system behaves at finite size as a set of harmonic oscillators with
discrete energies. As the size increases, the spectrum accumulates at high
energy where the eigenvalues $f_{k}(\alpha)$ of the coupling matrix become all
identical, making the system equivalent to a single quenched harmonic
oscillator. To characterize equilibration, we introduce the characteristic
function of any observable $A$, i.e. $\chi_{A}(t)=A(t)-\bar{A}$. This quantity
captures the dynamical fluctuations around the average value of the
observable. Equilibration of the observable $A(t)$ in closed quantum systems
occurs when the long-time Cesaro’s average of the squared fluctuation vanishes
[170, 171, 172]:

$\displaystyle\lim_{T\to\infty}Q_{A}(T)\equiv\lim_{T\to\infty}\langle|\chi_{A}(t)|^{2}\rangle_{T}=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}|\chi_{A}(t)|^{2}dt=0\
,$ (89)

while metastability shall be associated with a finite value
$\lim_{T\to\infty}Q_{A}(T)\neq 0$. The equilibration of the kinetic energy $K$
in Eq.(85) – or other physical observables – follows a similar argument as the
fidelity of a quantum system in the context of the spectrum of long-range
systems [52].
In _weak-long range interacting systems_ with translational invariance, the
spectrum becomes absolutely continuous in the thermodynamic limit. This
implies that $\lim_{t\to\infty}\chi_{K}(t)\to 0$ outside of the Cesaro’s
average due to the Riemann-Lebesgue lemma. For quantum systems with initial
states having no overlap with pure point portions of the spectrum,
equilibration as defined in equation (89) is ensured by Wiener’s theorem
[156]. Considering the thermodynamic limit of a _strong long-range
interacting_ system, where the system size $N$ increases, the eigenmodes
$f_{k}(\alpha)$ of the Hamiltonian tend to accumulate at high energy, near
$\omega_{k}|_{h_{f}}\sim 2h_{f}$. In fact, in the case of _flat interactions_
($\alpha=0$), the spectrum consists of a single infinite degenerate
eigenstate, resulting in dynamics that precisely correspond to a single
harmonic oscillator.

Summary: In the strong long-range regime, the discrete low-energy spectrum
yields quasi-stationary states. These manifest in persistent oscillations of
time-dependent observables, such as the kinetic energy, with an amplitude that
increases as the interaction becomes longer ranged $\alpha\to 0$.

####  3.2.3 Equilibration in presence of disorder

We would like to comment on the relation between these results and the
metastability that arises in disordered systems. From the perspective of
spectral discreteness, the metastability observed for quantum systems in the
strong long-range regime is fundamentally different. Indeed, ††margin:  Random
fully connected model when flat interactions ($\alpha=0$) are perturbed by
Gaussian distributed weak couplings $u_{ij}$, namely $J^{\rm
dis}_{\boldsymbol{r},\boldsymbol{r}^{\prime}}=J_{\boldsymbol{r},\boldsymbol{r}^{\prime}}+u_{ij}$
with

$\displaystyle P(u_{ij})\propto\exp\left(-N\,u_{ij}^{2}/2w^{2}\right)$ (90)

whose width $2w$ represents the disorder strength. The disordered couplings
lift the infinite ($\sim N-1$) degeneracy of the excited state at zero energy
and the spectrum becomes continuous apart from the single non-degenerate
ground-state at energy $J_{0}$, where $J_{0}>w$ is the strength of the flat
homogeneous interactions [173, 174, 175]. The density of states of the
continuous spectrum follows the celebrated Wigner semicircle law [176]. In
analogy with the non disordered case, we initialize the dynamics at
equilibrium for $h=2.2h_{\rm cr}$ at $t\leq 0$; then, at $t>0$ the magnetic
field suddenly switches at $h_{f}=1.1h_{\rm cr}$. The continuum nature of the
spectrum leads the spin-wave kinetic energy $K(t)$ to exponentially
equilibrate, see Fig. 8a.

(a) Dynamical Fluctuations

(b) Equilibration Time

Figure 8: Equilibration of the disordered long-range Ising model within spin-
wave approximation. Panel (a): Dynamical fluctuations decay as a function of
the disorder strength for the kinetic energy $K$, see Eq. (85). As the
disorder strength is decreased the decay rate also decreases until it vanishes
in the zero disorder limit (upper blue curve), where dynamical fluctuations
persist at all times and equilibration never occurs. Panel (b): the
equilibration time of the kinetic energy observable obtained by fitting the
curves in panel (a) via the exponential form in Eq. (91). The divergence in
the clean case ($2w\to 0$) is evident. Figures reproduced from Ref. [52].

Indeed, the amplitude of dynamical fluctuations in disordered systems decays
exponentially, allowing for the introduction of the equilibration time
$\tau_{\rm eq}$:

$\displaystyle Q_{K}(T)\sim e^{-T/\tau_{\rm eq}}\ ,$ (91)

Numerical analysis shows that the equilibration time $\tau_{\rm eq}$
monotonically decreases with increasing disorder strength $J$, as expected
(see Fig. 8b). The exponential decay of dynamical fluctuations and the
definition of the equilibration time provide insights into the equilibration
dynamics of disordered systems. Interestingly, the results obtained in Fig. 8
for quantum systems show remarkable similarities to those obtained for
classical spherical models with disordered couplings, as shown in Chapter 4 of
Ref. [177]. There, the Langevin dynamics of the disordered classical spherical
model it is shown not to exhibit metastability, as long as the initial state
is not magnetized.

The ††margin:  Effective dimensionality analysis presented here focuses on
characterizing the long-time equilibration dynamics of many-body quantum
systems, where the thermodynamic limit is generally taken before the long-time
limit. However, similar conclusions can be obtained by considering the long-
time limit of dynamical fluctuations in finite systems, which yields [170,
171, 172]:

$\displaystyle\lim_{T\to\infty}Q_{K}(T)\propto\frac{1}{d_{\rm eff}}$ (92)

Here, $d_{\rm eff}$ roughly represents the number of modes participating in
the dynamics. In the case of finite systems, where the entire spectrum is
discrete and only a finite number of modes exist, $d_{\rm eff}$ captures this
finite nature. As the thermodynamic limit is approached and the spectrum
becomes continuous, eigenvalues become dense in arbitrarily small energy
ranges. In the continuous limit, $d_{\rm eff}\to\infty$ for dynamics involving
initial states in the continuous spectrum. However, for long-range systems
with $\alpha<d$, where $d$ is the dimension of the system, a continuous theory
cannot be defined. This is because the only dense region in the spectrum
occurs around the energy maximum, where infinitely many degenerate eigenvalues
emerge. This violates the assumption of non-degenerate energy gaps underlying
equation (92) [178].

Summary: Metastability for strong long-range systems is fundamentally
different from the disordered one in that only a finite number of effective
degrees of freedom participate in the dynamics.

###  3.3 Kibble-Zurek mechanism

Within ††margin:  Landau-Zener problem the realm of quantum systems, the
Landau-Zener problem provides the earliest, and possibly simplest, example of
defect formation during a quasi-static drive [81, 179, 180]. The problem
describes a two-level system slowly driven across an avoided level crossing.
Although initialized in the ground state at the initial time $t_{i}<0$, the
system gradually populates the excited state, whose energy separation slowly
decreases during the dynamics. The energy gap starts from a minimum at $t=0$
and increases back until the time reaches the endpoint of the dynamics
$t_{f}>0$. The dynamical evolution is controlled by the rate parameter
$\delta$, i.e. $H(t)\equiv H(\delta\cdot t)$, so that the quasi-static limit
is reached as $\delta\to 0$. A straightforward criterion to establish whether
a quasi-static transformation remains adiabatic is to ensure that the rate of
change of the instantaneous minimal gap $\Delta=|E_{0}-E_{1}|$ remains smaller
than the square gap itself

$\displaystyle\dot{\Delta}(t)\ll\Delta^{2}(t).$ (93)

Latter criterion only involves equilibrium quantities since $E_{\ell}(t)$
represents the spectrum of the instantaneous Hamiltonian at the time $t$.

While Eq. (93) has been obtained by heuristic arguments a more rigorous
derivation of an adiabatic criterion and a discussion of how it compares with
Eq. (93) can be found in Ref. [181]. Eq. (93) has been introduced for the
Landau-Zener problem, but the argument in Ref. [181] applies to generic
quantum systems. In general, as long as Eq. (93) is satisfied, the excited
state population of a quasi-statically driven quantum system decreases with
the drive rate and the hypothesis for the quantum adiabatic theorem are
satisfied [182]. However, as the drive approaches a quantum critical point the
correlation length of a quantum system diverges and the instantaneous gap
vanishes $\Delta\to 0$. As a result, the dynamical scaling of the observables
close to the quantum phase transition is reminiscent of the thermodynamic
scaling at equilibrium. Yet, in order for such scaling to be displayed, the
drive has to be slow enough that the dynamical evolution actually occurs in
the vicinity of the equilibrium critical point.

Let us††margin:  Dynamical protocol focus on the concrete case of the
Hamiltonian in Eq. (1) whose internal control parameter is defined as
$\lambda=h(t)-h_{\rm cr}$, such that the ferromagnetic quantum critical point
occurs at $\lambda_{c}=0$. For a moment, let us imagine that the system is
finite, so that the spectrum remains gapped also at criticality. Then, the
hypothesis of the quantum adiabatic theorem [182] remains fulfilled and any
slow enough drive of internal parameters $\lambda(t)\sim\delta t$ only
generates adiabatic corrections $\sim\delta^{2}$ to the observables
expectations with respect to the equilibrium value, as it can be deduced by
simple thermodynamic arguments [183]. However, in the thermodynamic limit,
crossing the equilibrium critical point breaks down the conventional adiabatic
picture and the residual energy (heat) generated by the drive displays non-
analytic behaviour $E_{\rm res}\approx\delta^{\theta}$ with $\theta<2$ [184].
Our task within the present section is the determination of the universal
scaling index $\theta$ for quantum long-range systems.

####  3.3.1 State of the art

The ††margin:  Kibble-Zurek Kibble-Zurek mechanism allows to _relate the value
of the non-analytic exponent $\theta$ with the equilibrium critical exponents
of the model_. This ingenious theory relies on the adiabatic-impulse
approximation, which assumes that the dynamical evolution of a system starting
in its ground-state at $t=-\infty$ adiabatically follows the drive until the
freezing time $-\hat{t}$. Beyond the freezing time, the equilibration rate of
the system becomes too small with respect to the drive velocity, violating the
adiabatic condition in Eq. (93). Then, the freezing time satisfies the
condition $\dot{\Delta}(\hat{t})/\Delta(\hat{t})^{2}=1$. Due to the critical
scaling of the instantaneous gap $\Delta\propto\lambda^{z\nu}$, latter
relation leads to the freezing time inheriting the equilibrium critical
scaling $\hat{t}\propto\delta^{1/(1+z\nu)}$. For $t>-\hat{t}$, the state
dynamics is assumed to remain frozen for the entire interval
$t\in[-\hat{t},\hat{t}]$ until the unfreezing time $\hat{t}$, where
adiabaticity is restored (for simplicity we have assumed a symmetric gap).

Once the system has unfrozen the state evolution will resume on the opposite
side of the transition, where the Hamiltonian ground-state is supposed to
break the Hamiltonian symmetry. Then, the dynamics will induce a transition