between the symmetric and a symmetry-broken state. However, this transition
will occur at finite correlation length $\hat{\xi}$, since the process can
only start at the time $t\geq\hat{t}$ well within the symmetric phase of the
model. The dynamics has thus modified the character of the continuous phase
transition, making it rather similar to a first-order one, and the system will
likely form topological defects, whose size would be roughly proportional to
the (finite) correlation volume $\hat{\xi}^{d}$. Therefore, the total defect
density scales according to $N_{\rm
exc}\propto\hat{\xi}^{-d}\propto\delta^{d\nu/(1+z\nu)}$ [185].

Several ††margin:  Short-range scaling verifications of the Kibble-Zurek
scaling exist in local systems, both via numerical simulations, exact
theoretical studies and experiments [186]. In particular, first studies of
defect formation in quantum systems have been pursued on the Hamiltonian in
Eq. (1) in the $\alpha\to\infty$ limit, i.e. the nearest neighbour Ising
model, where finite size scaling arguments led to the prediction

$\displaystyle N_{\rm exc}^{\rm fss}\approx\delta^{\frac{1}{2z}}$ (94)

where the superscript fss stands for finite size scaling. Eq. (94) produces
$N_{\rm exc}\approx\sqrt{\delta}$ in agreement with the Kibble-Zurek
prediction $N_{\rm exc}\propto\delta^{d\nu/(1+z\nu)}$ since $z=\nu=1$ in this
case [187]. Soon after this seminal investigation, an exact solution to the
universal slow dynamics of the Ising model has been provided by mapping it to
a infinite sum of Landau-Zener problems, each representing the dynamics of a
single fermionic quasi-particle [188]. This exact solution provides a
different scaling theory for the defect density of the Ising model, which is
given by

$\displaystyle\int N_{\rm exc}(k)dk\approx\delta^{\frac{1}{2z_{\Delta}}}$ (95)

where we have defined $z_{\Delta}$ from the scaling of the dynamical gap. The
result in Eq. (95) has been also employed to prove validity of the Kibble-
Zurek argument in Kitaev chains with long-range pairing terms [189] where
$z_{\Delta}=z$, as well as in the perfect local case $\alpha=\infty$ [188].

Apart ††margin:  Fermi systems for the aforementioned results, which
explicitly refer to quadratic Fermi systems, the application of adiabatic
perturbation theory to slow quenches close to quantum critical points predicts
the scaling of the defect density to be in agreement with the classical
Kibble-Zurek prediction

$\displaystyle\int N^{\rm KZ}_{\rm
exc}(k)\,dk\approx\delta^{\frac{d\nu}{1+z\nu}},$ (96)

which also leads to the scaling exponent $\theta=z\nu/(1+z\nu)$ for the
residual energy [190]. Such prediction comes from the assumption that the
scaling form of the critical propagator reproduces the equilibrium critical
exponents. Since for $1d$ Fermi systems, one has $z\nu=1$, the perturbative
argument yields $d\nu/(z\nu+1)=1/2z$ in agreement with the finite size scaling
argument in Eq. (94). Interestingly, long-range anisotropic interactions with
different ranges depending on the type of coupling are known to violate the
perturbative assumption [191] even in the finite range case [192, 193, 185].

In summary, the applicability of the Kibble-Zurek result in quantum systems is
supported by two main arguments. The finite size scaling argument reported in
Eq. (94) and the perturbative argument, which reproduces the traditional
Kibble-Zurek scaling in Eq. (96). Both the arguments coincide for local
quantum many-body systems where the femionic quasi-particle description
applies. This is the case of the Ising Hamiltonian in Eq. (1) with $\alpha\gg
d+3$ as confirmed by the exact solution obtained at $\alpha=\infty$.

First ††margin:  Fully-connected quantum Ising model indications that the
scaling of the defect density in the $\alpha=0$ Ising model did not follow the
Kibble-Zurek prediction appeared in Ref. [194]. However, later investigations
were found to obey the Kibble-Zurek mechanism, at least for slow ramps
terminating at the critical point, i.e. $t\in[-1/\delta,0]$ [106]. This
apparent inconsistencies triggered more intensive numerical studies, which
unveiled a complicated landscape where the adiabatic crossing of the
equilibrium quantum critical point does not display any scaling with the ramp
rate $\delta$, but rather featured a novel form of dynamical universality as a
function of the scaled variable $\Lambda=N\,\delta$ [105].

(a) Slow quench to the quantum critical point

(b) Slow quench across the quantum critical point

Figure 9: Defect formation in the $\alpha=0$ Ising model during a slow quench:
Panel (a): Residual energy as a function of the drive rate $\delta$ for
different values of the final gap for a slow dynamics terminating exactly at
the critical point. Panel (b): Residual energy after a full ramp across the
quantum critical point for two different system sizes and three different
values of the universal scaling variable $\Lambda=N\delta=15,3.75,0.94$ from
top to bottom.

Summary: In locally-interacting systems the Kibble-Zurek mechanism is
supported by finite-size scaling and perturbative arguments, which however
break down in presence of long-range interactions.

####  3.3.2 Quasi-static dynamics for $\alpha=0$

The mosaic can be easily recomposed by the study of the slow drive dynamics
within the linear spin-wave theory in Eq. (46). This strategy coincides with
the one employed in Sec. 3.2, but with two important differences: first, since
we are considering a quasi-static drive, the dynamical evolution remain close
to the instantaneous equilibrium state and the only relevant source of
deviations from adiabaticity originates from the lowest energy mode.
Therefore, we can safely limit ourselves to the case $\alpha=0$, where just a
single spin-wave exist. Secondly, and more importantly, we are going to
consider a time dependent magnetic field of the form $h(t)=h_{\rm
cr}+\delta\,t$ for $t\in[-h_{\rm cr}/\delta,h_{\rm cr}/\delta]$, so that the
dynamics initiates in the ferromagnetic state before crossing the critical
point. A full treatment of the ferromagnetic state dynamics shall also include
the motion of the classical magnetization, whose coupling with the quantum
modes is suppressed by a factor $1/N$ in the thermodynamic limit. In the
following, we are going to discard the contribution of the classical mode to
the dynamics, since a classical variable in a bounded (singular) potential
generates a correction which scales at most as $\sim\delta^{2}$ and is,
therefore, negligible with respect to the contribution of the quantum mode
[81].

Within the aforementioned assumptions, ††margin:  Spin-wave evolution the
quasi-static dynamics of the $\alpha=0$ Ising model reduces to the evolution
of a single spin-wave. As for the sudden quench case, see Sec. 3.2, the
dynamics initialized in the ground state remains in a squeezed state at all
times [195, 196, 197]. Then, the dynamics generates any two particle states as
follows from Eq. (86). We focus on a cyclic transformation where the system is
initially in the ground state of the equilibrium Hamiltonian, thus, it is
convenient to rewrite the single spin-wave state as

$\displaystyle\psi_{0}(x,t)=\left(\frac{1}{2\pi\xi^{2}(t)}\right)^{\frac{1}{4}}e^{-W(t)\frac{x^{2}}{2}}e^{-i\frac{\varphi(t)}{2}}.$
(97)

with the effective time dependent frequency
$W(t)=-i\frac{\dot{\xi}(t)}{\xi(t)}+\frac{1}{2\xi^{2}(t)}$, and the
ininfluential phase
$\varphi(t)=\int^{t}\frac{dt^{\prime}}{2\xi^{2}(t^{\prime})}$. Thus, even in
the linear ramp case the entire dynamics is described by the differential Eq.
(88).

In order to determine the excitation density and the ground state fidelity
with respect to the instantaneous equilibrium solution of the problem, we
define the adiabatic basis $\psi_{n}^{\rm ad}(x,t)$, which is obtained taking
the equilibrium spin-wave eigenstates and replacing the constant frequency
with the time-dependent one [169]. Accordingly, one can expand the exact time-
dependent state in terms of the adiabatic basis $\psi(x,t)=\sum
c_{n}(t)\psi^{\rm ad}_{n}(x,t)$, leading to the following results for the
excitation density

$\displaystyle N_{\rm exc}(t)=\langle\hat{n}\rangle=\sum_{n\in
2\mathbb{N}}n|c_{n}|^{2}=\frac{\xi(t)^{2}}{2\omega(t)}\left[\left(\frac{1}{2\xi(t)^{2}}-\omega(t)\right)^{2}+\left(\frac{\dot{\xi}(t)}{\xi(t)}\right)^{2}\right],$
(98)

and the adiabatic ground-state fidelity

$\displaystyle
f(t)=|c_{0}|^{2}=\frac{1}{\xi(t)}\sqrt{\frac{2\omega(t)}{\left(\frac{1}{2\xi(t)^{2}}+\omega(t)\right)^{2}+\left(\frac{\dot{\xi}(t)}{\xi(t)}\right)^{2}}}.$
(99)

Interestingly, one can relate the former expressions to the squeezing
paramenter in Eq. (87) by the simple relation $\tanh(r)=\sqrt{N_{\rm
exc}(t)f(t)^{2}}$.

An ††margin:  Linear ramp analytic solution can be found for a linear ramp
across the quantum critical point with the resonant spin-wave having the
dynamical frequency

$\displaystyle\omega(t)^{2}=4h(t)(h(t)-2J_{0})\approx 8\delta|t|$ (100)

where the last expression on the r.h.s. has been obtained substituting
$h(t)=h_{\rm cr}-\delta t$ and expanding for small $\delta t$. The linear
scaling of $\omega(t)^{2}$ is the consequence of the gap scaling $z\nu=1/2$ of
the equilibrium problem, see Eq. (47). Eq. (101) represents the perfect
crossing of the quantum critical point, since the instantaneous spin-wave
frequency perfectly vanishes at $t=0$. In order to effectively incorporate
finite size effects, we shall introduce a small deviation from perfect
degeneracy and rewrite Eq. (101) as

$\displaystyle\omega(t)^{2}=8\delta|t|+\Delta_{N}^{2},$ (101)

where $\Delta_{N}$ is the minimal gap of the finite size system. Obviously,
$\lim_{N\to\infty}\Delta_{N}\to 0$ and the system attains perfect criticality
in the thermodynamic limit. Moreover, due to universality, the minimal gap
display power-law scaling of the form $\Delta_{N}^{2}\approx N^{-1/\nu_{*}}$
with $\nu_{*}=3/2$ as predicted by finite size scaling theory [198] and
confirmed by exact studies on the fully-connected quantum Ising model and
related flat interacting models [84, 199, 200, 201].

The model in Eq. (101) describes a cyclic transformation of the single
Hamiltonian mode and, in the limit $\delta\to 0$, it can be used to describe a
quasi-static cycle in quantum systems with infinitely degenerate spectrum.
According to the behaviour of the fidelity and excitation density in the
quasi-static limit $\delta\to 0$ the system presents three stages

  1. 1.

Perturbative regime ($N<\infty$).

  2. 2.

Kibble-Zurek regime ($N\to\infty$ and $t\in[-h_{\rm cr}/\delta,0]$).

  3. 3.

Non-adiabatic regime ($N\to\infty$ and $t\in[-h_{\rm cr}/\delta,h_{\rm
cr}/\delta]$).

Regime (1) occurs for a finite minimal gam $\Delta_{N}>0$: there adiabatic
perturbation is applicable and the dynamics remains adiabatic, i.e. $N_{\rm
exc}\propto\delta^{2}$ [202]. Regime (2) is realised for a thermodynamic
system ($\Delta_{N}\to 0$) whose dynamics terminates exactly at the quantum
critical point $t=\Delta_{\infty}=0$, where non-analytic corrections of the
form $\delta^{\theta}$ appear in the residual energy. As we are gonna see in
the following this regime is properly described by the Kibble-Zurek argument.
An actual crossing of the quantum critical point only occurs in regime (3) and
the system enters the non-adiabatic regime, where the residual energy and the
fidelity acquire dynamical correction which do not depend on the drive rate.

The latter result can be easily shown rephrasing Eq. (88) in a rate
independent form via the transformations

$\displaystyle
t=\delta^{-\frac{1}{3}}\tilde{t},\quad\xi=\delta^{-\frac{1}{6}}\tilde{\xi}$
(102)

which reduce Eq. (88) to the $\delta=1$ case. The expressions in Eqs. (98) and
(99) are invariant under the transformations in Eq. (102) in such a way that
the fidelity and excitation density at real times can be obtained by
$\tilde{\xi}_{\Delta}(\tilde{t})=\lim_{\delta\to 1}\xi_{\tilde{\Delta}}(t)$.
The subscript $\Delta$ has been introduced to explicit the dependence on
$\Delta$ of the solution of Eq.(88) In the new variables, the only dependence
of the dynamics on the rate $\delta$ remains in
$\tilde{\Delta}_{N}=\Delta_{N}/\delta^{1/3}$. While the transformations (102)
have been reported for the case of a linear quench, they can be easily
generalized to any non-linear drive $\lambda(t)=(\delta|t|)^{\tau}$, obtaining
results analogous to the one described in the present section [54]. Thus, the
invariance of Eq.(88) with respect to the rescaling in Eq. (102) is enough to
demonstrate that the dynamical evolution of the system only depends on the
combined variable $\Lambda=\delta N=\tilde{\Delta}_{N}^{-3}$ as first
evidenced by the numerical study in Ref. [105].

Summary: For a quasi-static drive terminating at the critical point Kibble-
Zurek scaling is observed, while for dynamical protocols crossing the critical
point the amount of defects is independent of the quench rate.

####  3.3.3 Adiabaticity breaking

However, in order to provide estimates for the defect density and fidelity in
the quasi-static limit one has to solve Eq. (88) exactly. In the following we
are going to drop all the $\sim$ superscripts over the rescaled variables, in
order to ease the notation. The crucial condition of adiabatic dynamics is for
the system to start in the ground-state at the beginning of the dynamics, i.e.
$\displaystyle{\lim_{t\to-\infty}}\psi(t)=\psi_{0}^{\rm ad}(t)$, leading to
the boundary conditions

$\displaystyle\lim_{t\to-\infty}\xi(t)^{2}=\frac{1}{2\omega(t)};\quad\lim_{t\to-\infty}\dot{\xi}(t)^{2}=0.$
(103)

First, ††margin:  Regime (2) it is instructive to consider the case of the
quasi-static dynamics terminating at the critical point, i.e. $t\in[-h_{\rm
cr}/\delta,0]$. As anticipated, the solution only depends on the combination
$\Lambda=N\delta$ and, then, the corrections observed in the dynamics
dramatically depend on what limit is taken first, see D. If one considers the
quasi-static limit ($\delta\to 0$) at finite size, the rescaled gap
$\Delta_{N}$ diverges and adiabatic corrections arise in all dynamical
quantities

$\displaystyle\lim_{\delta\to 0}N_{\rm
exc}(t_{0})=o\left(\delta^{2}\right);\qquad\lim_{\delta\to
0}f(t_{0})=1-o\left(\delta^{2}\right).$ (104)

More interestingly, if the thermodynamic limit is taken first, the rescaled
instantaneous gap vanishes $\lim_{N\to\infty}\Delta_{N}=0$. However, this is
not the case for the effective gap $1/\xi(t)^{2}$ nor for its derivative,
which attain the finite values

$\displaystyle\lim_{t\to 0^{-}}$
$\displaystyle\xi^{2}(0)=\frac{\Gamma(p)\Gamma(p+1)}{2\pi p^{2p}},$ (105)
$\displaystyle\lim_{t\to 0^{-}}$ $\displaystyle
2\dot{\xi}(0)\xi(0)=\frac{1}{\sqrt{3}}\ ,$ (106)

where $p=1/3$. The finiteness of the results in Eqs. (250) and (106)
corresponds to a vanishing fidelity in Eq. (99). Consequently, the defect
density diverges, see Eq. (98), but the excess energy remains finite

$\displaystyle\lim_{N\to\infty,\delta\to 0}E_{\rm res}(0)\simeq\lim_{t\to
0}\omega(t)N_{\rm exc}(t)\propto\delta^{\frac{1}{3}}.$ (107)

Eq. (107) defines regime (2) and is consistent with the outcome of the Kibble-
Zurek mechanis [202, 185] as well as with the numerical result in Ref. [106].

(a) Half ramp

(b) Full ramp

Figure 10:  Kibble-Zurek mechanism in the fully-connected model. Panel (a)
shows the residual energy as a function of the ramp speed in the case of a
half ramp $t\in[-h_{\rm cr}/\delta,0]$ for $\Lambda=\\{10^{9},3\cdot
10^{7},10^{6},3\cdot 10^{4},10^{3},3\cdot 10,1\\}$ from top to bottom. The
crossover between the Kibble-Zurek scaling (black dashed line) at large
$\Lambda$ and the analytic scaling (gray dashed line) at small $\Lambda$ is
evident. Panel (b) displays the residual energy after a quasi-static drive
obtained by spin-wave theory. The result is obtained within regime (3) and
perfectly reproduces the slow-drive universality numerically found in Ref.
[105]. Each colour represents a different value of $\Lambda=N\delta$ with
$N=2^{9}$ and $N=2^{11}$ (dashed and solid lines in panel b), i.e. the same
values displayed for the exact numerical study in Fig. 9b. The curves at
different sizes perfectly collapse when drawn as a function of the scaling
variables. Moreover, the agreement between the spin-wave theory and the
numerical study for the different values of $\Lambda$ is rather remarkable.

Regime (3) ††margin:  Regime (3) is obtained considering directly the
thermodynamic limit case $\Delta_{N}=0$ and taking the dynamics in $t=h_{\rm
cr}/\delta\approx\infty$ limit, which yields the $\delta$-independent results

$\displaystyle\lim_{t\to\infty}N_{\rm exc}(t)$ $\displaystyle=\frac{1}{3}$
(108) $\displaystyle\lim_{t\to\infty}f(t)$ $\displaystyle=\frac{\sqrt{3}}{2}\
,$ (109)

which characterise the non-adiabatic dynamics as they remain finite in the
$\delta\to 0$ limit. The analytical results in Eqs. (108) and (109) are
universal in the traditional of Kibble-Zurek mechanis result. So, they
faithfully reproduce the slow drive limit $\delta\to 0$ of any dynamical
protocol which crosses the critical point. The universality phenomenon is
analysed in details in Ref. [53]. On a less accurate, but perhaps more
quantitative level, it can be numerically verified that the analytic solution
in D.3 accurately describes any drive $\omega^{\prime}(t_{*})$ such that
$|\omega^{\prime}(\hat{t})-\omega(\hat{t})|^{2}\ll 1$, where $\omega(t)$ is
given in Eq. (101) [54]. It is worth noting that the non-adiabatic regime
described by Eq. (109) and Eq. (108) is profoundly different from the one
described in Ref. [203] for low-dimensional systems. There, the spin-wave
description is applied to the case $\alpha\to\infty$. Then, differently from
our case on has to integrate over the continuous spin-wave spectrum and non-
adiabaticity may arise also for a quench to the critical point due to the flat
density of states of 2d systems.

In summary, dynamical corrections to a quasi-static drive solely depend on the
universal variable $\Lambda=N\delta$. Indeed, the instantaneous minimal gap of
a finite system scales as $\Delta_{N}\propto\Lambda^{-1/3}$, it follows that
the thermodynamic limit ($N\to\infty$) and the adiabatic one ($\delta\to 0$)
do not commute. Accordingly, the observable expectations are universal when
displayed as a function of the universal variable $\Lambda$, see Fig. (10b).
This is in perfect agreement with the numerical findings of Ref. [105, 53].

Summary: Non-adiabatic corrections depend only on the product of the system
size and the adiabatic parameter. This leads to a novel form of universality,
where no Kibble-Zurek scaling is observed in the thermodynamic limit.

####  3.3.4 Full counting statistics of defects

Recently, ††margin:  Defect distribution interest has raised around the
universality of the higher cumulants of the defect statistics following a
quasi-static ramp. In general, the process of defect formation in finite local
systems has been argued to follow a Binomial distribution [204], making the
process of defect formation across a conventional quantum phase transition
akin to the classical process of a coin toss [205]. Approaching the
thermodynamic limit the probability to generate $n$ defects becomes normal and
reads

$\displaystyle P_{\rm local}(n)\approx\frac{1}{\sqrt{2\pi(1-p)\langle
n\rangle}}\exp\left(\frac{(n-\langle n\rangle)^{2}}{2(1-p)\langle
n\rangle}\right)\ ,$ (110)

where the average number of defects follows Kibble-Zurek scaling $\langle
n\rangle\propto\delta^{\frac{d\nu}{1+z\nu}}$ and $p$ is the probability for
the formation of a single defect. The above theory can be exactly verified in
the nearest neighbour transverse field Ising model, whose full counting
statistics can be calculated exactly. While for a finite quench rate $\delta$
all moments of the distributions remain finite in the slow drive limit one
recovers Eq. (110) with $(1-p)=3/\pi^{2}$ and $\langle
n\rangle=\frac{N}{2\pi}\sqrt{\frac{\delta}{2J_{0}}}$. These findings have been
also demonstrated on different quantum computing platforms [206, 207].

The phenomenology of local systems is certainly rich, but does not present any
peculiar features due to quantum fluctuations. Indeed, the same theoretical
framework can be applied to describe the statistics of the defects generated
across a classical and a quantum phase transition [204]. As we have already
argued in previous section, long-range interaction radically alter this
picture as they suppress the defect contribution arising from semi-classical
critical dynamics and promote the single quantum model as the leading source
of non-adiabatic corrections. Interestingly, using the same methodology
employed to derive Eq. (108), one can obtain the full counting defect
statistics. Indeed, the probability to generate $n$ defects in our problem is
just given by the $|c_{n}|^{2}$ coefficient in Eq. (98) [208]. ††margin:
Long-range result Some simple manipulations of Eq. (234) yield

$\displaystyle P_{\rm
LR}(n)\approx{{n+k-1}\choose{n}}\mathrm{sech}(r)^{2k}\tanh(r)^{n},$ (111)

which is a negative binomial distribution. The parameter $\tanh(r)$ defines
the squeezing parameter, see Eq. (87) of the single spin-wave of the system

In light of Eq. (111) the parallel between the local and long-range full
counting statistics becomes rather striking. The negative binomial
distribution described the probability to obtain $n$ failures before a given
(non-deterministic) number $k$ of successes occurs. So, the probability for a
defect (actually a pair since $n\in 2\mathbb{Z}$) to arise in the quantum
long-range system coincides with the probability to observe $n/2$ failures
before the $k$-th success. However, $k=1/2$ in the present problem. Negative
binomials of fractional order are often dubbed Polya distributions and they do
not have any equivalent in classical processes but they naturally emerge in
the defect formation of quantum long-range systems due to the quantum nature
of the problem. This has to be contrasted with the case of local critical
theories where defect formation is a purely classical process [205]. The
implications of these fundings to the quantum thermodynamics of the systems
are discussed in Ref. [209, 208].

Summary: Long-range interactions suppress the creation of defect, leading to a
negative binomial defects distribution, without a classical counterpart.

###  3.4 Dynamical quantum phase transitions - Loschmidt echo

Up to this point, our discussion focused on the most traditional examples of
dynamical critical phenomena, but, recently, experimental advancements in
quantum simulations with cold atoms [210, 211, 212, 213, 214, 215, 22] and
trapped ions [45] raised the interest on novel form of dynamical criticality
[11]. This is the case of dynamical phase transitions. On one side, the name
referred to the study of the out-of-equilibrium behavior of order parameters
[216, 217, 218, 219, 220, 221], which we refer to _dynamical phase transition
(DPT) in the order parameter_ and that will be discussed in details in Section
4.1.2. On the other, a novel form of dynamical criticality was discussed,
where _nonanalytic cusps in the Loschmidt echo_ rate function appear after a
quantum quench [222, 223, 55]. We refer to the latter here as _dynamical
quantum phase transitions (DQPT) in the Loschmidt echo_.

During ††margin:  Loschmidt amplitude and echo a quench, the system initially
prepared in an initial state $\ket{\Psi_{0}}$ is evolved through a time
independent final Hamiltonian $H$, i.e.
$\ket{\Psi(t)}=\exp(-\mathrm{i}Ht)\ket{\Psi_{0}}$. The Loschmidt amplitude
describes the amplitude of the system returning in its initial state at time
$t$ and reads

$\displaystyle\mathcal{G}(t)=\langle\Psi_{0}|\Psi(t)\rangle=\bra{\Psi_{0}}\exp(-\mathrm{i}Ht)\ket{\Psi_{0}}$
(112)

whose expression closely resembles the classical finite-temperature partition
function $Z(\beta)=\mathrm{Tr}\exp(-\beta H)$. The Loschmidt echo is simply
obtained by squaring the amplitude in Eq. (113) yielding

$\displaystyle\mathcal{L}(t)=|\mathcal{G}(t)|^{2}.$ (113)

The Loschmidt amplitude and Loschmidt echo play central roles in the theory of
DQPTs and appear in various contexts in quantum many-body theory [224, 225,
226, 227, 228, 229]. They exhibit a functional dependence on the system size
$N$, and in the limit of large $N$, they can be described by rate functions
that capture their scaling behavior [222, 230]. DQPTs are defined as
nonanalytic behaviors of the Loschmidt amplitude or Loschmidt echo as a
function of time. They can be considered as phase transitions in time,
analogous to equilibrium phase transitions associated with nonanalytic
structures of the free energy [231, 222]. A DQPT is characterized by a sudden
qualitative change in the dynamics of the system, typically accompanied by a
kink or nonanalyticity in the rate function of the Loschmidt amplitude or
Loschmidt echo. This nonanalytic behavior can vary depending on the system and
dimensionality, including power-law singularities, logarithmic singularities,
and other forms [222, 232, 233].

Theoretical evidences of DQPTs in the Loschmidt echo return rate was found in
numerous quantum systems [222, 223, 234, 235, 232, 236, 237, 238, 239, 240]
and connected with the behaviour of different local observables [241],
including different definitions of the order parameter [242, 243]. Given the
many successful experimental realizations of this kind of DQPTs, especially in
trapped ion systems with long-range interactions [244, 21], it is not
surprising that long-range interacting models were also a privileged tool for
the theoretical characterization of DQPTs [245, 242, 246, 247, 248, 249, 232,
189, 250, 243].

Summary: Dynamical quantum phase transitions (DQPTs) in the Loschmidt echo are
defined by non-analytic behavior of this quantity as a function of time.

####  3.4.1 Spherical spin-wave theory

In the following, we are going to show how the emergence of DQPTs in the long-
range Ising model can be captured by the study of the harmonic Hamiltonian in
Eq. (46). However, most DQPTs occur when a parameter of the Hamiltonian is
quenched across an underlying equilibrium phase transition and the equilibrium
spin-wave theory is not capable to capture both side of the transition within
the same formalism. It is worth noting that there are cases where DQPTs can
arise independently of conventional phase transitions [235, 251, 252, 236,
245].

††margin:  The spherical model

In order to capture a sudden quench across the dynamical critical point we are
going to consider a slightly modified version of the linear theory Eq. (46),
namely the quantum spherical model [253]. The idea consists in modifying the
dispersion of the linear spin-waves as where we have used the transformation

$\displaystyle\omega_{k}=\sqrt{h(\mu+2J_{0}f_{k}(\alpha))},$ (114)

where once again we have set $\gamma=1$. While in Eq. (47) the excitation
spectrum only depends on the magnetic field and the coupling, the novel
dispersion relation contains an addition parameter $\mu$. This has to be
calculated self-consistently by imposing a constraint on the spin-wave
potential energy, namely

$\displaystyle\sum_{k}\frac{\left\langle\left(\beta^{\dagger}_{k}+\beta_{k}\right)^{2}\right\rangle}{2h\omega_{k}}=\frac{N}{4}\
.$ (115)

The spherical constraint in Eq. (115) may induce a quantum critical point in
the quadratic model at a field value $h_{\rm cr}^{\rm sph}$ depending on the
features of the function $f_{k}(\alpha)$. Below the critical field strength
$h_{\rm cr}^{\rm sph}$ the constraint condition in Eq. (115) causes the spin-
wave to form a condensate state at $k=0$, whose formation follows the mean
field universality. In its classical version, the spherical model has been
introduced to mimic the finite temperature free-energy of $O(n)$-symmetric
spin systems in the $n\rightarrow\infty$ limit [254].

Summary: The spherical model corresponds to the low-energy long-range
Hamiltonian equipped with an additional parameter, determined by a self-
consistent equation. Its critical behavior belongs to the mean-field
universality.

####  3.4.2 Step approximation

The DQPT occurs in the spherical model following a sudden quench. Thus, we
have to consider the average in Eq. (115) over a time dependent state
$\ket{\Psi(t)}=\prod_{k}\ket{\psi_{k,0}(t)}$, where $\ket{\psi_{k,0}(t)}$ is
the single-spin wave state given by Eq. (229) with $n=0$. The explicit
calculation leads to the dynamical constraint equation ††margin:  State
evolution after a quench

$\displaystyle\frac{h}{N}\sum_{k}\frac{\xi_{k}^{2}(t)}{2}=1/4,$ (116)

where $\xi_{k}^{2}(t)$ is the solution of the Ermakov Eq. (88). Similarly to
Sec. 3.1, we are going to consider a sudden quench of the ferromagnetic
coupling $J_{0}^{\rm i}\to J_{0}^{\rm f}$, but now the final coupling value
$J_{0}^{\rm f}$ shall drive the system across the quantum phase boundary. In
this case, a solution of the differential Eq. (88) together with the dynamical
constraint in Eq. (116) is rather complicated and, up to our knowledge, has
not been attempted yet. On the other hand, a convenient simplifying assumption
consists in assuming that, as the ferromagnetic coupling $J_{0}$ is quenched,
the parameter $\mu$ also undergoes a discontinues jump between two constant
values $\mu_{0}$ at $t<0$ and $\mu_{f}$ for $t\geq 0$. This procedure goes
under the name of _step approximation_ and has been introduced in Refs. [163,
167]. In order for this approximation to be sensible, one should choose the
final value $\mu_{f}$ in order to reproduced its expected long-time
(equilibrated) value. As long as $\alpha>d$ the system can be safely assumed
to equilibrate at long time as shown in Ref. [167, 109],

Within the framework of the step-approximation, the frequency suddenly change
from their initial value $\omega_{k,i}$ to a final value $\omega_{k,f}$,
leading to the following sudden ††margin:  Solution quench solution for the
effective length of each spin-wave

$\displaystyle\xi_{k}(t)=\sqrt{\frac{1+\epsilon_{k}\sin^{2}(\omega_{k,f}t)}{2\omega_{k,i}}}$
(117)

with the quench parameter
$\epsilon_{k}=\left(\frac{\omega_{k,i}}{\omega_{k,f}}\right)^{2}-1$. In the
thermodynamic limit $N\rightarrow\infty$ the sum in Eq. (116) can be turned
into an integral. Then, taking into account the explicit solution in Eq. (117)
one obtains

$\displaystyle\int\\!\frac{dk}{2\pi}\frac{h}{2\omega_{k,i}}\left[\frac{\epsilon_{k}}{2}(1-\cos
2\omega_{k,f}t)\right]=0.$ (118)

Latter equation cannot be fulfilled at all times due to the oscillatory time
term. Yet, in the limit $t\to\infty$ the dephasing between the different modes
washes away the time dependence in Eq. (118), making the solution in Eq. (117)
exact also for the constraint problem as long as the final value of $\mu$ is
chosen in order to satisfy the following expression

$\displaystyle\int\frac{dk}{2\pi}\frac{\epsilon_{k}}{\omega_{k,i}}=0.$ (119)

This implicit equation determines the long-time asymptotic value of $\mu_{f}$
through the $\mu$ dependence of $\omega_{k,f}$. The consistency of the
equilibration assumption and, overall, of the step approximation can be
verified by the inspection of the numerical solution of the exact problem, see
Ref. [109].

Eq. (119) can be used to determine the dynamical critical coupling $J_{0}^{\rm
c,dyn}$ at which the dynamical excitation become gapless and the constraint
parameter $\mu_{f}$ in Eq. (119) approaches its critical value $\mu_{c}$.
Using these definitions, Eq. (119) can be rewritten as

$\displaystyle\frac{1}{2}=\sqrt{h}\int\frac{dk}{2\pi}\frac{\sqrt{2\mu_{0}+2J_{0}^{i}f_{k}(\alpha)}}{2\mu_{c}+2J_{0}^{\rm
c,dyn}f_{k}(\alpha)},$ (120)

where $\mu_{c}$ is the equilibrium critical value. The existence of a finite
value $J_{0}^{\rm c,dyn}$ satisfying Eq. (120) depends both on the parameters
$h,\mu_{0}$ and on the value of $\sigma$. The dynamical phase diagram of the
model is reported in Ref. [109]. In the present section, for the sake of
simplicity, we are going to assume that $J_{0}^{i}$ lies above its equilibrium
critical value $J_{0}^{i}>J_{0}^{\rm c}$, i.e. in the condensate phase, and
consider the case $\alpha<d+2$.

Given the quadratic nature of the spherical model, the overlap function can be
calculated analytically

$\displaystyle\begin{split}\mathcal{G}(t)=&\prod_{k}\Biggl{\\{}(8\omega_{k,i})^{1/4}e^{-i\varphi_{k}(t)}\biggl{(}2\omega_{k,i}\xi_{k}(t)+\frac{1}{\xi_{k}(t)}-\mathrm{i}2\dot{\xi}_{k}(t)\biggr{)}^{-1/2}\Biggr{\\}}\
,\end{split}$ (121)

where the unessential phase $\varphi(t)$ is defined below Eq. (88). The
††margin:  The Loschmidt echo rate Loschmidt echo rate function is obtained by
taking the logarithm of the squared overlap, yielding

$\displaystyle r(t)=-\lim_{N\rightarrow\infty}\frac{1}{N}\log\lvert
G(t)\rvert^{2}=-\log 2+\int\frac{dk}{2\pi}\log\lvert X_{k}(t)\rvert,$ (122)

where

$\displaystyle
X_{k}(t)=\frac{1}{\sqrt{8\omega_{k,i}}}\biggl{(}2\omega_{k,i}\xi_{k}(t)+\frac{1}{\xi_{k}(t)}-\mathrm{i}2\dot{\xi}_{k}(t)\biggr{)}$
(123)

As long as $\Omega_{k,i}$ is gapped, the expression in Eq. (123) remains
smooth and no cusp appears at finite time for the rate function defined in Eq.
(122). The non-analytic cusps characterizing DQPTs will only appear for a
sudden quench from the broken phase where $\Omega_{k,i}$ is gapless. This
result demonstrates how one can observe and characterize DQPTs by just
analyzing the quasi-particle spectrum. Therefore, as already mentioned, we are
going to consider a sudden quench of the coupling $J_{0}^{\rm i}\to J_{0}^{\rm
f}$, with the initial coupling within the ferromagnetic phase and the final
one above the dynamical critical threshold $J_{0}^{\rm f,dyn}$. Since the
dynamics is initiated in the broken phase, a complete treatment of the problem
shall also include the evolution of the classical mode representing the
condensate fraction of spin-waves as it was done in Ref. [238]. However, in
the present description we have discarded this contribution as it is not
necessary to observe the DQPTs.

The signature of DQPT in the Loschmidt echo dynamics is reported in Fig. 11,
the Loschmidt echo is shown for a quench from $J_{0}^{i}=2J_{0}^{c}$ to
$J_{0}^{f}=J_{0}^{c}/2$ for different values of $\sigma$. The rate function
clearly shows non-analyticities at the critical times:

$\displaystyle t^{*}_{m}=\frac{m\pi}{\omega_{k,\mathrm{f}}}\quad
m\in\mathbb{N}$ (124)

which appear due to logarithmic divergences in the integrand in Eq. (122).
Since the critical time scale is set by the post-quench gap, we do not expect
to see nonanalytic cusps in the Loschmidt echo for a quench into the gapless
phase, as previously mentioned.

Figure 11: DPT in the spherical model. The Loschmidt echo rate function of the
spin-wave theory (in presence of the spherical constraint) after a sudden
quench of the ferromagnetic coupling $J_{0}^{i}=2J_{0}^{c}$ to
$J_{0}^{f}=J_{0}^{c}/2$. The cusps in the rate function are clearly evident
and the second derivative is divergent since the dynamics considered is for
$\sigma=1.5$, see the discussion in the text. Figure adapted from Ref. [109].

Upon differentiating Eq. (122) with respect to time $n$ times, we encounter
terms that are proportional to $1/\omega_{k,i}^{n}$ when $t=t^{*}_{m}$. These
terms diverge as $k^{-n(\alpha-d)/2}$. However, since the integration over $k$
still has to be performed, the $n$-th derivative of the rate function only
diverges if $n(\alpha-d)/2>1$, or equivalently, $n(\alpha-d)>2$. Thus, the
smaller $\sigma$ the larger the order of the derivative for which the cusps
are expected. This analysis holds true throughout the entire region where
$d<\alpha<d+2$. On the other hand, it does not apply when $\alpha\geq d+2$
since there is no gapless phase to initiate the calculation. The discussion
presented above offers further evidence that the emergence of nonanalytic
cusps is not solely a feature of the step approximation. Instead, it is a
consequence of the initial conditions, specifically starting in the gapless
phase, as well as the specific form of the function $\xi_{k=0}(t)$ and its
time derivatives, which remain unchanged in the exact calculation.
Furthermore, the structure of these cusps remains unaltered even in the long-
time limit when $\mu(t)$ reaches equilibrium and the step approximation
becomes exact.

Summary: Loschmidt echo DQPTs in the spherical model appear for quenches from
the broken phase, as results from the exact solution.

####  3.4.3 Strong long-range regime

The aforementioned analysis cannot be extended to the regime $\alpha<d$ as the
system does not equilibrate and the assumptions at the root of the step
approximation outlined in Sec. 3.4.2 fail. Moreover, due to the gapped nature
of the spectrum the importance of the condensate motion of the system at
$\alpha<d$ is more prominent and cannot be easily discarded. Several numerical
simulations and analytical arguments have been used to show the existence of
the DQPT in the Loschmidt echo also for the actual Ising Hamiltonian for
different values of $\alpha$ [234, 255, 248, 247]. In particular, extensive
numerical studies have been devoted to investigate the connection between the
DQPTs occurring in the Loschmidt echo and the DPTs defined via the dynamical
scaling of the order parameter [242, 246], see also Section 4.1.2. Also, the
relation between the two different notions and the quasi-particle properties
of the model have been largely investigated, but mostly close to the local
limit [249, 250, 243]. In the next section, we are going to introduce the
dynamical Holstein-Primakoff transformation, which represents the proper
formalism to describe the motion of spin-wave coupled to the classical order
parameter and their feedback effect. However, we are only going to use it to
describe DPTs in the order parameter leaving aside further comments on
singularities of the Loschmidt echo.

Summary: In the strong long-range regime, the leading effect in the Loschmidt
echo DQPT comes from the “high-energy” dynamics of the order parameter, and it
can be related to other forms of dynamical criticality.

##  4 Dynamics in highly excited states

In this Section we will discuss the treatment of out-of-equilibrium dynamics
involving arbitrarily high energy initial states, as in standard quantum
quench protocols.

We will begin in Sec. 4.1 by reviewing mean-field dynamical phenomena for
$\alpha=0$. In Sec. 2.3 we showed that the fully-connected limit of quantum
spin systems reduces to the physics of a single collective degree of freedom.
In this Subsection we will discuss how this statement applies to non-
equilibrium dynamics as well.

Secondly, in Sec. 4.2, we will discuss how finite-range interactions with
$\alpha>0$ affect mean-field dynamical phenomena. Long-range interacting
system can be formally viewed as a perturbation of the mean-field limit, as
reviewed in Sec. 2.4. The perturbation term couples the collective degree of
freedom to many spin-fluctuation modes with various wavelengths, resulting in
a genuine many-body problem, which can be addressed via the non-equilibrium
spin-wave formalism developed in Refs. [56, 57]. The coupling strength to mode
$\mathbf{k}$ strongly depends on the interaction range governed by the
exponent $\alpha$, as encoded by the function $f_{\mathbf{k}}(\alpha)$. As a
result of this tunable decoupling, long-range interactions give rise to
anomalous non-equilibrium many-body phenomena. Collective spin ordering is
remarkably resilient out of equilibrium, generating long pre-thermal stages of
dynamics characterized by long-lived oscillating collective spin polarization.
This behaviour was observed in numerical simulations performed with a range of
techniques [242, 125, 127] and theoretically understood via the aforementioned
approach [56, 57]. This analysis shows that the duration of the prethermal
stage increases as $\alpha$ is decreased, and diverges with the system size
when $\alpha<d$ [256, 58].

###  4.1 Quench dynamics of fully-connected spin systems ($\alpha=0$)

This Section is devoted to the non-equilibrium dynamics of fully connected
spin systems. We study the time-evolution starting from ground states
$\ket{\psi_{0}}$ of a pre-quench Hamiltonian $\hat{H}(h_{0})$ evolving with a
different post-quench Hamiltonian $\hat{H}(h_{f})$. For the sake of
definiteness we will mostly consider the Ising model, Eq. (1) with $\gamma=1$,
and quenches in the transverse field from $h_{0}$ to $h_{f}$. As described in
Section 2.3.2, ground states of long-range Hamiltonians as (1) can be though
as _coherent states_ in some directions. Their dynamical behavior is
determined by a classical mean-field description emerging in the thermodynamic
limit when initialized in fully polarized states, as described in Section
4.1.1. The resulting dynamics of collective observables can give rise to new
forms of dynamical criticality, such as dynamical phase transitions, discussed
in Section 4.1.2. The semiclassical framework also allows us to describe the
growth of quantum fluctuations, which coincides with the flow of linearized
shifts around classical trajectories and is thus related to the standard
quantifiers of classical chaos, reviewed in Section 4.1.3. When the quantum
fluctuations become comparable to the typical length of the phase space, this
description breaks down, defining an Ehrenfest time that diverges with $N$ for
this class of systems. Remarkably, such semiclassical framework grasps crucial
aspects of quantum dynamics, such as the dynamics of scrambling or
entanglement as reviewed in Section 4.1.4 and 4.1.5 respectively.

####  4.1.1 Mean-field classical limit

The dynamics of a system with unbroken full permutational symmetry take place
in the totally-symmetric subspace (TSS) of the many-body Hilbert space
simultaneously invariant under all permutations. Such dynamics is amenable to
an exact representation in terms of few collective degrees of freedom,
characterized by an effective Planck constant $\hbar_{\rm eff}\sim 1/N$
suppressed with system size [219, 257]. We refer to B for a general
discussion.

For systems of interacting quantum spins the limiting semiclassical
description may be formulated more directly and intuitively in terms of states
with maximal collective spin $S=Ns$ — the so-called Dicke manifold. As
discussed in Sec. 2.3.1, the collective spin approaches a classical limit for
large $N$.131313For $s=1/2$, the TSS coincides with the Dicke manifold. For
larger $s$ there are more permutationally invariant states with lower $S$
(${\rm dimTSS}\sim N^{2s}$ for large $N$). However, for Hamiltonians without
spin self-interactions, one may always consider dynamics within the Dicke
manifold. For the infinite-range XY Hamiltonian in Eq. (5) the classical limit
$\hat{H}_{\alpha=0}/N\to\mathcal{H}_{\rm cl}$ is given by Eq. (11) in Sec.
2.3.1, where we discussed equilibrium properties. In this Section, we will use
the same approach to discuss _out-of-equilibrium_ properties. For
definiteness, throughout this Section, we will set $\gamma=1$ (quantum Ising
model).

The ††margin:  Classical dynamics non-equilibrium evolution
$\braket{\hat{\vec{S}}(t)}/N$ generated by a sudden change (“quench”) of a
Hamiltonian parameter is described by a classical trajectory
$\vec{\mathcal{S}}(t)$ on the unit sphere governed by
$\mathcal{H}_{\text{cl}}$, i.e.,

$\dot{\vec{\mathcal{S}}}=\big{\\{}\vec{\mathcal{S}},\mathcal{H}_{\text{cl}}\big{\\}}\
,$ (125)

with the canonical Poisson brackets
$\\{\mathcal{S}^{\mu},\mathcal{S}^{\nu}\\}=\epsilon_{\mu\nu\rho}\mathcal{S}^{\rho}$.
Evolution can be recast in terms of the spherical angles
$\theta(t),\,\phi(t)$. In the case of the Hamiltonian (5), the non-linear
precession of the collective spin is described by the classical equations of
motion141414For convenience, we rescale time by a factor $s$

$\displaystyle\begin{dcases}\dot{\theta}=2J_{0}\sin\theta\cos\phi\sin\phi\,,\\\
\dot{\phi}=-h+2J_{0}\cos\theta\cos^{2}\phi\,.\end{dcases}$ (126)

As the Hamiltonian governs a single degree of freedom, the classical limit is
trivially integrable and characterized by regular periodic trajectories in
phase space. Such behaviour corresponds to persistent spin oscillations after
a quench, whose period depends on the initial state. For $|h|<2J_{0}$, the
phase-space also features a separatrix with a diverging classical period,
terminating at the saddle point $\theta=0$ and characterized by an exponential
instability rate

$\lambda=\sqrt{h(J_{0}-h)}\ $ (127)

(i.e. the eigenvalue of the stability matrix at the saddle point).

While ††margin:  Chaotic dynamics in other all-to-all spin models such fully-
connected spin models generically exhibit periodic orbits, semiclassical
chaotic behavior can occur in a number of relevant situations. A standard
example comes from introducing time-dependent driving, thus breaking energy
conservation: The quantum kicked top [258, 259], corresponding to a step-wise
driving protocol applied to the model above, gives a paradigmatic regular-to-
chaotic crossover as a function of the driving parameters. Another source of
chaoticity comes from coupling the spins to other degrees of freedom, such as
a cavity mode, which gives rise to the textbook Dicke model [260, 261].
Finally, chaotic behaviour can arise from self-interactions of higher spins
$s>1/2$, which are generally described by $n=2s>1$ collective degrees of
freedom: In the absence of additional symmetries, self-interactions will break
classical integrability. These extended possibilities can be addressed with
the method summarized in B. In the rest of this Section, we will refer to them
when discussing the impact of chaos on the quantum dynamics of fully-connected
systems.

Summary: In the thermodynamic limit the quench dynamics of fully-connected
spin systems is described by classical periodic trajectories of a single
collective degree of freedom.

####  4.1.2 Dynamical phase transitions - Order parameter

Figure 12: Equilibrium configurations and possible instances of non-
equilibrium dynamics in the fully-connected quantum Ising model. (a-c)
Pictorial representation on the Bloch sphere of the collective spin for the
post-quench Hamiltonian. (a-b) For $h_{f}<|J|$, the energy possesses two
minima characterized by non-vanishing, opposite magnetizations along $x$. (c)
For $h_{f}>|J|$, the system is paramagnetic with a single equilibrium
configuration in the direction of the field. Initial fully polarized states at
$t=0$ are pictured as a point on the Bloch sphere, surrounded by a small grey
circle representing their transverse quantum fluctuations. Labels (a1-3)*
represent possible instances of such initial conditions. (a1-a3) Semiclassical
phase portrait of the ferromagnetic post-quench Hamiltonian, where the initial
states move along a nontrivial nonequilibrium trajectory, corresponding to the
initial conditions (a1-3)* respectively. (a4-a5) Associated dynamics of the
classical magnetization. Labels (a) refer to initial ferromagnetic initial
states $h_{0}<|J|$. Their time-evolution is characterized by ferromagnetic
periodic (green) trajectories [see (a1) and (a4)] or paramagnetic (blue) ones
[see (a3) and (a6)] with $\overline{S_{x}(t)}\neq 0$ and
$\overline{S_{x}(t)}=0$, respectively. These are separated by the unstable
(red) trajectory occurring at $h_{\rm dcr}$ [see (a2) and (a5)]. Labels (b-c)
refer to an initial paramagnetic state $h_{0}=\infty$ evolved with two
different Hamiltonians: (b) quench performed to a ferromagnetic Hamiltonian
$h_{f}<|J|$, the initial state lies on the unstable trajectory, (c) quench
performed to a different $h_{f}>|J|$ paramagnetic configuration. Images
adapted from Ref.[57, 58].

The ††margin:  Magnetic ordering in the asymptotic state non-equilibrium
evolution described above may or may not result in collective spin ordering at
long times. An abrupt change of dynamical ordering properties as a function of
driving control parameters is referred to as a dynamical phase transition
(DPT) [262, 218, 219, 263, 264, 265, 257, 266, 267, 268, 245, 248, 221]. In
particular, when a system is quenched from a symmetric state across the
equilibrium critical point, dynamical scaling properties associated with aging
or coarsening may appear [167, 269, 270, 164]. Conversely, when a system
undergoes a sudden quench from a broken-symmetry state, the resulting out-of-
equilibrium dynamics may display two different phases. One can define a non-
equilibrium order parameter by time-averaging the corresponding equilibrium
order parameter. This quantity may be vanishing or not depending on whether
the symmetry is dynamically restored after the quench. The associated
_dynamical critical point_ is believed to have a universal character. A
special interest was placed on systems that fail to rapidly approach thermal
equilibrium after the quench, as their dynamical universality may have no
equilibrium counterpart [257, 263].

Fully-connected ††margin:  Fully connected systems spin systems provide the
simplest instance of genuinely dynamical phase transitions. To illustrate this
we consider the infinite-range quantum Ising model [Eq.(5) with $\gamma=1$].
The character of non-equilibrium dynamics is encoded in the classical
trajectories of the collective spin, which may have paramagnetic or
ferromagnetic character. Here one studies quenches in the transverse field
$h_{0}\to h_{f}$ for which DPTs have been extensively studied [245, 219, 248,
271]. The two non-equilibrium phases are distinguished by the time average of
the equilibrium order parameter,

$\overline{\mathcal{S}^{x}}=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}dt\frac{\braket{\hat{S}^{x}(t)}}{Ns}$
(128)

which serves a non-equilibrium order parameter: It is finite
$\overline{\mathcal{S}^{x}}\neq 0$ for shallow quenches in the dynamical
ferromagnetic phase and it vanishes abruptly $\overline{\mathcal{S}^{x}}=0$ at
the dynamical critical point $h_{\rm dcr}=(h_{0}+J_{0})/2$, associated with
the critical trajectory corresponding to the phase-space separatrix; for
deeper quenches $h_{f}>h_{\rm dcr}$ the system lies in the dynamical
paramagnetic phase. See Fig. 12(a-b) for an illustration. This kind of DPT has
been realized experimentally with cold atoms in optical cavities [22] or in a
superconducting quantum simulator [272].

The spectral counterparts of these DPTs are given by excited-state quantum
phase transitions (ESQPT) [273, 274, 275, 276, 277, 278, 279]. This
corresponds to singularities of the density of states at some finite energy
density which distinguishes eigenstates with ferromagnetic nature from those
with paramagnetic nature: See Ref. [280] for a recent review.

The notion of DPT discussed here is in general distinct from that of DQPT
discussed in Sec. 3.4, and therefore they may even not occur concomitantly in
the same model. However, a connection has been pointed out whenever both
phenomena are present [242] (see also Refs. [238, 239, 57]). Below, in Section
4.2.3, we will discuss how this out-of-equilibrium phenomenon is affected by
decreasing the interaction range.

Let us mention that, due to the intrinsic semiclassical nature of dynamics in
this class of models, observables can be efficiently simulated using phase-
space numerical techniques [281, 282, 283], such as the Truncated Wigner
Approximation (TWA) [284, 285, 282] or its discrete [286, 129] or clustered
[287] versions. These methods have been intensively used also to explore
dynamics with finite $\alpha$ interactions, which we discuss in the next
section [129, 126, 288, 289, 290, 291, 292, 293, 294, 295, 296].

Summary: The asymptotic state after quantum quenches in fully connected
systems may or may not display ordering. These dynamical phase transitions are
associated with a separatrix in the classical phase-space.

####  4.1.3 Semiclassical dynamics of quantum fluctuations

The classical description of dynamics outlined above is exact in the
thermodynamic limit $N\to\infty$. In finite systems, however, it has a limited
time scale $T_{\text{Ehr}}(N)$ of validity, known as Ehrenfest time scale: At
long times $t\gtrsim T_{\text{Ehr}}(N)$, quantum fluctuations around the
classical limit will dominate the behavior of time-dependent local observables
and entanglement quantifiers. $T_{\text{Ehr}}$ can be estimated as the time at
which the size of quantum fluctuations becomes comparable with a
characteristic phase-space scale. This depends on the initial state and on the
nature of the underlying classical dynamics. In this Subsection, we discuss
the semiclassical dynamics of quantum fluctuations.

To ††margin:  Time-dependent Holstein-Primakoff compute the evolution of
quantum spin fluctuations it is convenient to generalize the Holstein-
Primakoff approach introduced in Sec. 2.4, above, to the non-equilibrium
context [56, 57]. When the system is driven out of equilibrium, the direction
of the collective spin configuration [parametrized by $\theta(t)$ and
$\phi(t)$] moves along the corresponding classical trajectory on the unit
sphere. We thus let the adapted frame of reference
$(\mathbf{X},\mathbf{Y},\mathbf{Z})$ in Eq. (12) vary in time, in such a way
that the $\mathbb{Z}$-axis follows the evolution of
$\braket{\mathbb{\hat{S}}(t)}\propto\mathbf{Z}(t)$. This way, the collective
spin components along $\mathbf{X}$ and $\mathbf{Y}$ are associated with
quantum fluctuations and will be mapped to canonical bosonic variables.

The time-dependent spin rotation described above is implemented by the time-
dependent unitary operator

$\hat{V}(\theta(t),\phi(t))=e^{-i\phi(t)\,\hat{S}^{z}}\,e^{-i\theta(t)\,\hat{S}^{y}}.$
(129)

where the time-dependence of the angles is for the moment unspecified. The
Heisenberg equations for spin components $\hat{S}^{\mu}$ with $\mu=X,Y,Z$ in
the mobile frame will then read

$\frac{d}{dt}\,\hat{S}^{\mu}=\frac{d}{dt}\,V\hat{S}^{\mu}V^{\dagger}=\frac{1}{i}[\hat{S}^{\mu},\tilde{H}]\,,\qquad\text{where}\quad\tilde{H}(t)\equiv\hat{V}\,\hat{H}\,\hat{V}^{\dagger}+i\hat{V}\dot{\hat{V}}^{\dagger}\,.$
(130)

The effective time-dependent Hamiltonian $\tilde{H}(t)$ includes inertial
forces arising from the time dependence of $\hat{V}$. A direct calculation
shows

$i\hat{V}\dot{\hat{V}}^{\dagger}=-\,\vec{\omega}(t)\cdot\hat{\vec{S}}\qquad\text{with}\quad\vec{\omega}(t)=\left(-\sin\theta\,\dot{\phi},\dot{\theta},\cos\theta\,\dot{\phi}\right).$
(131)

The time-dependent Hamiltonian $\tilde{H}(t)$ is then transformed to a bosonic
Hamiltonian via the Holstein-Primakoff transformation, cf. Eq. (19). This
yields an expression of the form151515Here we rescaled time by a factor $s$
for convenience.

$\begin{split}\tilde{H}(t)\approx&+(Ns)^{1}\;\mathcal{E}\left(\theta(t),\phi(t)\right)\\\
&+(Ns)^{1/2}\;\Big{(}\tilde{h}^{(1)}_{Q}(t)\hat{q}+\tilde{h}^{(1)}_{P}(t)\hat{p}\Big{)}\\\
&+(Ns)^{0}\left(\tilde{h}^{(2)}_{QQ}(t)\,\frac{\hat{q}^{2}}{2}+\tilde{h}^{(2)}_{PP}(t)\,\frac{\hat{p}^{2}}{2}+\tilde{h}^{(2)}_{QP}(t)\,\frac{\hat{q}\hat{p}+\hat{p}\hat{q}}{2}\right)\\\
&+\mathcal{O}\Big{(}(Ns)^{-1/2}).\end{split}$ (132)

Compared to the “static” rotated-frame Hamiltonian (obtained by just rotating
the spins and mapping to bosons) [see e.g. Eq. (20)], the additional inertial
Hamiltonian modifies the linear terms as $\tilde{h}^{(1)}_{Q}(t)\equiv
h^{(1)}_{Q}(\theta(t),\phi(t))+\sin\theta(t)\;\dot{\phi}(t)$ and
$\tilde{h}^{(1)}_{P}(t)\equiv h^{(1)}_{P}(\theta(t),\phi(t))-\dot{\theta}(t)$,
while quadratic ones are modified as $\tilde{h}^{(2)}_{QQ,PP}(t)\equiv
h^{(2)}_{QQ,PP}\big{(}\theta(t),\phi(t)\big{)}-\cos\theta(t)\;\dot{\phi}(t)$
and $\tilde{h}^{(2)}_{QP}(t)\equiv
h^{(2)}_{QP}\big{(}\theta(t),\phi(t)\big{)}$.

The evolution of $\theta(t)$ and $\phi(t)$ is fixed by the vanishing of the
linear terms $\tilde{h}^{(1)}(t)$, ensuring
$\braket{\hat{S}^{X}(t)}=\braket{\hat{S}^{Y}(t)}=0$. This yields the classical
mean-field equations of motion governed by $\mathcal{H}_{\text{cl}}$, i.e. Eq.
(126) for our model.

On the other hand, the number of collective excitations
$\hat{n}_{0}=(\hat{q}^{2}+\hat{p}^{2}-1)/2$ [see e.g. Eq. (19)] non-trivially
evolves in time. Its dynamics are governed by the time-dependent quadratic
Hamiltonian parametrized by $\tilde{h}^{(2)}(t)$ above. In order to evaluate
them, one computes the Heisenberg equations of motion ††margin:  Dynamics of
quantum fluctuations

$\begin{dcases}\dot{\hat{q}}=+\widetilde{h}^{(2)}_{QP}(t)\;\hat{q}+\widetilde{h}^{(2)}_{PP}(t)\;\hat{p}\\\
\dot{\hat{p}}=-\widetilde{h}^{(2)}_{QQ}(t)\;\hat{q}-\widetilde{h}^{(2)}_{QP}(t)\;\hat{p}\end{dcases}\
,$ (133)

with solution $\begin{pmatrix}\hat{q}(t)\\\
\hat{p}(t)\end{pmatrix}=U(t)\begin{pmatrix}\hat{q}(0)\\\
\hat{p}(0)\end{pmatrix}$, where the $2\times 2$ propagator $U(t)$ can be
formally written as the time-ordered exponential of the matrix defined by the
right-hand side of Eq. (133). One can collect the dynamical fluctuations (or
“correlations”) $G^{QQ}(t)\equiv\langle\hat{q}^{2}(t)\rangle$,
$G^{PP}(t)\equiv\langle\hat{p}^{2}(t)\rangle$ and
$G^{QP}(t)\equiv\frac{\langle\hat{q}(t)\hat{p}(t)+\hat{p}(t)\hat{q}(t)\rangle}{2}$
in the $2\times 2$ _correlation matrix_

$G(t)=\begin{pmatrix}G^{QQ}(t)&G^{QP}(t)\\\
G^{QP}(t)&G^{PP}(t)\end{pmatrix}=U(t)\,G(t=0)\,U^{T}(t)\,.$ (134)

The number of dynamically generated excitations can be expressed as

$\langle\hat{n}_{0}(t)\rangle=\frac{G^{QQ}(t)+G^{PP}(t)-1}{2}=\frac{1}{2}\text{Tr}\bigg{[}G(t)-\frac{\mathbb{1}}{2}\bigg{]}\
.$ (135)

Note that $\det G(t)\equiv 1/4$, which is an exact property of _pure_ Gaussian
states preserved by Hamiltonian evolution. For our fully-connected Ising
model, the equations of motion for the correlation matrix read

$\displaystyle\begin{dcases}\dot{G}^{QQ}=2J_{0}\cos\theta\sin\phi\cos\phi\,G^{QQ}+2J_{0}\left(\cos^{2}\phi-\sin^{2}\phi\right)\,G^{QP}\\\
\dot{G}^{PP}=-2J_{0}\cos\theta\sin\phi\cos\phi\,G^{PP}-2J_{0}\cos^{2}\phi\sin^{2}\theta\,G^{PQ}\\\
\dot{G}^{PQ}=-J_{0}\cos^{2}\phi\sin^{2}\theta\,G^{QQ}+J_{0}\left(\cos^{2}\phi-\sin^{2}\phi\right)\,G^{PP}\\\
\end{dcases}\ .$ (136)

Crucially, ††margin:  Relation to the linearized classical flow because we
obtained these equations by expanding the Hamiltonian in powers of
$\hbar_{\text{eff}}$ and because classical and quantum evolution generated by
quadratic Hamiltonians coincide, _the semiclassical dynamics of quantum
fluctuations_ — characterized by the time-dependent correlation matrix $G(t)$
— _obeys the same equation of motion as the linearized flow of displacements
from the classical trajectories._ This statement actually applies to arbitrary
semiclassical systems with $n$ degrees of freedom, where the correlation
matrix $G(t)$ of the quantum fluctuations becomes a $2n\times 2n$ matrix. We
refer to Ref. [297] or B for a complete discussion. The correlation matrix
$G(t)$ is equivalent to the monodromy matrix whose eigenvalues define _the
finite-time classical Lyapunov spectrum $\\{\lambda_{k}(t)\\}$_ [298].

When the classical dynamics is integrable, nearby initial conditions
generically separate linearly in time, as it becomes manifest via action-angle
variables [58]. Thus, the temporal growth of the quantum correlations is
polynomial, $\langle\hat{n}_{0}(t)\rangle\sim t^{2}\ $. Isolated unstable
trajectories like the separatrix discussed in Sec. 4.1.2 are characterized by
exponential sensitivity, and hence $\langle n_{0}(t)\rangle\sim e^{2\lambda
t}$, where $\lambda$ is the largest eigenvalue of the saddle point that
controls the instability. The asymptotic growth also depends on the initial
conditions for systems with a mixed regular-chaotic phase space, e.g.
resulting from integrability breaking within a Kolmogorov-Arnold-Moser
scenario [298]. On the other hand, in systems with fully developed chaos in
phase space, the Lyapunov spectrum is uniform and nonvanishing. This implies
an asymptotic exponential growth of quantum fluctuations,
$\langle\hat{n}_{0}(t)\rangle\sim e^{2\lambda t}\ $. The classification is
concluded by the case of stable equilibrium configurations, the linearized
dynamics of which are equivalent to that of coupled harmonic oscillators.
Accordingly, all the quantities of interest perform bounded (periodic or
quasiperiodic) oscillations. This classification is summarized in the first
row of Table 1.

Classical trajectory | Stable | Regular | Chaotic
---|---|---|---
|  |  | (Unstable)
Collective fluctuations | oscillations | $t^{2}$ | $e^{2\lambda t}$
Ehrenfest time scale | $\mathcal{O}(\sqrt{N})$ | $\mathcal{O}(\sqrt{N})$ | $\mathcal{O}(\ln N)$
entanglement entropy | oscillations | $\ln t$ | $\Lambda_{K}\,t$
square commutator | oscillations | $t^{2}$ | $e^{2\lambda t}$
Table 1:  Summary of the dynamical behaviour of entanglement and chaos
quantifiers of $N$-particle collective systems in the semiclassical regime.
The growth of the entanglement quantifiers and the square commutator depends
on the nature of the limiting classical trajectory in the $2n$-dimensional
phase space (stable configuration, regular or chaotic), up to the Ehrenfest
time. Here, $\lambda\equiv\lambda_{1}$ is the maximum finite time Lyapunov
exponent, and $\Lambda_{K}=\sum_{k=1}^{2K}\lambda_{k}$ is the sum of the $2K$
largest Lyapunov exponents, where $K$ is the number of degrees of freedom
associated with the considered subsystem. For $K$=$n/2$, one has the classical
Kolmogorov-Sinai entropy rate
$\Lambda_{\text{KS}}=\sum_{k\,:\,\lambda_{k}>0}\lambda_{k}$.

The ††margin:  Ehrenfest time formalism outlined in this Section is
quantitatively accurate as long as the number of collective excitations does
not grow too large $\langle\hat{n}_{0}\rangle\ll N$ compared to the system
size. As shown in Sec. 2.3 this assumption is generically valid for ground
states, even at the quantum critical points [84, 299]. Out of equilibrium,
this condition defines the Ehrenfest time scale, given by

$\langle\hat{n}_{0}(T_{\text{Ehr}})\rangle\sim N\,.$ (137)

On this time scale the quadratic truncation of the bosonic representation
loses accuracy. The non-linear corrections generally lead to saturation of the
growth of quantum fluctuations and to revivals on much longer times. Putting
everything together, we have

$\begin{dcases}\text{regular trajectories}&\langle\hat{n}_{0}(t)\rangle\sim
t^{2}\qquad T_{\text{Ehr}}\sim\hbar_{\text{eff}}^{-1/2}\sim\sqrt{N}\\\
\text{unstable (chaotic) trajectories}&\langle\hat{n}_{0}(t)\rangle\sim
e^{2\lambda t}\quad\,T_{\text{Ehr}}\sim\ln\hbar_{\text{eff}}^{-1/2}\sim\ln
N\end{dcases}\ .$ (138)

The dynamical growth of collective quantum fluctuations goes hand in hand with
the scrambling of quantum information and the dynamics of quantum
entanglement. As we will discuss in the next two sections, the approach
described here allows us to derive an exact relation between
$\langle\hat{n}_{0}(t)\rangle$, scrambling and entanglement.

Summary: Collective quantum fluctuations evolve as the linearized flow of
displacements around the classical trajectory before the Ehrenfest time scale.
The latter is defined as the time for which the number of quantum fluctuations
becomes comparable with the system size and consequently, it depends on the
classical phase-space.

####  4.1.4 Scrambling dynamics

Scrambling has been recently proposed as a pathway to characterize chaos in
many-body dynamics. Generically identified as delocalization of quantum
information, scrambling is commonly quantified by the dynamics of the _square-
commutator_ :

$c(t)=\langle\left|[\hat{A}(t),\hat{B}]\right|^{2}\rangle\ ,$ (139)

or of the closely related out-of-time order correlators (OTOC)
$\langle\hat{A}(t)\hat{B}\hat{A}(t)\hat{B}\rangle$, where the expectation
value is taken in a quantum state $\hat{\rho}$, i.e.,
$\langle\cdot\rangle=\text{Tr}(\cdot\hat{\rho})$. The ††margin:  Square-
commutator exponential growth square commutator was originally introduced by
Larkin and Ovchinnikov [300] to semi-classically describe the exponential
sensitivity to initial conditions.161616The heuristics goes as follows: for
$\hat{A}=\hat{x}$ and $\hat{B}=\hat{p}$ in the limit $\hbar\to 0$, upon
canonical quantization one has
$c(t)\simeq\hbar^{2}\\{x(t),p(0)\\}_{\text{PB}}^{2}\simeq\hbar^{2}\left|\frac{\partial
x(t)}{\partial x(0)}\right|^{2}\ .$ (140) Hence, $c(t)$ encodes the square of
the derivatives of the classical trajectory to respect to the initial
conditions. Thus, whenever the underlying classical limit is chaotic, $c(t)$
is expected to grow exponentially in time as

$c(t)\simeq\hbar_{\text{eff}}^{2}\,\,e^{\tilde{\lambda}t}\ ,$ (141)

with a rate $\tilde{\lambda}$ which may be related to the classical Lyapunov
exponent (but it is in principle distinct). This holds at intermediate times
before the Ehrenfest scale $t<T_{\text{Ehr}}\sim\ln\hbar_{\text{eff}}^{-1}$,
in this context also referred to as _scrambling time_. Interest in the square-
commutator was revived after Kitaev’s proposal to use it to characterize many-
body dynamics [301]. In this context, it was shown that the rate
$\tilde{\lambda}$ is upper bounded by quantum effects as
$\tilde{\lambda}\leq\frac{2\pi T}{\hbar}$ [302], as a consequence of the
quantum fluctuation-dissipation theorem [303, 304]. This constraint – now
known as “the bound to chaos” – is saturated by models of black holes,
including the Sachdev-Ye-Kitaev model (SYK) [305, 301], a system of fully
interacting disordered Majorana fermions where $\hbar_{\text{eff}}\sim\ln N$.

In the present case of fully-connected systems with a classical limit (Section
4.1.1), scrambling before the Ehrenfest time thus directly probes the
sensitivity of the classical trajectories to infinitesimal perturbations. One
can study the square commutator in Eq.(139) by taking the expectation value in
pure quasiclassical initial states and by looking at the square commutator
between two collective spin projections, namely

$c_{\alpha\beta}(t)=-\bigg{(}\frac{1}{Ns}\bigg{)}^{2}\bra{\psi_{0}}\left[\hat{S}^{\alpha}(t),\hat{S}^{\beta}(0)\right]^{2}\ket{\psi_{0}}\
,$ (142)

where $\alpha,\beta=x,y,z$ and $\ket{\psi_{0}}$ is a fully polarized spin-
coherent initial state. Using the expansion of the quantum fluctuations
elaborated in Section 4.1.3, we can compute the semiclassical evolution of the
out-of-time-order square commutator. By plugging the expansion of the rotated
spin operators (19) into the definition (142), one then substitutes the formal
solution for the spin fluctuations at time $t$, i.e.,
$\hat{Q}(t)=U_{qq}(t)\,\hat{Q}(0)+U_{qp}(t)\,\hat{P}(0)$ and
${\hat{P}(t)=U_{pq}(t)\,\hat{Q}(0)+U_{pp}(t)\,\hat{P}(0)}$. The initial
fluctuations for coherent states are
$\langle\hat{Q}(0)^{2}\rangle=\langle\hat{P}^{2}(0)\rangle=1/2$ and
$\langle\hat{Q}(0)\hat{P}(0)\rangle=\langle\hat{P}(0)\hat{Q}(0)\rangle=0$. The
resulting out-of-time square commutator in Eq.(142) thus reads ††margin:
Semiclassical result

$\displaystyle\begin{split}c_{\alpha\beta}(t)&=\Big{[}{X}_{\alpha}(t)\big{(}\,U_{qq}(t)\,\,{Y}_{\beta}(0)-U_{qp}(t)\,\,{X}_{\beta}(0)\,\big{)}+{Y}_{\alpha}(t)\,\big{(}\,U_{pq}(t)\,Y_{\beta}(0)-U_{pp}(t)\,X_{\beta}(0)\,\big{)}\Big{]}^{2}\\\
&\quad\quad+\mathcal{O}(\hbar_{\text{eff}})\ .\end{split}$ (143)

This expresses a quantitative relation between the square-commutator and the
formal evolution $U(t)$ [cf, below Eq. (133)] of the quantum fluctuations,
which encodes the of the evolution of linearised displacements and the _finite
time_ Lyapunov exponent spectrum $\\{\lambda_{k}(t)\\}$, as described in
Section 4.1.3. Hence, when the classical limit is integrable, the square-
commutator will grow as $c(t)\simeq t^{2}$, within the Ehrenfest time scale
$T_{\text{Ehr}}$. On the other hand in the presence of exponential sensitivity
associated with a phase-space separatrix (cf. Sec. 4.1.1 above) or with chaos,
the square-commutator $c(t)\simeq e^{2\lambda_{1}t}$ grows exponentially
before $T_{\text{Ehr}}$ with $\lambda_{1}=\lambda_{1}(t)$ the maximal finite-
time Lyapunov exponent of the underlying semiclassical trajectory. The
different scenarios are summarized in Table 1.

Results ††margin:  Comparison with numerics for the fully-connected quantum
Ising model are shown in Fig. 13, where we consider $c_{zz}(t)$ and compare
analytical result (black full line) with numerical exact diagonalization
results for finite system sizes. Parameters are specified in the caption. The
plot highlights the relation between entanglement entropy (discussed below)
and scrambling before the Ehrenfest time. As the considered model has an
integrable classical limit, the square-commutator only grows exponentially
$c_{zz}(t)\sim e^{2\lambda t}$ for quenches at the dynamical critical point
$h_{\rm dcr}$, associated with a classical separatrix with instability rate
$\lambda$ given in Eq. (127).

The exponential ††margin:  Chaotic fully-connected models sensitivity of the
square-commutator in the presence of the underlying separatrix has been
explored in a number of fully connected mean-field models [292, 306, 307, 308]
and recently probed experimentally [309]. When the initial state in Eq. (139)
is a random permutationally invariant state, the growth rate of the square-
commutator corresponds to the average of the finite-time Lyapunov over the
whole phase-space. In the presence of an instability, this leads to a modified
exponential growth $c{(t)}\sim e^{\lambda t}$ (rather than $e^{2\lambda t}$)
[310]. As we discussed above, fully-connected spin systems may exhibit
classically chaotic evolution when driven periodically (e.g. quantum kicked
top), coupled to other degrees of freedom (e.g. Dicke models), or for larger
individual spins $s>1/2$. Analysis as above predicts exponential growth of the
square-commutator for underlying classical chaos, as reported in the
literature for the quantum kicked top [311, 292, 312, 313, 307, 297, 313, 314,
315, 316, 317], the Dicke model [318, 319, 320, 321] and other spin models
[314, 322]. Recently, it has been pointed out that scrambling may become
super-exponential in fully connected models when the average in Eq.(139) is
done over at infinite temperature state [323].

Figure 13: Entanglement and scrambling dynamics after a quench of the
transverse field in the fully-connected Ising model, from $h_{0}=0$ to
$h_{f}>0$. The analytical prediction (black lines) is compared with exact
numerical results (colours) at finite $N=50,\,200,\,800$. We study the growth
of the entanglement entropy (145)(top panel) and the square commutator (142)
(bottom panel) quenching above, below, and at the dynamical phase transition
(DPT) at $h_{f}=h_{\rm dcr}$, as pictorially shown in Figure 12 (a). (a1-b1)
Quench above the DPT: $h_{f}=2J_{0}>h_{\rm dcr}$. (a2-b2) Quench below the
DPT: $h_{f}=0.2J_{0}<h_{\rm dcr}$. (a3-b3) Quench at the DPT:
$h_{f}=J_{0}/2=h_{\rm dcr}$. Time is measured in units of $J_{0}$. Plots
adapted from Ref.[324].

Summary: In fully-connected models, the square-commutator growth is
quantitatively determined by the semiclassical phase space before the
Ehrenfest time. Hence it encodes finite times Lyapunov exponents.

A similar statement also applies to other quantifiers of quantum information
spreading, in particular to entanglement entropy, which we turn to analyze in
the next Subsection.

####  4.1.5 Entanglement dynamics

It is by now well established that a large body of information about many-body
dynamics, their thermalization properties, and the complexity of their
numerical simulations, can be inferred from the evolution of bipartite
entanglement entropies. For ††margin:  Definitions a composite system with
Hilbert space ${\mathcal{H}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}}$ in a pure
state $\hat{\rho}=\ket{\psi}\bra{\psi}$, the bipartite entanglement between
subsystems $A$ and $B$ is encoded in the reduced density matrix
${\hat{\rho}_{A}=\text{Tr}_{B}\,{\hat{\rho}}}$.171717The nonvanishing
eigenvalues of ${\hat{\rho}_{B}=\text{Tr}_{A}\,{\hat{\rho}}}$ are equal to
those of $\hat{\rho}_{A}$. The system is entangled with respect to the
bipartition $(A,B)$ if $\hat{\rho}_{A}$ (equivalently, $\hat{\rho}_{B}$) is
not pure. The amount of bipartite entanglement can be quantified by the Renyi
entropies

$S_{A}^{\alpha}=-\frac{1}{1-\alpha}\ln\text{Tr}\,\hat{\rho}_{A}^{\alpha}\ ,$
(144)

parameterized by $\alpha>1$. The von Neumann entropy is obtained as their
limit for $\alpha\to 1$, i.e.,

${S_{A}=-\text{Tr}\big{(}\hat{\rho}_{A}\ln\hat{\rho}_{A}\big{)}}\ .$ (145)

As far ††margin:  Pictures for local models as local Hamiltonians are
concerned, the baseline features of entanglement entropy growth of pure states
out-of-equilibrium are well understood. As a broad consequence of the light-
cone spreading of quantum correlations (see discussion in Sec.3.1), in
thermalizing local systems the entanglement entropy $S_{A}(t)$ grows linearly
in time before saturating to a value proportional to its volume [112, 113,
325]. The underlying mechanism is well understood in integrable systems, where
it has been explained via a semiclassical picture based on quasi-particle
pairs propagation [326]. Analytical insights in chaotic many-body systems came
from the study of random unitary circuits [114]. On the other hand, the
presence of localized integrals of motion (exact or approximate) causes a
slowdown of entanglement growth, with a distinguished logarithmic increase for
many-body localized systems [327]. Neither of these scenarios is adequate for
long-range interacting systems, where numerical results exhibited slow
logarithmic entanglement growth even in the absence of quenched disorder [125,
127, 292].

A successful picture to capture entanglement growth in fully-connected systems
was only achieved more recently [58, 297]. The semiclassical growth of quantum
fluctuations (described in the section above) allows us to analytically relate
the dynamics of the entanglement entropy $S_{A}(t)$ to the quantifiers of
chaos, leading to a general unifying picture. This formalism yields a clean
prediction of logarithmic growth in the absence of semiclassical chaos. This
constitutes the origin of the slowdown of entanglement growth, which was first
observed in numerical simulations of dynamics in long-range interacting
quantum spin chains [127, 125, 292]: We will complete this discussion in
Section 4.2.5 below. We will illustrate how _$S_{A}(t)$ asymptotically
coincides with the logarithm of the phase space volume spanned by the quantum
fluctuations of the subsystem degrees of freedom_, as originally identified in
a seminal work by Zurek and Paz [328]; see Refs. [329, 330, 331, 58, 297] for
more recent literature.

Figure 14:  Entanglement dynamics in infinite range spin-chains. (a) The
system is partitioned into two blocks of $N_{A}$ and $N_{B}$ spins$1/2$,
initially fully polarized. (b) Collective spins of the two blocks. (c)
Collective spin in the factorized initial state, represented on the Bloch
sphere. The shaded area represents the quantum uncertainty of transverse
components. (d) Nonlinear interactions determine spin squeezing, which makes
the two blocks increasingly correlated (entangled). The rate of squeezing is
governed by the separation of nearby semiclassical trajectories, and by Eq.
(2) it determines the rate of growth of entanglement entropy. Right panels:
(e) For generic (noncritical) quenches, nearby trajectories separate linearly
in time, leading to a polynomially fast squeezing. (f) For a critical quench,
the collective spin lies on the stable manifold of an unstable fixed point in
phase space. In this case, nearby trajectories separate exponentially fast in
time at a rate $\lambda$ set by the eigenvalue of the linearized flow.

In ††margin:  Entanglement of all-to-all spin models the case of fully-
connected $N$-particle systems considered here, one considers a bipartition
between any two sets of spins, where the only relevant parameter is the number
$N_{A}=f_{A}N$ of particles in subsystem $A$ (with
$N_{B}=N-N_{A}=f_{B}N$).181818Due to permutational symmetry, spatial
bipartitions have no meaning. The collective spin $\hat{\vec{S}}$ can be
correspondingly decomposed as
$\hat{\vec{S}}=\hat{\vec{S}}_{A}+\hat{\vec{S}}_{B}$ (see Fig. 14). Within the
semiclassical description, the bipartite system can be represented by bosonic
operators $(\hat{q}_{A},\hat{q}_{B})$ and $(\hat{q}_{B},\hat{p}_{B})$,
associated with the quantum fluctuations of the two spins $\hat{\vec{S}}_{A}$
and $\hat{\vec{S}}_{B}$, respectively, via Holstein-Primakoff mapping. These
quantum fluctuations are characterized by the correlation matrix $G(t)$
defined in Eq. (134). It is convenient to define the subsystem’s _reduced
correlation matrix_ $G_{A}(t)$ as the $2\times 2$ matrix of quantum
fluctuations built out of the variables of subsystem $A$ alone, i.e.,

$\begin{split}G_{A}&=\begin{pmatrix}\langle\hat{q}_{A}^{2}\rangle&\frac{\langle\hat{q}_{A}\hat{p}_{A}+\hat{p}_{A}\hat{q}_{A}\rangle}{2}\\\
\frac{\langle\hat{q}_{A}\hat{p}_{A}+\hat{p}_{A}\hat{q}_{A}\rangle}{2}&\langle\hat{p}_{A}^{2}\rangle\end{pmatrix}\equiv\begin{pmatrix}G^{q_{A}q_{A}}&G^{q_{A}p_{A}}\\\
G^{q_{A}p_{A}}&G^{p_{A}p_{A}}\end{pmatrix}.\end{split}$ (146)

In the semiclassical regime of small $\hbar_{\text{eff}}$, the reduced density
matrix $\hat{\rho}_{A}(t)$ is asymptotically Gaussian to leading order, and
thus fully determined by $G_{A}(t)$. The entanglement properties can thus be
computed via standard techniques [329], see also Refs. [330, 331]. The Von
Neumann and the second Renyi entropies of a single boson
$(\hat{q}_{A},\hat{p}_{A})$ in such a Gaussian state can be expressed in terms
of the determinant of $G_{A}$ as [332] ††margin:  Entanglement from quadratic
fluctuations

$\displaystyle S_{A}=2\sqrt{\det G_{A}}\text{arccoth}$
$\displaystyle\left(2\sqrt{\det G_{A}}\right)+\frac{1}{2}\log\left(\det
G_{A}-\frac{1}{4}\right)\ ,$ (147a) $\displaystyle S^{(2)}_{A}(t)$
$\displaystyle=\frac{1}{2}\ln\;\det\big{(}2G_{A}(t)\big{)}\ .$ (147b)

On the other hand, the matrix $G_{A}$ can be directly related to the
correlation matrix $G$ of collective excitations $(\hat{q},\hat{p})$.191919One
can perform a linear canonical transformation to the collective
$(\hat{q},\hat{p})$ and relative $(\delta\hat{q},\delta\hat{p})$ fluctuation
modes:
$\begin{dcases}\hat{q}=+\sqrt{f_{A}}\;\hat{q}_{A}+\sqrt{f_{B}}\;\hat{q}_{B}\\\
\delta\hat{q}=-\sqrt{f_{B}}\;\hat{q}_{A}+\sqrt{f_{A}}\;\hat{q}_{B}\end{dcases}\begin{dcases}\hat{p}=+\sqrt{f_{A}}\;\hat{p}_{A}+\sqrt{f_{B}}\;\hat{p}_{B}\\\
\delta\hat{p}=-\sqrt{f_{B}}\;\hat{p}_{A}+\sqrt{f_{A}}\;\hat{p}_{B}\end{dcases}\
.$ (148) Since the Hamiltonian is a function of the collective spin only, the
latter bosonic mode is frozen in the vacuum.  The explicit computation shows
that the determinant can be expressed as

$\det G_{A}=\frac{1}{4}+f_{A}f_{B}\;\langle\hat{n}_{0}\rangle$ (149)

where $\hat{n}_{0}=(\hat{q}^{2}+\hat{p}^{2}-1)/2$ represents the number of
bosonic excitations of the collective spin [cf. Eq. (135)]. While the global
evolution preserves the total volume, i.e., $\det\big{(}2G(t)\big{)}\equiv 1$,
the information loss generated by projecting the collective quantum
fluctuations onto a subsystem yields an increase of entropy. By Eq. (147b),
this increase may be visualized as an enhancement of the projected volume
spanned by the reduced quantum fluctuations within the subsystem phase space,
due to the progressive stretching of the phase-space volume spanned by the
quantum fluctuations, see Fig. 14. The ††margin:  Entanglement dynamics and
classical phase-space growth of entanglement entropy out of equilibrium is
thus completely determined by the dynamical generation of the collective
excitations
$\langle\hat{n}_{0}(t)\rangle=\frac{1}{2}\text{Tr}\left[G(t)-\frac{\mathbb{1}}{2}\right]$.
As discussed in Section 4.1.3, $G(t)$ describes the flow of linearized
displacements around the classical trajectories. This connection highlights
that the entanglement growth in the semiclassical regime is determined by the
chaoticity properties of the underlying classical phase space. In fact, the
qualitative time-dependence of $\langle\hat{n}_{0}(t)\rangle$ depends on the
nature of the classical trajectories. See Table 1 for a summary.

The classical dynamics of fully-connected spin-$1/2$ systems are generically
integrable (as discussed in Section 4.1.1). Hence the temporal growth of the
quantum correlations is at most polynomial, $\langle n_{exc}(t)\rangle\sim
t^{2}$ [see Fig. 5(e)], leading to

$S_{A}(t)\underset{t\gg 1}{\thicksim}S_{A}^{(2)}(t)\sim\ln t\ .$ (150)

In the (non-generic) case of quenches to dynamical critical points (see
discussion in Section 4.1.2), the collective spin moves along an isolated
unstable trajectory called separatrix [see Fig. 5(f)]. Out-of-equilibrium
generation of collective excitations is thus exponentially fast in such
critical quenches, $\langle\hat{n}_{\text{exc}}(t)\rangle\sim e^{2\lambda t}$,
leading to a linear growth of entanglement entropy with a predicted slope

$S_{A}(t)\underset{t\gg 1}{\thicksim}S_{A}^{(2)}(t)\sim\lambda t\ .$ (151)

This phenomenology fully describes the entanglement entropy dynamics of the
fully-connected quantum Ising model. This is shown in Fig. 13, where
$S_{A}(t)$ is studied for quenches in the transverse field from $h_{0}=0$ to
different $h_{f}$ below above and at the dynamical critical point $h_{\rm
dcr}$, discussed in Sec. 4.1.2.

This analysis can be extended to general fully-connected systems whose
classical limit has $n>1$ degrees of freedom [297] and may exhibit chaos. In
this case, quantum fluctuations grow as
$\langle\hat{n}_{\text{exc}}(t)\rangle\sim e^{2\lambda t}$, and the growth of
the entanglement entropy $S_{A}(t)$ is generically linear in time with a rate
set by the sum of the largest $2n_{A}$ Lyapunov exponents [333, 330, 297]:

$S_{A}(t)\underset{t\gg
1}{\thicksim}S^{(2)}_{A}(t)\thicksim\Lambda_{A}t=\bigg{(}\sum_{k=1}^{2n_{A}}\lambda_{k}\bigg{)}t\
.$ (152)

For $n_{A}=n/2$, this rate coincides with the classical Kolmogorov-Sinai
entropy rate $\Lambda_{KS}=\sum_{\lambda_{k}\,:\,\lambda_{k}>0}\lambda_{k}$
[298]. A linear growth of entanglement entropy thus occurs in chaotic fully-
connected spin systems, such as the quantum kicked top [334, 335, 336, 337,
338, 339, 340, 341, 342, 343, 344, 345, 346, 292], the Dicke model [347, 339,
348, 349, 350, 351, 352, 353, 354, 355, 319, 297], and larger-$s$ fully-
connected systems.

The classification above is concluded by the case of near-equilibrium dynamics
around stable equilibrium configurations. In this case, the linearized
dynamics are equivalent to that of coupled harmonic oscillators, leading to
persistent oscillations in entanglement dynamics. See Table 1 for a summary.

It ††margin:  Squeezing induced picture is worth noting that the stretching of
collective quantum fluctuations in phase space, which lies at the origin of
bipartite entanglement growth, is very explicitly connected with important
witnesses of multipartite quantum entanglement, such as _spin squeezing_ [356]
and _the Quantum Fisher Information_ [357]. A popular quantifier of spin-
squeezing is defined by the minimal transverse variance of collective quantum
spin fluctuations [358, 359]:

$\xi^{2}\equiv\min_{\lvert\mathbf{u}\rvert=1,\mathbf{u}\perp\mathbf{Z}}\frac{\Big{\langle}\big{(}\mathbf{u}\cdot\mathbf{\hat{S}}\big{)}^{2}\Big{\rangle}}{N/4}.$
(153)

The squeezing parameter $\xi^{2}$ is equal to one $1$ for coherent (separable)
states, while it is smaller for squeezed states $\xi^{2}<1$. It has long been
known [360, 361] that collective spin squeezing is a witness of many-body
quantum entanglement, which can be generated by fully connected interactions
[358, 359]. Indeed, over the timescale $t\ll T_{\rm Ehr}$, $\xi(t)$, squeezing
is explicitly related to $\langle n_{0}(t)\rangle$ in Eq. (149) as

$\xi^{2}(t)=1+2\langle\hat{n}_{\text{0}}(t)\rangle-2\sqrt{\langle\hat{n}_{\text{0}}(t)\rangle(\langle\hat{n}_{\text{0}}(t)\rangle+1)}\;.$
(154)

Equations (147a, 149) and (154) express the quantitative link — pictorially
illustrated in Figure 14 — between the entanglement entropy $S_{A}$ and and
the spin squeezing parameter $\xi$, in collective spin models in the
semiclassical regime in and out of equilibrium. Following this relation, we
will refer to the mechanism of dynamical entanglement entropy growth outlined
here as _spin-squeezing picture_ [58]. For definiteness in this Report, we
focus on bipartite entanglement entropies; we refer the readers to Ref. [297]
for further details on the dynamics of multipartite entanglement.

In conclusion, we recall that the analysis presented above is valid before the
Ehrenfest time scale defined in Eq. (137). Over longer times entanglement
entropies saturate, $S^{\infty}_{A}\propto\log N_{A}$.

Summary: In fully connected systems, entanglement dynamics follow a
semiclassical picture (before Ehrenfest time), which relates it to the
squeezing of quantum fluctuations in phase space. This relation generically
predicts a logarithmic entanglement growth, in the absence of semiclassical
chaos.

###  4.2 Quench dynamics of long-range interacting spin systems ($\alpha>0$)

As soon as the interaction range is finite, the full permutational symmetry of
the problem is broken and, in principle, the system could thermalize. The
purpose of this Section is to describe the persistent collective character of
dynamics in systems with $0<\alpha<d$. By formulating a non-equilibrium spin-
wave theory [56, 57], (Sec. 4.2.1), we will be able to develop a physical
picture in terms of a semiclassical collective degree of freedom coupled to
excitations with finite wavelengths. Analysis of the resulting non-linear
coupled equations allows us to demonstrate freezing of finite-wavelength modes
for long times, resulting in lower bounds to thermalization time scales (Sec.
4.2.2). This scenario affects DPTs (Sec. 4.2.3), scrambling dynamics (Sec.
4.2.4), and the characteristic slow growth of entanglement entropy (Sec.
4.2.5).

####  4.2.1 Dynamics of quantum fluctuations with finite interaction range

To study the impact of finite range interactions on the mean-field dynamics,
we resort to a _non-equilibrium spin-wave theory_ developed in Refs. [56, 57].
The goal is to refine the time-dependent formalism of Sec. 4.1.3 to the full
momentum-space representation of the finite-range spin Hamiltonian, Eqs. (28)
and (29). Similarly to our discussion of equilibrium properties in Sec. 2.4,
one can single out a collective $\alpha$-independent part of the Hamiltonian
involving $\mathbf{k}=\mathbf{0}$ spin operators only, and an
$\alpha$-dependent perturbation which activates spin fluctuation modes at
$\mathbf{k}\neq\mathbf{0}$.

We aim at expanding the individual spins as bosonic fluctuations around a yet
unspecified time-dependent quantization axis ${\mathbf{Z}}$ — which we will
later self-consistently require to coincide with the instantaneous direction
of the collective spin $\langle\vec{\hat{S}}(t)\rangle\propto{\mathbf{Z}}(t)$.
To achieve this, one performs the time-dependent rotation generated by
$\hat{V}(t)$ in Eq. (129). The spin components in this time-dependent frame
are governed by the inertial Hamiltonian
$\tilde{\hat{H}}(t)=\hat{V}\,\hat{H}\,\hat{V}^{\dagger}+i\hat{V}\dot{\hat{V}}^{\dagger}\
,$ as in Sec. 4.1.3. One then applies Holstein-Primakoff transformations to
the individual rotating spins, as in Eq. (60). The resulting transformed
Hamiltonian can be organized in the usual form ††margin:  Time-dependent
Holstein Primakoff

$\tilde{H}(t)=\frac{1}{s}\bigg{[}(Ns)^{1}\tilde{\mathcal{E}}_{0}(t)+(Ns)^{1/2}\tilde{H}_{1}(t)+(Ns)^{0}\tilde{H}_{2}(t)+(Ns)^{-1/2}\tilde{H}_{3}(t)+(Ns)^{-1}\tilde{H}_{4}(t)+\dots\bigg{]}\
.$ (155)

As we already observed in the discussion of equilibrium properties, comparison
to the expansion for $\alpha=0$ shows that
$\hat{q}\equiv\tilde{q}_{\mathbf{k}=\mathbf{0}}\
,\,\hat{p}\equiv\tilde{p}_{\mathbf{k}=\mathbf{0}}\
,\,\hat{n}_{0}\equiv\hat{n}_{\mathbf{k}=\mathbf{0}}$, and that the total
occupation number of spin-wave excitations $\hat{n}_{\text{sw}}$ is given by
the sum of bosonic occupation numbers of all the other spin-wave modes at
finite wavelength, cf. Eq. (39)

$\hat{n}_{\text{sw}}=\sum_{\mathbf{k}\neq\mathbf{0}}\hat{n}_{\mathbf{k}}\
,\quad\text{with}\quad\hat{n}_{\mathbf{k}}\equiv\frac{\tilde{q}_{\mathbf{k}}\tilde{q}_{-\mathbf{k}}+\tilde{p}_{\mathbf{k}}\tilde{p}_{-\mathbf{k}}-1}{2}\
.$ (156)

The individual occupation numbers quantities
$\hat{n}_{\mathbf{k}\neq\mathbf{0}}$ are exactly conserved by the collective
part of the Hamiltonian $\tilde{H}_{\alpha=0}(t)$, which only depends on
$\mathbf{k}\neq\mathbf{0}$ bosons through the collective spin length [i.e.
through $\hat{n}_{\text{sw}}$].

The equations of motion of the classical angles $\theta(t)$, $\phi(t)$ is once
again found by imposing the condition that
$\braket{\hat{S}^{X}(t)}=\braket{\hat{S}^{Y}(t)}=0$ [56, 57], namely that the
collective bosonic mode describes fluctuations around the instantaneous
average spin polarization. This amount to setting the coefficients of
$\tilde{q}_{\mathbf{k}=\mathbf{0}}$ and $\tilde{p}_{\mathbf{k}=\mathbf{0}}$
equal to zero. Taking into account the leading term $\tilde{H}_{1}(t)$ only,
one retrieves the usual classical mean-field equations of motion (126).

However, as first demonstrated in Refs. [56, 57], in the presence of finite-
range interactions the collective spin trajectory may get modified by quantum
fluctuations. This effect is the non-equilibrium counterpart of the
corrections to the equilibrium spin polarization arising from a finite spin-
wave density in the quantum ferromagnetic phase, cf. Sec. 2.4.3. The
corrections arise from the terms in $\tilde{H}_{3}(t)$ involving a bosonic
operator with $\mathbf{k}=\mathbf{0}$ (such that the remaining two operators
have momenta $\pm\mathbf{k}$). In physical terms, these interactions describe
the scattering of a collective spin excitation into a pair of spin waves with
opposite finite momenta and vice versa. Taking into account this “feedback”
from quantum fluctuations, one obtains a pair of modified equations for the
angles $\theta(t),\phi(t)$: ††margin:  Equation of motion with the feedback

$\displaystyle\frac{d}{dt}\theta=$
$\displaystyle+2J_{0}(1-\epsilon)\sin\theta\cos\phi\sin\phi$
$\displaystyle-2J_{0}\bigg{(}\frac{1}{Ns}\sum_{\mathbf{k}}f_{\mathbf{k}}(\alpha)\;\left\langle\tilde{p}_{\mathbf{k}}\tilde{p}_{-\mathbf{k}}\right\rangle\bigg{)}\sin\theta\cos\phi\sin\phi$
(157a)
$\displaystyle+2J_{0}\Bigg{(}\frac{1}{Ns}\sum_{\mathbf{k}}f_{\mathbf{k}}(\alpha)\frac{\left\langle\tilde{q}_{\mathbf{k}}\tilde{p}_{-\mathbf{k}}+\tilde{p}_{\mathbf{k}}\tilde{q}_{-\mathbf{k}}\right\rangle}{2}\Bigg{)}\cos\theta\sin\theta\cos^{2}\phi,$
$\displaystyle\frac{d}{dt}\phi=$
$\displaystyle-h+2J_{0}(1-\epsilon)\cos\theta\cos^{2}\phi$
$\displaystyle-2J_{0}\bigg{(}\frac{1}{Ns}\sum_{\mathbf{k}\neq
0}f_{\mathbf{k}}(\alpha)\;\left\langle\tilde{q}_{\mathbf{k}}\tilde{q}_{-\mathbf{k}}\right\rangle\bigg{)}\cos\theta\cos^{2}\phi$
(157b) $\displaystyle+2J_{0}\Bigg{(}\frac{1}{Ns}\sum_{\mathbf{k}\neq
0}f_{\mathbf{k}}(\alpha)\frac{\left\langle\tilde{q}_{\mathbf{k}}\tilde{p}_{-\mathbf{k}}+\tilde{p}_{\mathbf{k}}\tilde{q}_{-\mathbf{k}}\right\rangle}{2}\Bigg{)}\sin\phi\cos\phi,$

where we introduced the time-dependent _spin-wave density_

$\epsilon(t)\equiv\frac{\langle\hat{n}_{\text{tot}}(t)\rangle}{Ns}=\frac{1}{Ns}\sum_{\mathbf{k}}\frac{\langle\tilde{q}_{\mathbf{k}}\tilde{q}_{-\mathbf{k}}\rangle+\langle\tilde{p}_{\mathbf{k}}\tilde{p}_{-\mathbf{k}}\rangle-1}{2}\,.$
(158)

As usual, we observe that the impact of quantum fluctuations on the classical
trajectory is suppressed in the classical limit $s\to\infty$, and grows with
$\alpha$ at fixed $s$. Similarly to what we found in equilibrium (Sec. 2.4),
the properties of $f_{\mathbf{k}}(\alpha)$ imply that all the quantum feedback
terms in the right-hand sides of the above equations of motion are vanishingly
small in the thermodynamic limit for $0<\alpha<d$ and parametrically small for
$\alpha\searrow d$. In the latter case, we can make the usual replacement
$(1/N)\sum_{\mathbf{k}}\mapsto\int d^{d}\mathbf{k}/(2\pi)^{d}$ in the
thermodynamic limit. However, in finite systems those corrections could be
expected to play a role for arbitrary $\alpha$ at sufficiently long times.

In turn, ††margin:  Spin-waves Hamiltonian the evolution of quantum
fluctuations is regulated by the spin-wave Hamiltonian $\tilde{H}_{2}(t)$ to
the same order of approximation as above. It is instructive to report its
explicit expression for the considered variable-range quantum Ising model:

$\hat{H}_{2}(t)=-2J_{0}\sum_{\mathbf{k}}f_{\mathbf{k}}(\alpha)\,\times\\\
\left(\cos^{2}\theta\cos^{2}\phi\,\frac{\tilde{q}_{\mathbf{k}}\tilde{q}_{-\mathbf{k}}}{2}+\sin^{2}\phi\;\frac{\tilde{p}_{\mathbf{k}}\tilde{p}_{-\mathbf{k}}}{2}-\cos\theta\sin\phi\cos\phi\,\frac{\tilde{p}_{\mathbf{k}}\tilde{q}_{-\mathbf{k}}+\tilde{q}_{\mathbf{k}}\tilde{p}_{-\mathbf{k}}}{2}\right)\\\
+2J_{0}\cos^{2}\phi\,\sum_{\mathbf{k}}\hat{n}_{\mathbf{k}}\,.$ (159)

This Hamiltonian is equivalent to a set of _externally driven_ quantum
harmonic oscillators, labelled by the momentum $\mathbf{k}$, where the driving
is given by the motion of $\theta(t)$ and $\phi(t)$ and controlled by the
couplings $f_{\mathbf{k}}(\alpha)$.

To close ††margin:  Equations of motion for quantum fluctuations the system of
equations of motion it is convenient to define the momentum-resolved
correlation functions

$\displaystyle
G^{qq}_{\mathbf{k}}(t)\equiv\langle\,\tilde{q}_{\mathbf{k}}(t)\tilde{q}_{-\mathbf{k}}(t)\,\rangle,\quad$
$\displaystyle
G^{pp}_{\mathbf{k}}(t)\equiv\langle\,\tilde{p}_{\mathbf{k}}(t)\tilde{p}_{-\mathbf{k}}(t)\,\rangle,\quad$
(160a) $\displaystyle
G^{qp}_{\mathbf{k}}(t)\equiv\frac{1}{2}\langle\,\tilde{q}_{\mathbf{k}}(t)\tilde{p}_{-\mathbf{k}}(t)$
$\displaystyle+\,\tilde{p}_{\mathbf{k}}(t)\tilde{q}_{-\mathbf{k}}(t)\rangle\
.$ (160b)

Starting from the Heisenberg equations
$\frac{d}{dt}\tilde{q}_{\mathbf{k}}=i[\tilde{H}_{2}(t),\tilde{q}_{\mathbf{k}}]$
and
$\frac{d}{dt}\tilde{p}_{\mathbf{k}}=i[\tilde{H}_{2}(t),\tilde{p}_{\mathbf{k}}]$
we compute

$\begin{dcases}\dot{G}^{qq}_{\mathbf{k}}=\,4J_{0}f_{\mathbf{k}}(\alpha)\,\cos\theta\cos\phi\sin\phi\;G^{qq}_{\mathbf{k}}+4J_{0}\left(\cos^{2}\phi-
f_{\mathbf{k}}(\alpha)\,\sin^{2}\phi\right)\,G^{qp}_{\mathbf{k}}\,,\\\
\dot{G}^{pp}_{\mathbf{k}}=-4J_{0}\left(\cos^{2}\phi-
f_{\mathbf{k}}(\alpha)\,\,\cos^{2}\theta\cos^{2}\phi\right)G^{qp}_{\mathbf{k}}-4J_{0}f_{\mathbf{k}}(\alpha)\,\,\cos\theta\cos\phi\sin\phi\,G^{pp}_{\mathbf{k}}\,,\\\
\dot{G}^{qp}_{\mathbf{k}}=-2J_{0}\left(\cos^{2}\phi-
f_{\mathbf{k}}(\alpha)\,\,\cos^{2}\theta\cos^{2}\phi\right)G^{qq}_{\mathbf{k}}+2J_{0}\left(\cos^{2}\phi-
f_{\mathbf{k}}(\alpha)\,\,\sin^{2}\phi\right)G^{pp}_{\mathbf{k}}\,.\end{dcases}$
(161)

Like Eqs. (136), these equations are also not independent due to the relation
$4(G^{pq}_{\mathbf{k}})^{2}=4\;G^{pp}_{\mathbf{k}}\;G^{qq}_{\mathbf{k}}-1$,
which is an exact property of Gaussian pure states, and which is then
satisfied at all times and for all ${\mathbf{k}}$’s to the considered level of
approximation.

The _general physical picture_ is now clear:

  * •

To lowest order the collective spin follows the classical mean-field
trajectory;

  * •

This collective spin motion acts as an external drive for the spin-wave
excitations, whereby the couplings $f_{\mathbf{k}}(\alpha)$ control the
_driving amplitude_ ;

  * •

The dynamically populated non-equilibrium spin-wave “bath” may in turn back-
react to modify the collective spin dynamics via the quantum feedback terms.

To quadratic order of approximation in the spin waves, the quantum many-body
dynamics of the system is described by the closed set of coupled non-linear
evolution equations (157) and (161), together with suitable initial conditions
(which may be a ground or thermal state of a pre-quench Hamiltonian — see the
discussion on equilibrium states in Sec. 2.4). This effective decoupling
between the dominant zero mode and the finite $k$-suppressed spin waves has
recently been exploited in Refs.[362, 363]. There, the quadratic zero mode is
replaced by a (fully quantum) rotor, while the $k$ spin waves are kept at the
quadratic level. This approach allows to reproduce the dynamics of quantum
fluctuations beyond Ehrenfest time.

The formalism derived above is expected to provide an accurate description of
the time-evolving many-body wave function whenever the dynamically generated
spin-wave density $\epsilon(t)$ remains small, i.e. $\epsilon(t)\ll 1$. This
diluteness condition allows us to theoretically describe the many-body
dynamics as a self-consistent driven weakly-interacting bosonic gas. The
neglected higher-order terms account for the non-linear scattering among spin
waves, which is expected to contribute to the dynamics only over time scales
parametrically long in $1/\epsilon$. Physically, the diluteness condition
above corresponds to the requirement that the time-evolving collective spin
magnitude $S$ remains close to its maximal value $Ns$, as showcased by Eq.
(39). The quality of this approximation is significantly impacted by the
interaction range via $f_{\mathbf{k}}(\alpha)$.

Below we will use the formalism outlined above to review how the finite-range
of interactions impacts the mean-field DPT, scrambling, and entanglement
dynamics.

Summary: Finite-wavelength fluctuations can be modeled as a set of driven
bosonic modes, where the drive is given by the collective spin motion. The
quantum fluctuations, in turn, may back-react on the collective spin dynamics.
The full dynamics is thus described by a set of coupled non-linear evolution
equations.

####  4.2.2 Prethermal freezing of spin-wave excitations

The dynamical generation of spin-wave excitations with non-vanishing momenta
for $\alpha>0$ is responsible for modififcations to the mean-field dynamics.
As is manifest in Eqs. (157), the impact of the quantum feedback is controlled
by the finite-range perturbation via $f_{\mathbf{k}\neq\mathbf{0}}(\alpha)$.
In turn, the same quantities $f_{\mathbf{k}\neq\mathbf{0}}(\alpha)$ bound the
rate itself of dynamical generation of spin waves: From Eq. (159) it is clear
that $\big{[}\hat{n}_{\mathbf{k}\neq\mathbf{0}},\tilde{H}(t)\big{]}\propto
f_{\mathbf{k}\neq\mathbf{0}}(\alpha)$, consistently with the sum of the first
two equations in (161).

To formulate a more precise bound we can proceed as follows: First, we can
eliminate the free spin-wave precession, i.e. the last term in Eq. (159), by
switching to the “interaction picture”:
$\tilde{q}_{\mathbf{k}}\mapsto\cos\Phi(t)\tilde{q}_{\mathbf{k}}+\sin\Phi(t)\tilde{p}_{\mathbf{k}}$,
$\tilde{p}_{\mathbf{k}}\mapsto-\sin\Phi(t)\tilde{q}_{\mathbf{k}}+\cos\Phi(t)\tilde{p}_{\mathbf{k}}$,
where $\Phi(t)=\int_{0}^{t}ds\,2J_{0}\cos^{2}\phi(s)$.202020While the axis
$\mathbb{Z}$ is fixed by a self-consistency requirement, the orientation of
the transverse axes $\mathbb{X}$, $\mathbb{Y}$ is a “gauge freedom” within the
formalism. The transformation above amounts to performing a time-dependent
rotation of the transverse spin components by angle $\Phi(t)$. With this
choice of co-moving frame, the spin-wave Hamiltonian fully vanishes for
$\alpha=0$: spin fluctuations look frozen. Such a “gauge” transformation does
not change the dynamical population of spin waves. All terms of the right-hand
side of the modified Eq. (161) are now proportional to
$f_{\mathbf{k}\neq\mathbf{0}}(\alpha)$. We can then apply Gronwall’s lemma
[364] to this linear system of differential equations to bound the growth of
the $G_{\mathbf{k}}(t)$’s: In particular, ††margin:  Spin waves suppression

$\langle\hat{n}_{\mathbf{k}}(t)\rangle=\frac{1}{2}\left(G^{qq}_{\mathbf{k}}+G^{pp}_{\mathbf{k}}(t)-1\right)\leq\frac{1}{2}\left[\exp\Big{(}c|f_{\mathbf{k}}(\alpha)|J_{0}t\Big{)}-1\right]\,,$
(162)

where $c$ is a constant related to the norm of the monodromy matrix on the
right-hand side of Eq. (161). Thus, spin-wave excitations at momentum
$\mathbf{k}$ can only be dynamically generated over time scales $J_{0}t\gg
1/f_{\mathbf{k}}(\alpha)$.

Of course, the bound in Eq. (162) is only useful for $\alpha\lesssim d$. As
discussed in Sec. 2.4 and in C, the couplings
$f_{\mathbf{k}\neq\mathbf{0}}(\alpha)$ are suppressed in the thermodynamic
limit $N\to\infty$ when $0<\alpha<d$ and are parametrically small as
$\alpha\searrow d$; see Fig. 3 for an illustration. This straightforwardly
provides us with a bound on the rate of growth of the population of bosonic
excitations for $0<\alpha<d$ (see C):

$|f_{\mathbf{k}}(\alpha)|\leq\;\frac{\text{const}}{(\lvert\mathbf{k}\rvert
L)^{\beta}}\
,\quad\text{with}\quad\beta\equiv\min\big{(}d-\alpha,(d+1)/2\big{)}\,.$ (163)

Therefore, there exists a long time scale

$T_{\text{sw}}\sim N^{\beta/d}\ ,$ (164)

during which the dynamical excitation of spin waves with _finite wavelengths_
is suppressed.212121Note that this further a posteriori justifies the
Holstein-Primakoff approach, as the density of spin waves remains small over a
long time window. We can further easily bound the total density of spin-wave
excitations _at finite times_ :

$\epsilon(t)\leq\frac{1}{Ns}\sum_{\mathbf{k}}\frac{1}{2}\left[\exp\Big{(}c|f_{\mathbf{k}}(\alpha)|J_{0}t\Big{)}-1\right]\sim\frac{J_{0}t}{N}\sum_{\mathbf{k}}|f_{\mathbf{k}}(\alpha)|$
(165)

For $0<\alpha<d$ we have [recall
$\mathbf{k}\equiv\mathbf{k}_{\boldsymbol{\ell}}=2\pi\boldsymbol{\ell}/L$, with
$\boldsymbol{\ell}=(\ell_{1},\dots,\ell_{d})$ integers]

$\epsilon(t)\leq\frac{J_{0}t}{L^{d}}\,\sum_{|\boldsymbol{\ell}|<L/2}\frac{1}{|\boldsymbol{\ell}|^{\beta}}\sim\frac{J_{0}t}{L^{\beta}}$
(166)

This bound proves that the spin-wave formalism described here is
asymptotically exact in this regime and that the mean-field description of the
collective spin polarization dynamics becomes exact in the thermodynamic
limit. (We stress that this conclusion requires that $N\to\infty$ must be
taken before $t\to\infty$.)

It ††margin:  Stability analysis of the spin waves must be noted that the
bound (162) above is overly pessimistic, as it considers exponential
instability for _all_ the bosonic modes (and with the worst possible rates).
In reality, our system can be approximately viewed as a set of quantum
harmonic oscillators subject to periodic parametric drive at a frequency given
by the classical collective spin precession. Such driven oscillators may or
may not meet a resonance; for a given quench, resonances can be detected by
performing a stability analysis of the stroboscopic (Floquet) evolution
operator for each spin-wave mode [58]:

$\begin{pmatrix}\tilde{q}_{\mathbf{k}}(t_{0}+T_{\text{cl}})\\\
\tilde{p}_{\mathbf{k}}(t_{0}+T_{\text{cl}})\end{pmatrix}=U_{\mathbf{k}}(T_{\text{cl}})\cdot\begin{pmatrix}\tilde{q}_{\mathbf{k}}(t_{0})\\\
\tilde{p}_{\mathbf{k}}(t_{0})\end{pmatrix}$ (167)

over the period $T_{\text{cl}}$ of the motion of the angles
$\theta(t),\phi(t)$. The eigenvalues
$e^{\pm\lambda_{\mathbf{k}}T_{\text{cl}}}$ of the $2\times 2$ matrix
$U_{\mathbf{k}}(T_{\text{cl}})$ directly give information on the _Floquet
quasi-frequency_ $\lambda_{\mathbf{k}}$ (see, e.g., Refs. [81, 365]) of the
driven oscillator, which may be purely real (resonance at mode $\mathbf{k}$)
or purely imaginary $\lambda_{\mathbf{k}}=i\omega_{\mathbf{k}}$ (non-
resonance). In the latter case, the mode population performs bounded
quasiperiodic oscillations, whereas in the former case, it blows up
exponentially. The real parts of the Floquet quasi-frequencies play the role
of finite-time Lyapunov exponents for our system. To quantify the global
stability of the system’s dynamics for a given quantum quench, one can compute
the sum of all _positive_ Lyapunov exponents, which gives the Kolmogorov-Sinai
entropy rate $\Lambda_{\text{KS}}$. It is convenient to inspect this quantity
on varying the fully polarized initial configurations, parametrized by
spherical angles:

$\Lambda_{\text{KS}}(\theta_{0},\phi_{0})=\sum_{{\mathbf{k}}:\lambda_{\mathbf{k}}>0}\text{Re}[\lambda_{\mathbf{k}}(\theta_{0},\phi_{0})]\
.$ (168)

This quantity detects whether resonances are present for the considered
quench, and quantifies the actual (initial) instability. While the stability
analysis described here can be done for arbitrary $\alpha$, it is only really
meaningful for $\alpha\lesssim d$, as for larger $\alpha$ the density of spin-
waves becomes finite in a finite time, and non-linear effects in the full
many-body dynamics cannot be neglected. For $0<\alpha<d$ the spin waves
effectively reduce to a _discrete_ set of periodically driven quantum
oscillators associated with long-wavelength modes with ${\mathbb{k}}\propto
1/L$, cf. the discussion in Sec. 2.4.

The stability ††margin:  Spin waves stability in the Ising chain the analysis
described above was performed for a long-range quantum Ising chain with
$0<\alpha<1$ in Ref. [58]. For a large set of initial conditions, spin waves
are found to be _stable_ (i.e. non-resonant), and consequently, their
population remains bounded in time. While there is no simple rule to predict
the existence of resonances, numerical exploration suggests that quenches
_near_ dynamical criticality typically give rise to resonant excitation of
long-wavelength spin waves. In other words, the classical separatrix of the
mean-field dynamics for $\alpha=0$ typically broadens to a _finite_ layer of
instability (chaoticity) for $\alpha>0$. In the left panel of Fig. 15 we
report the $k$-resolved spin-wave dynamical population for $\alpha=0.7$,
obtained by numerically solving Eq. (161). The occurrence of resonances is
systematically illustrated in the right panel of Fig. 15, which displays the
value of $\Lambda_{\text{KS}}(\theta_{0},\phi_{0})$ as a function of the
initial configuration. These results show that, at least in this case, only
quenches near dynamical criticality give rise to resonant excitation of spin
waves.

For quenches that do exhibit resonances at certain non-vanishing momenta, the
bound (162) captures the qualitative time-dependence of the population of the
corresponding modes (the actual rate $\lambda_{\mathbf{k}}$ is, however,
generally lower). Thus, the rapidly growing population of excitations is
expected to generate non-linear effects in the full many-body dynamics over
comparatively short time scales of order

$T_{\rm Ehr}\sim\log N\,.$ (169)

As anticipated above such resonances involve selected modes with small momenta
$\mathbf{k}\propto 1/L$, as finite-momentum modes are driven very weakly.
While all the other modes remain weakly excited (at least) for much longer
times $T_{\rm sw}$, the long-wavelength resonant modes form, together with the
collective spin, a full-fledged non-linear dynamical system for intermediate
time scales $T_{\rm Ehr}\ll t\ll T_{\rm sw}$, which will generally feature
semiclassical chaotic behaviour. As the target temporal window lies beyond the
Ehrenfest time scale, the self-consistent approximation expounded here is not
expected to provide a quantitatively accurate description of dynamics in this
regime, and one must resort to numerical simulations to probe this conjectured
chaotic behaviour.

A closely related prethermal regime was discussed by Mori in ††margin:  Quasi-
conservation of permutation operators Ref. [291] from the point of view of the
quasi-conservation of permutation (swap) operators. Even though permutational
symmetry is dynamically broken for $\alpha>0$, Ref. [291] proves the following
bound for long-range interacting quantum Ising chains:

$\left|\langle\hat{P}_{ij}(t)\rangle-\langle\hat{P}_{ij}(0)\rangle\right|\leq|i-j|\,c_{\alpha}\frac{J_{0}t}{N^{1-\alpha}}$
(170)

where $\hat{P}_{ij}$ is the swap operator between spins at sites $i,j$ of the
chain and $c_{\alpha}$ is a $N$-independent constant. [This result can be
easily generalized to arbitrary long-range interacting quantum systems.] As a
result, starting from a permutationally-invariant state with
$\langle\hat{P}_{ij}(0)\rangle=1$ (e.g. fully polarized along a given
direction), one deduces that the symmetry under permuting spins _at a finite
distance_ is approximately conserved, i.e. $\langle\hat{P}_{ij}(t)\rangle\sim
1$, up to time scales $\sim N^{1-\alpha}$. This time scale coincides with
$T_{\rm sw}$ defined above upon setting $d=1$. The underlying physical reason
is that permutation at finite distances probe finite-wavelength fluctuations,
which are consistently frozen over the same time scale. On the other hand,
long-wavelength fluctuations are associated with permutations of spins over
distances proportional to the length $L$ of the spin chain, for which no long-
time quasi-conservation can be guaranteed. Indeed, a severe breakdown of
permutational symmetry may occur over a comparatively short time $T_{\rm
Ehr}\sim\log N$ over such large length scales [cf. the discussion on
resonances above]. For the quench considered in Ref. [291], it can be shown
that spin-wave resonances exist [58]. This observation points to the onset of
a chaotic semiclassical regime over a time scale $T_{\text{Ehr}}\sim\ln N$, in
agreement with the detection of semiclassical chaotic behaviour within the
numerical simulations of prethermal dynamics in Ref. [291].

Finally, it is important to remark that for finite systems, the full self-
consistent system of equations (157) and (161) allows us to investigate non-
linear effects arising from dynamical changes in the driving frequency
triggered by the quantum feedback. Long-time numerical computations suggest
that the stability diagram based on the resonance picture outlined above is
stable to the inclusion of the quantum feedback. It is an open problem to
characterize conditions for long-time stability in long-range interacting spin
systems — both within and beyond the spin-wave approximation discussed here.
An exciting possibility would be to show Kolmogorov-Arnold-Moser-type
stability for this class of systems, which would provide a kind of
interpolation between scenarios in few-body and many-body physics.

Summary: In the strong long-range regime, there exists a long-time scale
during which the finite wavelength excitations are suppressed. A stability
analysis shows that they are stable for generic quenches, while instabilities
for finite $\alpha$ can appear in the proximity of dynamical phase
transitions.

Figure 15: Time-dependent $k$-resolved spin-wave population for $\alpha=0.7$
(right panel) after a quench from $h_{0}=0$ to $h_{f}=2J$. The blue color
gradient for the spin-wave populations in Fourier modes follows the
quasimomentum $\lvert k\rvert$ from the darkest ($k=\pm 2\pi/L$) down to
smaller-wavelength modes with larger $\lvert k\rvert$ (only the first $20$
modes out of $N=500$ are shown). (Right) Density plot of the Kolmogorov-Sinai
entropy rate $\Lambda_{\text{KS}}(\theta_{0},\phi_{0})$ for different initial
conditions $(\theta_{0},\phi_{0})$ on the Bloch sphere for $\alpha=0.7$,
$h=0.5J$. The picture is converged with respect to refining the $k$-space
discretization (here $N=100$). Plots adapted from Ref. [58].

####  4.2.3 Impact of finite-range interactions on dynamical phase transitions

In the last Subsection, we established that for quenches sufficiently far away
from dynamical critical points the dynamical generation of spin-wave
excitations is non-resonant. This observation leads to a long time window
$0<t<T_{\text{sw}}\sim N^{\beta/d}$ (at least) where the spin-wave population
remains low, and the collective spin evolution remains close to the
unperturbed (mean-field) persistent oscillations. In this scenario, the fate
of the mean-field dynamical phase transitions (DPT, discussed above in Sec.
4.1.2) stands out as a naturally prominent question.

This issue has been largely studied as a function of the interaction range
parameter $\alpha$ using various numerical approaches [242, 234, 239, 366,
294, 293]. For the standard long-range quantum Ising chain, a DPT is found in
the thermodynamic limit for $0\leq\alpha<2$, while it is absent for
$\alpha>2$, in qualitative agreement with the equilibrium phase diagram. This
was shown using matrix-product state dynamical simulations in Ref. [242] [see
Fig. 16(a-c)], where the relation between these DPTs and the singularities in
the time-dependence of the Loschmidt echo (DQPTs) has been elucidated, see
also Ref. [234, 239]. This transition has been studied via the semiclassical
truncated Wigner approximation in Ref. [293], where it was found that the
critical exponents for $\alpha\lesssim 0.5$ are the same as the mean-field
DPT. These DPTs have been experimentally observed with trapped-ion quantum
simulators, which simulate versions of the long-range quantum Ising chain [21,
367, 272, 368].

Figure 16:  Dynamical phase transition with finite $\alpha$ in the Ising
Hamiltonian (1) with $\gamma=1$. (a-c) Asymptotic value of the order parameter
in Eq. (128) for different values of the post-quench field $h_{f}/J$ for
quenches from $h_{0}=-$ with (a) $\alpha=0.1$ (b) $\alpha=1.5$ and (c)
$\alpha=3$. Panels adapted from Ref. [242]. (d-f) Experimental data of the
spin-magnetization $\langle\hat{S}^{x}(t)\rangle$ dynamics on a trapped ion
quantum simulator. Data with $L=16$ and $\alpha\simeq 0.8$ and transverse
field $h_{f}/J=B_{z}/J_{0}=0.6,0.,1.6$ in (d,e,d) respectively. Panels adapted
from Ref. [21].

Finite-range ††margin:  New dynamical phases perturbations can also have a
strong impact of the qualitative aspect of the mean-field dynamical phase
diagram by inducing new exotic dynamical phases. Refs. [56, 57] studied the
impact of short-range interactions on top of the fully-connected Ising model,
showing the emergence of a _chaotic dynamical phase_ within which the
asymptotic magnetic ordering is characterized by strong sensitivity to the
parameters and initial conditions. It is found that nonequilibrium
fluctuations can significantly affect the critical dynamics, resulting in a
pseudo-aleatory collective evolution, reminiscent of a particle moving in a
multiple-well potential, with a large initial energy, and subject to friction.
The nonequilibrium phase diagram universally acquires the basic
characteristics of a “coin toss” experiment. This result is confirmed by
matrix-product-state numerical simulations away from the perturbative regime
[56, 57], and a similar scenario was observed in the Dicke model [369].

Summary: Consistent with the prethermal scenario discussed above, dynamical
phases on the two sides of a DPT persist for $0<\alpha\leq 2$, as demonstrated
by numerical results. Semiclassical chaos appears in correspondence of the
dynamical critical point.

####  4.2.4 Scrambling dynamics with variable-range interactions

Let us briefly discuss how a finite value of $\alpha>0$ impacts the scrambling
of quantum information and in particular the OTOC (139) dynamics, introduced
in Sec. 4.1.4. We recall that the square-commutator has been initially put
forward due to its exponential growth, i.e.

$C(t)=\braket{}{[\hat{A}(t),\hat{B}(0)]}{{}^{2}}\simeq\hbar_{\text{eff}}^{2}\,\,e^{\lambda
t}\ ,$

valid before the Ehrenfest or scrambling time $t\ll
T_{\text{Ehr}}\sim\ln\hbar_{\text{eff}}^{-1}\sim\ln N$ for systems with a
semiclassical chaotic limit. This kind of behaviour characterizes as well
large-$N$ all-to-all disordered interacting models [370, 371, 372, 373, 374,
375], despite the absence of an obvious semiclassical limit, including the SYK
which saturates the bound to chaos [302]. However, ††margin:  Fast scrambling
such exponential growth — also known as _fast scrambling_ [376] — is
challenged by finite-range interactions. The square-commutator $C(t)$ was
proved to grow at most polynomially in locally interacting systems [377] and
in long-range interacting systems with $\alpha>d$ [378]. This led to several
proposals suggesting fully-connected interactions as a resource for fast
scrambling [379, 380, 381, 382].

Finite-range interactions lead to a well-defined spatial structure, which
allows investigating the space-dependent square commutator ††margin:  Space
resolved quare-commutator

$C(r,t)=\langle\left|[\hat{A}_{\mathbb{x}}(t),\hat{B}_{\mathbb{x}_{0}}(t)]\right|^{2}\rangle\
,$ (171)

where $r=|{\mathbb{x}}-{\mathbb{x}_{0}}|$ is the distance between the location
of the two considered operators. For locally interacting systems the square
commutator $C(r,t)$ becomes appreciable at times $t\sim x/v_{B}$ [383, 384],
where $v_{B}$, referred to as the _butterfly velocity_ [385], is generally
smaller than the Lieb-Robison one $v_{B}\leq v_{LR}$ [cf. (71)]. By contrast,
long-range interactions are found to induce a non-linear light-cone effect,
whereby information can spread super-ballistically. This occurrence has been
studied numerically in quantum spin chains with variable range interactions
[386, 387, 292, 388] and established via effective hydrodynamics descriptions
in disordered models [389, 390, 391]. All these studies indicate the absence
of ballistic spreading of $C(r,t)$ for $\alpha\leq d$. For quantum spin
systems in the strong long-range regime $\alpha\lesssim d$, one may use the
non-equilibrium spin-wave theory reviewed in Sec. 4.2.1 above to study the
dynamics of the space-resolved square commutator $C(r,t)$. In agreement with
the onset of semiclassical chaos for near-critical quantum quenches, discussed
above in Sec. 4.2.2, one may expect a concomitant exponential growth of the
square commutator in that regime.

Summary: Strong long-range interactions lead to violations of the ballistic
growth of the square-commutator and, in addition, have been proposed as a
resource for fast (exponential) scrambling.

####  4.2.5 Entanglement entropy dynamics: Spin-squeezing vs Quasiparticle
picture

Long-range ††margin:  Counterintuitive slow growth interacting systems exhibit
a conceptually different dynamics of the entanglement entropy (147a) with
respect to locally interacting systems [112, 113, 325, 326, 114]. On one side,
their non-local interactions allow quantum correlations between distant
degrees of freedom to build up very quickly. As discussed above, this leads to
violations of the Lieb-Robinson bound (71) and nonlinear light-cone spreading
of quantum correlations, see Sec. 3.1 and 4.2.4. On the other hand, the
bipartite entanglement entropy growth after a quench with the Hamiltonian (1)
was found to exhibit a counterintuitive dramatic slowdown as the range of
interactions is increased: It becomes logarithmically slow for algebraically-
decaying couplings with $\alpha$ smaller than the spatial dimensionality $d$
[125, 127, 292, 58], see also Fig. 17. Such numerical results can be
rationalized using the semiclassical techniques introduced above, which leads
to a complete picture of entanglement growth for long-range interacting
systems [58]. Reviewing this framework is the goal of the present Subsection.
For completeness, we mention in passing that multipartite entanglement
associated with algebraically decaying interactions has been studied in depth,
e.g. in the form of dynamical spin-squeezing [392, 393, 394, 395, 396, 397] or
via its relation to dynamical susceptibilities in equilibrium [398].

Figure 17: Logarithmic growth of the entanglement entropy after a quench with
the Ising Hamiltonian (1) with $\gamma=1$, $d=1$. (a) Exponential of
$S_{A}(t)$ as a function of dimensionless time, for different transverse
fields $h=0.7J,1J,1.3J$. For each $h$, the results for $L=30,40,50$ and
$\alpha=0.8,0.9,1$ are plotted. (b) $S_{A}(t)$ in logarithmic scale for
different initial states
$|\psi_{0}\rangle=|\uparrow\uparrow\dots\uparrow\rangle$ and the ones
generated by applying single site Pauli operators. Simulation with $L=50$,
$\alpha=0.5$ and $h=J$. Image adapted from Refs. [127, 125] for (a), (b)
respectively.

In the fully-connected limit $\alpha=0$, the growth of entanglement is
determined by the squeezing of the collective fluctuations stemming from the
underlying classical trajectory, see Fig. 5. The general framework illustrated
in Sec. 4.1.5 predicts logarithmic growth in the absence of semiclassical
chaos (Table 1), which is generic in fully-connected spin systems222222without
self-interactions if $s>1/2$.

For finite $\alpha$ the behavior of $S_{A}(t)$ can be understood at
intermediate times by accounting for spin-wave excitations with non-vanishing
momentum $\mathbf{k}$ [cf. Eq.(160)] on top of the entanglement dynamics
arising from collective spin excitations (or spin squeezing), discussed in
Sec. 4.1.5 above. As ††margin:  Spin waves entanglement entropy described in
Sec. 4.2.1, the time-evolving state of the spin-wave excitations is encoded in
the correlations in Eq.(160), i.e.
$({G}^{qq}_{\mathbf{k}}(t),{G}^{qp}_{\mathbf{k}}(t),{G}^{pp}_{\mathbf{k}}(t))$
defined by
${G}^{\alpha\beta}_{\mathbf{k}}(t)=\frac{1}{2}\Big{\langle}\tilde{\alpha}_{\mathbf{k}}(t)\tilde{\beta}_{-\mathbf{k}}(t)+\tilde{\beta}_{\mathbf{k}}(t)\tilde{\alpha}_{-\mathbf{k}}(t)\Big{\rangle}$
for $\alpha,\beta=q,p$. Within the linear spin-wave analysis, the state of a
subsystem composed of $M=f_{A}N<N$ spins contained in a region $A$ of the
lattice is a Gaussian bosonic state determined by the instantaneous
correlations

$\bigg{\\{}G^{\alpha\beta}_{\mathbf{r},\mathbf{r^{\prime}}}(t)=\big{\langle}\alpha_{\mathbf{r}}(t)\beta_{\mathbf{r^{\prime}}}(t)+\beta_{\mathbf{r}}(t)\alpha_{\mathbf{r^{\prime}}}(t)\big{\rangle}\bigg{\\}}_{\begin{subarray}{c}\mathbf{r},\mathbf{r^{\prime}}\in
A\\\ \alpha,\beta=q,p\end{subarray}}$ (172)

within $A$, which can be expressed in terms of
$\tilde{G}^{\alpha\beta}_{\mathbf{k}}(t)$ via Fourier antitransform. This set
of correlations uniquely identifies the reduced density matrix
$\hat{\rho}_{A}(t)$. The von Neumann entropy of this Gaussian bosonic state
can be computed via standard techniques [331], namely

$S_{A}=\sum_{i=1}^{M}\,S(\nu_{i})\ ,\quad\text{with}\quad
S(\nu_{i})=\frac{\nu_{i}+1}{2}\,\ln\frac{\nu_{i}+1}{2}-\frac{\nu_{i}-1}{2}\,\ln\frac{\nu_{i}-1}{2}\
,$ (173)

where $\nu_{i}$ are the symplectic eigenvalues of the correlation matrix.

Figure 18: Entanglement dynamics after a quench from the ferromagnetic ground
state $h_{0}=0$ with a long-range Ising Hamiltonian (1) with $\gamma=1$ and
$h_{f}=2$. Comparison between finite-size MPS-TDVP numerical data (light-to-
dark blue curves for increasing $N$), the spin-squeezing contribution (grey)
and full spin-wave entanglement (black), for $\alpha=0.1$ (left panel) and
$0.7$ (right panel), for the quench $h_{0}=0\to h_{f}=2J$, with $N=500$.
Figure adapted from Ref. [58].

For long-range interactions with $0<\alpha<d$, the growth of $S_{A}(t)$ turns
out to be determined by the stability of the discrete set of long-wavelength
excitations, expressed by the Floquet quasi-frequencies $\lambda_{\mathbf{k}}$
with $\lvert\mathbf{k}\rvert\propto 1/L$: see the dedicated discussion in Sec.
4.2.2 and Fig. 15. In particular, one can apply the general semiclassical
description entanglement discussed in Sec. 4.1.5 and summarized in Table 1. If
all the modes are stable (i.e., non-resonant), then $S_{A}(t)\sim\ln t$
exhibits a slow growth dominated by the collective spin fluctuations with
$\mathbf{k}=\mathbf{0}$ only. This is indeed the case for typical quenches
away from dynamical criticality, as discussed in Sec. 4.2.2. This observation
underlies and rationalizes the previous numerical findings of logarithmic
growth of the von Neumann entanglement entropy [125, 127, 292], reported at
the beginning of this Subsection. On the other hand, if some mode is unstable
(i.e., resonant), then ${S_{A}(t)\sim\Lambda_{\text{KS}}\,t}$ exhibits a fast
growth dominated by the instabilities, with $\Lambda_{\text{KS}}$ in Eq. 168.
This is what may happen for quenches in the proximity of dynamical critical
points, discussed in Sec. 4.1.2.

The ††margin:  Long-range picture _physical picture_ for the long-range
entanglement dynamics is now clear before the Ehrenfest time:

  * •

The leading contribution comes from the semi-classical squeezing of the
collective spin, which grows logarithmically in the absence of classical
chaos;

  * •

In the strong long-range regime, the suppressed long-wavelength spin waves
provide a subleading contribution to the entanglement growth.

The above analysis shows that slow logarithmic growth of the entanglement
entropy should be generally expected in quench dynamics of spin systems with
strong long-range interactions starting from a state with large spin
polarization232323Subject to the usual caveat of the absence of individual
spin self-interactions for $s>1/2$..

One ††margin:  Ising model example can solve the spin-wave equation of motion
in Eqs (161) and compute the resulting time-dependent entanglement entropy via
Eq. (173). The results for a typical quench in the long-range quantum Ising
chain in a transverse field [cf. Eq. (1) with $\gamma=1$, $d=1$] are shown in
Fig. 18, where the exact numerical $S_{A}(t)$ for finite system size is
compared with fully-connected “spin-squeezing” contribution and with the
result obtained with the inclusion of spin waves. This analysis applies to a
wide variety of spin models and quenches, as shown in E.1 where we study the
long-range Ising Hamiltonian with transverse and longitudinal field for a
quench near the critical point.

We remark that the underlying mechanism crucially relies on the _discreteness_
of the set of excitation modes (the long-wavelength spin waves) which result
in a bounded, subleading contribution to entanglement growth. This property
††margin:  Other types of perturbations is characteristic of strong long-range
interactions ($\alpha<d$) and generically does not occur for other types of
perturbations. If, for instance, a _finite-range_ perturbation is added on top
of a fully-connected model, one can still have stable excitations. However,
the presence of the continuous spectrum of excitations results in light-cone
spreading of quantum correlations and linear growth of entanglement according
to a standard quasiparticle picture [111], see e.g. [57] and E.2 for an
example. In the weak long-range regime $d<\alpha<d+2$, the growth of
entanglement has been related to the nonlinear dispersion relation of
quasiparticles [104].

We finally reiterate that the picture of entanglement dynamics for long-range
interacting spin systems reviewed here, based on the semiclassical dynamics of
quantum spin fluctuations, covers setups not encompassed by other theoretical
pictures such as quasi-particles [111], spacetime membranes [114] or local
integrals of motion [399].

Summary: For $0\leq\alpha<d$, a semi-classical picture predicts that the
entanglement growth is dominated by the collective spin squeezing (logarithmic
for generic quenches) and decorated by a discrete set of stable excitations.

##  5 Dynamical phases induced by periodic driving

In this Section we will expand our analysis to non-autonomous, coherently
driven systems. We will show how the previously introduced ideas allow to
characterize nonequilibrium phases of spin systems with novel kinds of
collective order dynamically stabilized by a periodic drive, which would not
be possible in equilibrium [400]. Here, long-range interactions play the
twofold role of protecting long-range order in highly excited states and
hindering heating. We explain how this basic mechanism also protects
spatiotemporal order such as time-crystalline behavior [401, 26, 402].

###  5.1 Kapitza phases

As realized by Kapitza long ago, a rigid pendulum can be stabilized upside
down by periodically driving its suspension point with tuned amplitude and
frequency. While this dynamical stabilization is feasible in a variety of