# Spreading of Correlations and Entanglement in the Long-Range Transverse
Ising Chain

J. T. Schneider  CPHT, CNRS, Ecole Polytechnique, IP Paris, F-91128 Palaiseau,
France  J. Despres  CPHT, CNRS, Ecole Polytechnique, IP Paris, F-91128
Palaiseau, France  S. J. Thomson  CPHT, CNRS, Ecole Polytechnique, IP Paris,
F-91128 Palaiseau, France  L. Tagliacozzo  Department of Physics and SUPA,
University of Strathclyde, Glasgow G4 0NG, United Kingdom  Departament de
Física Quàntica i Astrofísica and Institut de Ciències del Cosmos (ICCUB),
Universitat de Barcelona, Martí i Franquès 1, 08028 Barcelona, Catalonia,
Spain  L. Sanchez-Palencia  CPHT, CNRS, Ecole Polytechnique, IP Paris, F-91128
Palaiseau, France

(August 29, 2024)

###### Abstract

Whether long-range interactions allow for a form of causality in non-
relativistic quantum models remains an open question with far-reaching
implications for the propagation of information and thermalization processes.
Here, we study the out-of-equilibrium dynamics of the one-dimensional
transverse Ising model with algebraic long-range exchange coupling. Using a
state of the art tensor-network approach, complemented by analytic
calculations and considering various observables, we show that a weak form of
causality emerges, characterized by non-universal dynamical exponents. While
the local spin and spin correlation causal edges are sub-ballistic, the causal
region has a rich internal structure, which, depending on the observable,
displays ballistic or super-ballistic features. In contrast, the causal region
of entanglement entropy is featureless and its edge is always ballistic,
irrespective of the interaction range. Our results shed light on the
propagation of information in long-range interacting lattice models and pave
the way to future experiments, which are discussed.

Long-range interactions may dramatically impact the dynamics of correlated
systems [1]. In the quantum regime, a number of basic concepts, such as the
equivalence of the thermodynamic ensembles [2, *kastner2010b], the Mermin-
Wagner-Hohenberg theorem [4, 5, 6] or the area law for entanglement entropy
[7, 8, 9, 10], break down. Long-range interactions may also be responsible for
negative heat capacities and anomalous response functions [11]. The
paradigmatic model considers interactions that fall off algebraically with the
distance $R$, as $1/R^{\alpha}$. Such long-range interactions can now be
emulated in quantum simulators [12, 13, 14, *NaturePhysicsInsight2012bloch,
*NaturePhysicsInsight2012blatt, *NaturePhysicsInsight2012aspuru-guzik,
*NaturePhysicsInsight2012houck, 19, 20, *tarruell2018, *aidelsburger2018,
*lebreuilly2018, *LeHur2018, *bell2018, *alet2018] with artificial ion
crystals [27, 28, 29, 30, 16, 31], Rydberg gases [32, 33, 34, 35], magnetic
atoms [36, 37, 38, 39, 40], polar molecules [41, 42, 43], nonlinear optical
media [44], and solid-state defects [45, 46, 47]. A major asset of these
systems is that they are free of screening effects and the interaction is
truly long range. Moreover, the exponent $\alpha$ can be controlled [27, 35].

Another dramatic effect of long-range interactions is the breakdown of the
notion of causality in non-relativistic quantum models. For a wide class of
lattice models with short-range interactions, the Lieb-Robinson bound implies
that correlations decay exponentially beyond a limit set by the system’s
maximum group velocity [48, 49, 50]. This effective light-cone effect has been
demonstrated for various models in both experiments [51, 52, 53] and numerics
[54, 55, 56, 57, 58, 59, 60]. In contrast, for sufficiently long-range
interactions, causality breaks down and information can propagate arbitrarily
fast [61, 62]. This is consistent with the absence of known Lieb-Robinson
bounds [49, 63, 64] and the vanishing of the characteristic dynamical time
scale in the thermodynamic limit [65]. Conversely, algebraic interactions that
decay fast enough are effectively short range, and one recovers ballistic
propagation of information [66, 67, 68]. The intermediate regime is, however,
strongly debated, and whether a form of causality emerges remains an open
question. Known bounds in principle allow for super-ballistic propagation [49,
63, 64] but they are challenged by numerical simulations, which point towards
a significantly slower propagation [61, 10, 69, 70, 71]. The latter, however,
reported different propagation scaling laws. Microscopic mean-field theory
suggests that these apparent contradictions may be attributed to the
coexistence of several signals governed by different dynamical scaling laws
[72]. This prediction, however, relies on a generic but non-universal form of
the correlation functions and ignores beyond-mean-field effects. Experiments
performed with trapped ions have reported bounded propagation [73, 74] but
they are limited to very small systems, which prevents extraction of the
scaling laws and closure of the debate.

The aim of this work is to characterize the spreading of quantum correlations
in the intermediate regime of a long-range spin system numerically.
Specifically, we determine the scaling laws for the propagation of a variety
of observables in the long-range transverse Ising (LRTI) chain using matrix
product state simulations. We find that a weak form of causality emerges,
characterized by generic algebraic scaling forms $t\sim R^{\beta}$, where the
specific value of the exponent $\beta$ depends on the observable. The spin-
correlation and local-spin causal edges are both sub-ballistic, with the same
dynamical exponent ($\beta_{\textrm{\tiny CE}}=\beta_{\textrm{\tiny SE}}>1$).
In the vicinity of the edge, however, the local maxima propagate differently,
i.e. ballistically ($\beta_{\textrm{\tiny m}}=1$) for spin-spin correlations
and super-ballistically ($\beta_{\textrm{\tiny m}}<1$) for local spins. In
contrast, the Rényi entanglement entropies always propagate ballistically,
irrespective of the range of interactions, in both the local and quasi-local
regime. The analytic quasi-particle picture, based on linear spin wave theory
(LSWT), accurately reproduces the numerics and provides a clear interpretation
of the numerical results. The different algebraic space-time patterns of
correlation functions provide an unambiguous fingerprint of the dynamical
regimes of the model, suggesting the emergence of a dynamical phase-diagram.
These correlation patterns can be directly measured in state-of-the-art
experiments.

#### Model and approach.—

The dynamics of the LRTI chain is governed by the Hamiltonian

$\hat{H}=\sum_{R\neq
R^{\prime}}\frac{J}{|R-R^{\prime}|^{\alpha}}\hat{S}^{x}_{R}\hat{S}^{x}_{R^{\prime}}-2h\sum_{R}\hat{S}^{z}_{R},$
(1)

where $\hat{S}_{R}^{j}$ ($j=x,y,z$) are the spin-$1/2$ operators on lattice
site $R\in[0,N-1]$, $N$ is the system size, $J>0$ is the coupling energy, and
$h$ is the transverse field. It can be realized on various quantum simulation
platforms, including cold Rydberg gases [75, 76, 35] and artificial ion
crystals, where the exponent $\alpha$ can be controlled via light-mediated
interactions [27, 29, 30, 74, 73]. At equilibrium, the phase diagram of the
LRTI chain comprises two gapped phases separated by a second order quantum
phase transition, see Fig. 1 and Ref. [8]. For low fields $h$ and rather
short-range couplings (high values of $\alpha$), the nearest-neighbor anti-
ferromagnetic couplings dominate and the system forms a staggered Néel-ordered
phase along the $x$ direction. For a large field $h$ and long-range couplings
(low values of $\alpha$), the spin-field interaction is favoured and a
$z$-polarized phase is formed.

Figure 1:  Phase diagram and dynamical properties of the LRTI chain in the
$(h/J,\alpha)$ plane. The $x$-Néel phase (shades of blue) is separated from
the $z$-polarized phase (shades of orange) by a critical line in black. The
three dynamical regimes are shaded differently. As shown in the insets, the
local and quasi-local regimes are characterized by distinct algebraic scaling
laws for the correlation and spin edges (CE and SE, respectively), their
maxima in space-time (max), and the entanglement edge (EE), after global or
local quench. The different behaviors of these observables in different
regions of the phase diagram provide clear, experimentally accessible,
signatures of the local and quasi-local regimes.

Out of equilibrium, the LSWT predicts three dynamical regimes (shaded by
different colors in Fig. 1) [61, 69]: For $\alpha\geq 2$ (local regime), the
spin wave excitations are regular with bounded energies $E_{k}$ and group
velocities $V_{\textrm{g}}(k)=\partial_{k}E_{k}$ ($k$ is the momentum). This
regime is reminiscent of the short-range case where correlations spread at
finite speed, giving rise to a linear causality cone [48, 49]. For $\alpha<1$
(instantaneous regime), $E_{k}$ features an algebraic, infrared divergence.
There is no characteristic time scale and correlations spread arbitrarily
fast. Finally, for $1\leq\alpha<2$ (quasi-local regime), $E_{k}$ is bounded
but $V_{\textrm{g}}(k)$ diverges as $V_{\textrm{g}}(k)\sim 1/k^{2-\alpha}$.
Whether some form of causality emerges in the quasi-local regime remains
debated [61, 74, 73, 69, 77, 65, 71], and, in the following, we mainly focus
on this case.

Except where otherwise indicated, the results discussed below are obtained
using time-dependent variational principle (TDVP) simulations within a matrix-
product state (MPS) framework [78, 8, 61]. Convergence of the calculations
with the bond dimension has been systematically checked. Our results are
summarized on the dynamical phase diagram of Fig. 1. The main features of the
observable that we consider are all described by an algebraic space-time
dependence of their edges and maxima, shown in the insets.

#### Spin correlations and global quench.—

We first consider the spreading of the connected spin correlations,
$G_{j}(R,t)=G_{j}^{0}(R,t)-G_{j}^{0}(R,0)$, with
$G_{j}^{0}(R,t)=\langle\hat{S}^{j}_{R}(t)\hat{S}^{j}_{0}(t)\rangle-\langle\hat{S}^{j}_{R}(t)\rangle\langle\hat{S}^{j}_{0}(t)\rangle$,
along the directions $j=x,z$. Figure 2(a) shows a typical TDVP result for
$G_{z}(R,t)$, for a quench in the quasi-local regime, $\alpha=1.7$, from
$(h/J)_{\textrm{i}}=50$ to $(h/J)_{\textrm{f}}=1$, both in the $z$-polarized
phase. Similar results are found when scanning the values of
$(h/J)_{\textrm{i}}$, $(h/J)_{\textrm{f}}$, and $1<\alpha<2$. The correlation
pattern shows a series of maxima which propagate algebraically [notice the
log-log scale in Fig. 2(a)]. In the asymptotic long-time and long-distance
limits, they are well fitted by the scaling law $t\sim R^{\beta_{\textrm{\tiny
m}}}$ (dashed blue lines). The correlation edge (CE), which sets the causality
horizon, cannot, however, be inferred from the behavior of the maxima [58]. To
find the CE, we track the points of the $R-t$ plane where the correlations
reach a fraction $\epsilon$ of their maximum. Scanning
$\epsilon\in[0.01,0.1]$, we find that the CE is well fitted by the algebraic
scaling law $t\sim R^{\beta_{\textrm{\tiny CE}}}$ (solid line), with
$\beta_{\textrm{\tiny CE}}$ nearly independent of $\epsilon$ [79]. Figure 2(b)
shows the results of the fits (empty symbols) versus the exponent $\alpha$.

The numerical results for $G_{x}(R,t)$ have an extra checkerboard structure
inside the causal region due to the antiferromagnetic coupling in Eq. (1).
Once this structure is taken into account, we can identify a CE and local
maxima, and extract the corresponding exponents. The results are similar to
those found for $G_{z}(R,t)$, see filled symbols in Fig. 2(b) [79].

Figure 2:  Spreading of spin correlations in a global quench. (a) TDVP results
for $G_{z}$ and a quench from $(h/J)_{\mathrm{i}}=50$ to
$(h/J)_{\mathrm{f}}=1$ in the quasi-local regime at $\alpha=1.7$. The solid
green and dashed blue lines are fits to the CE and the extrema, respectively.
(b) Dynamical exponents $\beta_{\textrm{\tiny CE}}$ (green diamonds) and
$\beta_{\textrm{\tiny m}}$ (blue disks), fitted form results as in panel (a)
for $G_{z}$ (empty symbol) and $G_{x}$ (full symbol), and comparison to the
LSWT predictions (solid green and dashed blue lines). (c) Spreading velocities
$V_{\textrm{\tiny CE}}$ (green diamonds) and $V_{\textrm{m}}$ (blue disks), in
the local regime and comparison to the LSWT predictions (solid green and
dashed blue lines).

Comparing the fitted dynamical exponents $\beta_{\textrm{\tiny m}}$ and
$\beta_{\textrm{\tiny CE}}$ to the predictions of the LSWT [72], we find an
excellent agreement, as shown in Fig. 2(b). While the maxima spread
ballistically, $\beta_{\textrm{\tiny m}}\simeq 1$, the CE is sub-ballistic
with $\beta_{\textrm{\tiny CE}}\simeq\beta_{\textrm{\tiny LSWT}}=3-\alpha>1$.
This is characteristic of gapped long-range models, such as the LRTI model in
the z-polarized phase [72]. It confirms the emergence of a weak, slower-than-
ballistic form of causality for $1<\alpha<2$.

In the local regime, $\alpha\geq 2$, we find that both the maxima and the CE
spread ballistically, $\beta_{\textrm{\tiny m}}\simeq\beta_{\textrm{\tiny
CE}}\simeq 1$. The spreading velocities are, however, different from each
other, which is characteristic of a non-phononic excitation spectrum. The CE
velocity is $V_{\textrm{\tiny CE}}\simeq 2V_{\mathrm{g}}(k^{*})$, i.e. twice
the maximum group velocity of the spin waves, while that of the maxima is
$V_{\textrm{m}}\simeq 2V_{\varphi}(k^{*})$, i.e. twice the phase velocity
$V_{\varphi}(k)=E_{k}/k$ at the quasi-momentum $k^{*}$ where the group
velocity is maximum, see Fig. 2(c).

We have also performed calculations for large quenches. While the LSWT is well
justified only for weak quenches, we have found that the dynamical exponents
it predicts are extremely robust, even for quenches across the critical line.
We found dynamical scaling laws in very good agreement with those reported on
Figs. 2(b) and (c), although the signal is blurred as compared to weak
quenches owing to the proliferation of quasi-particles when the quench
amplitude increases [79].

#### Local magnetization and local quench.—

We now consider the dynamics of another quantity, namely the local
magnetization $\langle\hat{S}^{z}_{R}(t)\rangle$, and perform a local quench.
We initialize the system in the ground state of the $z$-polarized phase and
flip the central spin. Figure 3(a) shows a typical TDVP result for the
quantity $1/2-\langle\hat{S}^{z}_{R}(t)\rangle$ versus the time $t$ and the
distance $R$ from the flipped spin, in the quasi-local regime. The result
displays, as in the case of global quenches, a twofold algebraic structure.
Fitting the maxima (dashed blue line) and the spin edge (SE, solid green line)
as previously, we extract the dynamical exponents $\beta_{\textrm{\tiny m}}$
and $\beta_{\textrm{\tiny SE}}$ plotted in Fig. 3(b) (blue disks and green
diamonds, respectively). The SE follows the same sub-ballistic scaling law as
the CE in the previous case, $\beta_{\textrm{\tiny SE}}\simeq 3-\alpha$ (solid
green line). In contrast, the maxima of the local spin spreads faster than
those of the spin correlations, and here we find $\beta_{\textrm{\tiny m}}<1$,
corresponding to a super-ballistic propagation.

Figure 3:  Spreading of the local magnetization for a local quench. (a) TDVP
results for the quantity $1/2-\langle\hat{S}^{z}_{R}(t)\rangle$ versus the
time and the distance from the flipped spin for $h/J=50$ and $\alpha=1.8$. (b)
Dynamical exponents $\beta_{\textrm{\tiny SE}}$ (green diamonds) and
$\beta_{\textrm{\tiny m}}$ (blue disks), fitted form results as in panel (a)
and comparison to the LSWT predictions (solid green and dashed blue lines).
(c) Spreading velocity $V_{\textrm{\tiny SE}}$ (green diamonds) and comparison
to the LSWT prediction in the local regime (solid green line).

To understand this behavior, let us first recall that the standard LSWT
generically predicts a super-ballistic propagation of the maxima in gapless
systems [72]. Here, although the LRTI model is gapped, the initial state
$\ket{\Psi_{0}}\simeq\ket{\uparrow...\uparrow\,\downarrow\,\uparrow...\uparrow}$
of the local quench is orthogonal to the ground state
$\ket{\Psi_{\textrm{\tiny GS}}}\simeq\ket{\uparrow\,\uparrow...\uparrow}$, and
thus lives in the first excited manifold. The ground state and the gap are
thus irrelevant, and we may expect super-ballistic propagation of the maxima
consistent with the TDVP results. This argument applies to any observable for
such a local quench. More specifically, the LSWT applied to
$\langle\hat{S}^{z}_{R}(t)\rangle$ for the local quench yields
$\beta_{\textrm{\tiny m}}=\alpha-1$, in very good agreement with the TDVP
results, see Fig. 3(b) [79].

In the local regime ($\alpha>2$) we find that
$\langle\hat{S}^{z}_{\tilde{R}}(t)\rangle$ propagates ballistically. Both this
property and the SE velocity extracted from the TDVP calculations are in good
agreement with the LSWT analysis, see Fig. 3(c). Note that in this regime, we
do not observe maxima propagating at a different velocity. This is also
consistent with the LSWT analysis. In the local regime, the quantity
$\langle\hat{S}^{z}_{R}(t)\rangle$ is the sum of several contributions, each
with a twofold structure but that are in phase quadrature and cancel each
other. A similar effect has been found for density correlations deep in the
Mott insulator phase of the Bose-Hubbard model [58, 79].

#### Entanglement entropy.—

We finally study the spreading of quantum information after a local quench. It
may be measured via the Rényi entropies of the reduced density matrix of a
block,

$\mathcal{S}_{n}(R,t)=\frac{1}{1-n}\log\\!\big{\\{}\tr[\hat{\rho}^{n}(R,t)]\big{\\}}\,,$
(2)

with $n>0$ and
$\hat{\rho}(R,t)=\hat{\rho}_{A}=\tr_{B}(\outerproduct{\Psi(t)}{\Psi(t)})$ the
reduced density matrix at time $t$ of the subsystem $A=[R,R+1,\ldots,N/2]$.
The position $R$ is measured from the flipped central spin and $B$ denotes the
complementary subsystem, see Inset in Fig. 4(a).

Figure 4(a) shows a typical TDVP result for the von Neumann entropy,
$n\rightarrow 1$ 111The von Neumann entropy also reads as
$\mathcal{S}_{\textrm{vN}}=\mathcal{S}_{n\rightarrow
1}(R,t)=-\mathrm{Tr}\left\\{\hat{\rho}(R,t)\mathrm{log}\left[\hat{\rho}(R,t)\right]\right\\}$,
and the same local quench as in Fig. 3(a). Similar results are found for the
other Rényi entropies [79]. We find that $\mathcal{S}_{n}(R,t)$ is a monotonic
function of both position and time, and no local maxima such as those found
previously for the correlation function and the magnetization are observed.
This is consistent with a causal spreading of information and the expectation
that the entanglement entropy decreases with the size of the smaller partition
($A$). The entanglement edge (EE) is clearly visible in Fig. 4(a) and a fit to
the algebraic scaling law $t\sim R^{\beta_{\textrm{\tiny EE}}^{n}}$ allows us
to extract the exponent. Varying the threshold in a wide range,
$\epsilon\in[0.2,0.8]$, we find $\beta_{\textrm{\tiny EE}}^{n}\approx 1$
within error bars, irrespectively of the interaction range, see cyan points in
Fig. 4(b). The error bars correspond to the variation of $\beta_{\textrm{\tiny
EE}}^{n}$ with $\epsilon$, which is due to finite-size effects (see below).
Our results imply that the propagation of entanglement is close to ballistic
in both the local and quasi-local regimes. This contrasts with the behavior of
the one- and two-point observables considered thus far, which exhibit sub-
ballistic causal edges. It is a direct consequence of the fact that the
bipartite entanglement entropy is a highly non-local quantity, which takes
into account all entangled pairs on either side of the bipartition $R$ [50,
81, 82, 83].

More precisely, we have analytically computed the entanglement entropy within
LSWT [79]. The results are in good agreement with those extracted via TDVP for
the same system sizes, see Fig. 4(b). The ballistic spreading of entanglement
is further confirmed by the LSWT results obtained in much larger systems, with
significantly smaller error bars, see Fig. 4(c). Bipartitioning the quenched
initial state within LSWT yields an entanglement Schmidt rank bounded by $2$.
This is consistent with our TDVP results at all times and the saturation of
the entanglement entropies to
$\mathcal{S}_{n}(R,t\rightarrow\infty)\simeq\log(2)\simeq 0.69$, see Fig.
4(c). The two eigenvalues of the reduced density matrix $\hat{\rho}_{A}$,
$\lambda_{1}(R,t)$ and $\lambda_{2}(R,t)=1-\lambda_{1}(R,t)$, can then be
computed analytically. In the asymptotic limit and for not too small values of
$R/t$, we find $\lambda_{2}(R,t)\propto
t^{\frac{1}{2-\alpha}}\zeta(\frac{3-\alpha}{2-\alpha},R)\sim(t/R)^{\frac{1}{2-\alpha}}$,
with $\zeta$ the Hurwitz zeta function [84, 85]. Hence, the $n$-order Rényi
entropy is a function of the ratio $R/t$,

$\mathcal{S}_{n}(R,t)\simeq\frac{1}{1-n}\log\\!\big{\\{}\lambda_{1}\\!\left({R}/{t}\right)^{n}+\lambda_{2}\\!\left({R}/{t}\right)^{n}\big{\\}}\,.$
(3)

This confirms the ballistic propagation of the entanglement entropy
($\beta_{\textrm{\tiny EE}}^{n}=1$) consistent with the results of Fig. 4 for
$n=1$ and other Rényi orders $n$ [79].

Figure 4:  Spreading of the von Neumann entanglement entropy
$\mathcal{S}_{n=1}(R,t)$ for the same quench as in Fig. 3(a). (a) TDVP results
for $h/J=50$, $\alpha=1.8$, and a system size $N=96$. The solid green line
marks a power law fit to the EE with $\epsilon=0.5$, yielding
$\beta_{EE}^{n}=0.899\pm 0.005$. (b) Dynamical exponents for the EE obtain via
TDVP (cyan downwards triangles) and LSWT (magenta upwards triangles), with
error bars corresponding to the variations with respect to the threshold
$\epsilon$. (c) Dynamical exponents obtained via LSWT for a larger system,
$N=512$.

#### Conclusion.—

Our results show the emergence of a weak form of causality in the intermediate
regime of the long-range Ising model, characterized by algebraic propagation
laws with exponents that depend on the observables and the range of
interactions. While local spins and spin correlations both have a sub-
ballistic propagation edge, $t\propto R^{\beta}$, with $\beta>1$, the causal
region is characterized by local maxima propagating super-ballistically and
ballistically, respectively. The distinction between the causal edge and the
local maxima, which can show drastically different dynamical behaviors, is
thus pivotal in the characterization of causality in long-range quantum
systems. In contrast, the propagation of entanglement is ballistic in both the
local and quasi-local regimes, and the causal region is featureless.

These results call for future experimental and theoretical work. On the one
hand, our predictions are directly relevant to quantum simulators using for
instance trapped ions, where the interaction range can be controlled. While
analysis of spin and correlation spreading in first experiments have been
limited by finite-size effects [74, 73], systems with more than 50 ions,
comparable to the system size used in our simulations, are now accessible [86]
and Rényi entanglement entropies can now be measured [87, 88]. On the other
hand, it would be interesting to further test the robustness of the observed
algebraic scaling laws by quantitatively investigating the dependence (if any)
of the exponents on the strength of the quenches, the phases of the model and
their values in different models, such as the long-range XY and Hubbard
models, as well as in dimensions higher than one. This extended analysis could
allow the identification of dynamical universality classes, i.e. models which
share the same algebraic laws for correlations out-of-equilibrium.

We thank the CPHT computer team for valuable support. The numerical
calculations were performed using HPC resources from GENCI-CINES (Grant
2019-A0070510300). We acknowledge use of the ITensors.jl software package
[89]. L.T. acknowledges support from the Ramón y Cajal program RYC-2016-20594,
the “Plan Nacional Generación de Conocimiento” PGC2018-095862-B-C22 and the
State Agency for Research of the Spanish Ministry of Science and Innovation
through the “Unit of Excellence María de Maeztu 2020-2023” award to the
Institute of Cosmos Sciences (CEX2019-000918-M).

J.T.S. and J.D. contributed equally to this work.

## References

  * [1] T. Dauxois, S. Ruffo, E. Arimondo, and M. Wilkens, eds., _Dynamics and thermodynamics of systems with Long Range Interactions_ , Lecture Notes in Physics (Springer-Verlag, Berlin Heidelberg, 2002).
  * [2] M. Kastner, _Nonequivalence of ensembles for long-range quantum spin systems in optical lattices_ , Phys. Rev. Lett. 104, 240403 (2010).
  * [3] M. Kastner, _Nonequivalence of ensembles in the Curie-Weiss anisotropic quantum Heisenberg model_ , J. Stat. Mech. 2010(07), P07006 (2010).
  * [4] N. D. Mermin and H. Wagner, _Absence of ferromagnetism or antiferromagnetism in one- or two-dimensional isotropic Heisenberg models_ , Phys. Rev. Lett. 17, 1133 (1966).
  * [5] P. C. Hohenberg, _Existence of long-range order in one and two dimensions_ , Phys. Rev. 158, 383 (1967).
  * [6] D. Peter, S. Müller, S. Wessel, and H. P. Büchler, _Anomalous behavior of spin systems with dipolar interactions_ , Phys. Rev. Lett. 109, 025303 (2012).
  * [7] J. Eisert, M. Cramer, and M. B. Plenio, _Colloquium: Area laws for the entanglement entropy_ , Rev. Mod. Phys. 82, 277 (Feb 2010).
  * [8] T. Koffel, M. Lewenstein, and L. Tagliacozzo, _Entanglement entropy for the long-range Ising chain in a transverse field_ , Phys. Rev. Lett. 109, 267203 (2012).
  * [9] A. Cadarso, M. Sanz, M. M. Wolf, J. I. Cirac, and D. Pérez-García, _Entanglement, fractional magnetization, and long-range interactions_ , Phys. Rev. B 87, 035114 (2013).
  * [10] J. Schachenmayer, B. P. Lanyon, C. F. Roos, and A. J. Daley, _Entanglement growth in quench dynamics with variable range interactions_ , Phys. Rev. X 3, 031015 (2013).
  * [11] W. Thirring, _Systems with negative specific heat_ , Z. Phys. 235, 339 (1970).
  * [12] I. Buluta and F. Nori, _Quantum simulators_ , Science 326(5949), 108 (2009).
  * [13] I. M. Georgescu, S. Ashhab, and F. Nori, _Quantum simulation_ , Rev. Mod. Phys. 86, 153 (2014).
  * [14] J. I. Cirac and P. Zoller, _Goals and opportunities in quantum simulation_ , Nat. Phys. 8, 264 (2012).
  * [15] I. Bloch, J. Dalibard, and S. Nascimbène, _Quantum simulations with ultracold quantum gases_ , Nat. Phys. 8, 267 (2012).
  * [16] R. Blatt and C. F. Roos, _Quantum simulations with trapped ions_ , Nat. Phys. 8, 277 (2012).
  * [17] A. Aspuru-Guzik and P. Walther, _Photonic quantum simulators_ , Nat. Phys. 8, 285 (2012).
  * [18] A. A. Houck, H. E. Tureci, and J. Koch, _On-chip quantum simulation with superconducting circuits_ , Nat. Phys. 8, 292 (2012).
  * [19] C. Gross and I. Bloch, _Quantum simulations with ultracold atoms in optical lattices_ , Science 357, 995 (2017).
  * [20] L. Sanchez-Palencia, _Quantum simulation: From basic principles to applications_ , C. R. Phys. 19, 357 (2018).
  * [21] L. Tarruell and L. Sanchez-Palencia, _Quantum simulation of the Hubbard model with ultracold fermions in optical lattices_ , C. R. Phys. 19, 365 (2018).
  * [22] M. Aidelsburger, S. Nascimbene, and N. Goldman, _Artificial gauge fields in materials and engineered systems_ , C. R. Phys. 19, 394 (2018).
  * [23] J. Lebreuilly and I. Carusotto, _Quantum simulation of zero temperature quantum phases and incompressible states of light via non-Markovian reservoir engineering techniques_ , C. R. Phys. 19, 433 (2018).
  * [24] K. Le Hur, L. Henriet, L. Herviou, K. Plekhanov, A. Petrescu, T. Goren, M. Schiro, C. Mora, and P. P. Orth, _Driven dissipative dynamics and topology of quantum impurity systems_ , C. R. Phys. 19, 451 (2018).
  * [25] M. Bell, B. Douçot, M. Gershenson, L. Ioffe, and A. Petkovic, _Josephson ladders as a model system for 1d quantum phase transitions_ , C. R. Phys. 19, 484 (2018).
  * [26] F. Alet and N. Laflorencie, _Many-body localization: An introduction and selected topics_ , C. R. Phys. 19, 498 (2018).
  * [27] X.-L. Deng, D. Porras, and J. I. Cirac, _Effective spin quantum phases in systems of trapped ions_ , Phys. Rev. A 72, 063407 (2005).
  * [28] R. Islam, E. Edwards, K. Kim, S. Korenblit, C. Noh, H. Carmichael, G.-D. Lin, L.-M. Duan, C.-C. Joseph Wang, J. Freericks, _et al._ , _Onset of a quantum phase transition with a trapped ion quantum simulator_ , Nat. Comm. 2, 377 (2011).
  * [29] B. P. Lanyon, C. Hempel, D. Nigg, M. Müller, R. Gerritsma, F. Zähringer, P. Schindler, J. T. Barreiro, M. Rambach, G. Kirchmair, _et al._ , _Universal digital quantum simulation with trapped ions_ , Science 334(6052), 57 (2011).
  * [30] C. Schneider, D. Porras, and T. Schaetz, _Experimental quantum simulations of many-body physics with trapped ions_ , Rep. Prog. Phys. 75(2), 024401 (2012).
  * [31] A. Bermudez, L. Tagliacozzo, G. Sierra, and P. Richerme, _Long-range Heisenberg models in quasiperiodically driven crystals of trapped ions_ , Phys. Rev. B 95, 024431 (2017).
  * [32] V. Bendkowsky, B. Butscher, J. Nipper, J. P. Shaffer, R. Löw, and T. Pfau, _Observation of ultralong-range Rydberg molecules_ , Nature (London) 458, 1005 (2009).
  * [33] H. Weimer, M. Müller, I. Lesanovsky, P. Zoller, and H. P. Büchler, _A Rydberg quantum simulator_ , Nat. Phys. 6(5), 382 (May 2010).
  * [34] P. Schausz, M. Cheneau, M. Endres, T. Fukuhara, S. Hild, A. Omran, T. Pohl, C. Gross, S. Kuhr, and I. Bloch, _Observation of spatially ordered structures in a two-dimensional Rydberg gas_ , Nature (London) 491, 87 (2012).
  * [35] A. Browaeys, D. Barredo, and T. Lahaye, _Experimental investigations of dipole-dipole interactions between a few Rydberg atoms_ , J. Phys. B: At. Mol. Opt. Phys. 49(15), 152001 (2016).
  * [36] A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and T. Pfau, _Bose-Einstein condensation of chromium_ , Phys. Rev. Lett. 94, 160401 (Apr 2005).
  * [37] Q. Beaufils, R. Chicireanu, T. Zanon, B. Laburthe-Tolra, E. Maréchal, L. Vernac, J.-C. Keller, and O. Gorceix, _All-optical production of chromium Bose-Einstein condensates_ , Phys. Rev. A 77, 061601 (2008).
  * [38] M. Lu, N. Q. Burdick, S. H. Youn, and B. L. Lev, _Strongly dipolar Bose-Einstein condensate of dysprosium_ , Phys. Rev. Lett. 107, 190401 (2011).
  * [39] S. Baier, M. J. Mark, D. Petter, K. Aikawa, L. Chomaz, Z. Cai, M. Baranov, P. Zoller, and F. Ferlaino, _Extended Bose-Hubbard models with ultracold magnetic atoms_ , Science 352(6282), 201 (2016).
  * [40] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. Pfau, _The physics of dipolar bosonic quantum gases_ , Rep. Prog. Phys. 72(12), 126401 (2009).
  * [41] A. Micheli, G. K. Brennen, and P. Zoller, _A toolbox for lattice-spin models with polar molecules_ , Nat. Phys. 2(5), 341 (May 2006).
  * [42] B. Yan, S. A. Moses, B. Gadway, J. P. Covey, K. R. A. Hazzard, A. M. Rey, D. S. Jin, and J. Ye, _Observation of dipolar spin-exchange interactions with lattice-confined polar molecules_ , Nature (London) 501, 521 (2013).
  * [43] S. A. Moses, J. P. Covey, M. T. Miecnikowski, D. S. Jin, and J. Ye, _New frontiers for quantum gases of polar molecules_ , Nat. Phys. 13, 13 (2017).
  * [44] O. Firstenberg, T. Peyronel, Q.-Y. Liang, A. V. Gorshkov, M. D. Lukin, and V. Vuletic, _Attractive photons in a quantum nonlinear medium_ , Nat. Phys. 502, 71 (2013).
  * [45] L. Childress, M. V. Gurudev Dutt, J. M. Taylor, A. S. Zibrov, F. Jelezko, J. Wrachtrup, P. R. Hemmer, and M. D. Lukin, _Coherent dynamics of coupled electron and nuclear spin qubits in diamond_ , Science 314(5797), 281 (2006).
  * [46] G. Balasubramanian, P. Neumann, D. Twitchen, M. Markham, R. Kolesov, N. Mizuochi, J. Isoya, J. Achard, J. Beck, J. Tissler, _et al._ , _Ultralong spin coherence time in isotopically engineered diamond_ , Nat. Mater. 8, 383 (2009).
  * [47] F. Dolde, I. Jakobi, B. Naydenov, N. Zhao, S. Pezzagna, C. Trautmann, J. Meijer, P. Neumann, F. Jelezko, and J. Wrachtrup, _Room-temperature entanglement between single defect spins in diamond_ , Nat. Phys. 9, 139 (2013).
  * [48] E. H. Lieb and D. W. Robinson, _The finite group velocity of quantum spin systems_ , Comm. Math. Phys. 28, 251 (1972).
  * [49] M. B. Hastings and T. Koma, _Spectral gap and exponential decay of correlations_ , Comm. Math. Phys. 265, 781 (2006).
  * [50] P. Calabrese and J. Cardy, _Time dependence of correlation functions following a quantum quench_ , Phys. Rev. Lett. 96, 136801 (2006).
  * [51] M. Cheneau, P. Barmettler, D. Poletti, M. Endres, P. Schauss, T. Fukuhara, C. Gross, I. Bloch, C. Kollath, and S. Kuhr, _Light-cone-like spreading of correlations in a quantum many-body system_ , Nature (London) 481, 484 (2012).
  * [52] T. Fukuhara, P. Schausz, M. Endres, S. Hild, M. Cheneau, I. Bloch, and C. Gross, _Microscopic observation of magnon bound states and their dynamics_ , Nature (London) 502, 76 (2013).
  * [53] R. Geiger, T. Langen, I. E. Mazets, and J. Schmiedmayer, _Local relaxation and light-cone-like propagation of correlations in a trapped one-dimensional Bose gas_ , New J. Phys. 16, 053034 (2014).
  * [54] A. M. Läuchli and C. Kollath, _Spreading of correlations and entanglement after a quench in the one-dimensional Bose-Hubbard model_ , J. Stat. Mech. 2008(05), P05018 (2008).
  * [55] S. R. Manmana, S. Wessel, R. M. Noack, and A. Muramatsu, _Time evolution of correlations in strongly interacting fermions after a quantum quench_ , Phys. Rev. B 79, 155104 (2009).
  * [56] P. Barmettler, D. Poletti, M. Cheneau, and C. Kollath, _Propagation front of correlations in an interacting Bose gas_ , Phys. Rev. A 85, 053625 (2012).
  * [57] G. Carleo, F. Becca, L. Sanchez-Palencia, S. Sorella, and M. Fabrizio, _Light-cone effect and supersonic correlations in one- and two-dimensional bosonic superfluids_ , Phys. Rev. A 89, 031602(R) (2014).
  * [58] J. Despres, L. Villa, and L. Sanchez-Palencia, _Twofold correlation spreading in a strongly correlated lattice Bose gas_ , Sci. Rep. 9, 4135 (2019).
  * [59] L. Villa, J. Despres, and L. Sanchez-Palencia, _Unraveling the excitation spectrum of many-body systems from quantum quenches_ , Phys. Rev. A 100, 063632 (Dec 2019).
  * [60] L. Villa, J. Despres, S. J. Thomson, and L. Sanchez-Palencia, _Local quench spectroscopy of many-body quantum systems_ , Phys. Rev. A 102, 033337 (Sep 2020).
  * [61] P. Hauke and L. Tagliacozzo, _Spread of correlations in long-range interacting quantum systems_ , Phys. Rev. Lett. 111, 207202 (2013).
  * [62] J. Eisert, M. van den Worm, S. R. Manmana, and M. Kastner, _Breakdown of quasilocality in long-range quantum lattice models_ , Phys. Rev. Lett. 111, 260401 (2013).
  * [63] M. Foss-Feig, Z.-X. Gong, C. W. Clark, and A. V. Gorshkov, _Nearly linear light cones in long-range interacting quantum systems_ , Phys. Rev. Lett. 114, 157201 (2015).
  * [64] D. V. Else, F. Machado, C. Nayak, and N. Y. Yao, _Improved Lieb-Robinson bound for many-body Hamiltonians with power-law interactions_ , Phys. Rev. A 101, 022333 (2020).
  * [65] L. Cevolani, G. Carleo, and L. Sanchez-Palencia, _Spreading of correlations in exactly solvable quantum models with long-range interactions in arbitrary dimensions_ , New J. Phys. 18, 093002 (2016).
  * [66] C.-F. Chen and A. Lucas, _Finite speed of quantum scrambling with long range interactions_ , Phys. Rev. Lett. 123, 250605 (Dec 2019).
  * [67] M. C. Tran, C.-F. Chen, A. Ehrenberg, A. Y. Guo, A. Deshpande, Y. Hong, Z.-X. Gong, A. V. Gorshkov, and A. Lucas, _Hierarchy of linear light cones with long-range interactions_ , Phys. Rev. X 10, 031009 (Jul 2020).
  * [68] T. Kuwahara and K. Saito, _Strictly linear light cones in long-range interacting systems of arbitrary dimensions_ , Phys. Rev. X 10, 031010 (Jul 2020).
  * [69] L. Cevolani, G. Carleo, and L. Sanchez-Palencia, _Protected quasi-locality in quantum systems with long-range interactions_ , Phys. Rev. A 92, 041603(R) (2015).
  * [70] J. Schachenmayer, A. Pikovski, and A. M. Rey, _Dynamics of correlations in two-dimensional quantum spin models with long-range interactions: A phase-space Monte-Carlo study_ , New J. Phys. 17(6), 065009 (2015).
  * [71] A. S. Buyskikh, M. Fagotti, J. Schachenmayer, F. Essler, and A. J. Daley, _Entanglement growth and correlation spreading with variable-range interactions in spin and fermionic tunneling models_ , Phys. Rev. A 93, 053620 (2016).
  * [72] L. Cevolani, J. Despres, G. Carleo, L. Tagliacozzo, and L. Sanchez-Palencia, _Universal scaling laws for correlation spreading in quantum systems with short- and long-range interactions_ , Phys. Rev. B 98, 024302 (2018).
  * [73] P. Richerme, Z.-X. Gong, A. Lee, C. Senko, J. Smith, M. Foss-Feig, S. Michalakis, A. V. Gorshkov, and C. Monroe, _Non-local propagation of correlations in quantum systems with long-range interactions_ , Nature (London) 511, 198 (2014).
  * [74] P. Jurcevic, B. P. Lanyon, P. Hauke, C. Hempel, P. Zoller, R. Blatt, and C. F. Roos, _Quasiparticle engineering and entanglement propagation in a quantum many-body system_ , Nature (London) 511, 202 (2014).
  * [75] D. Jaksch, J. I. Cirac, P. Zoller, S. L. Rolston, R. Côté, and M. D. Lukin, _Fast quantum gates for neutral atoms_ , Phys. Rev. Lett. 85, 2208 (Sep 2000).
  * [76] M. D. Lukin, M. Fleischhauer, R. Cote, L. M. Duan, D. Jaksch, J. I. Cirac, and P. Zoller, _Dipole blockade and quantum information processing in mesoscopic atomic ensembles_ , Phys. Rev. Lett. 87, 037901 (Jun 2001).
  * [77] M. F. Maghrebi, Z.-X. Gong, M. Foss-Feig, and A. V. Gorshkov, _Causality and quantum criticality in long-range lattice models_ , Phys. Rev. B 93, 125128 (2016).
  * [78] J. Haegeman, B. Pirvu, D. J. Weir, J. I. Cirac, T. J. Osborne, H. Verschelde, and F. Verstraete, _Variational matrix product ansatz for dispersion relations_ , Phys. Rev. B 85, 100408 (2012).
  * [79] For details, see Supplemental Material. It discusses the $\epsilon$ method for determining the correlation, spin, and entanglement edges, the spreading of spin correlations along the $x$ direction as well as strong quenches, and the LSWT for the spreading of the local magnetization and entanglement entropies. It includes Refs. [90, 91].
  * [80] The von Neumann entropy also reads as $\mathcal{S}_{\textrm{vN}}=\mathcal{S}_{n\rightarrow 1}(R,t)=-\mathrm{Tr}\left\\{\mathaccentV{hat}05E{\rho}(R,t)\mathrm{log}\left[\mathaccentV{hat}05E{\rho}(R,t)\right]\right\\}$.
  * [81] V. Alba and P. Calabrese, _Entanglement and thermodynamics after a quantum quench in integrable systems_ , Proc Natl Acad Sci USA 114(30), 7947 (Jul. 2017).
  * [82] A. Bastianello and P. Calabrese, _Spreading of entanglement and correlations after a quench with intertwined quasiparticles_ , SciPost Phys. 5(4), 033 (Oct. 2018).
  * [83] P. Calabrese, (Aug. 2020).
  * [84] A. Hurwitz, _Einige Eigenschaften der Dirichletschen Functionen $F(s)=\sum(\frac{D}{n})\cdot\frac{1}{n}$, die bei der Bestimmung der Klassenanzahlen binärer quadratischer Formen auftreten_ (Springer Basel, Basel, 1932), pp. 72–88, ISBN 978-3-0348-4161-0.
  * [85] W. Magnus, F. Oberhettinger, and R. P. Soni, _The gamma function and related functions_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 1966), pp. 1–37, ISBN 978-3-662-11761-3.
  * [86] J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong, and C. Monroe, _Observation of a many-body dynamical phase transition with a 53-qubit quantum simulator_ , Nature (London) 551, 601 (2017).
  * [87] N. M. Linke, S. Johri, C. Figgatt, K. A. Landsman, A. Y. Matsuura, and C. Monroe, _Measuring the Rényi entropy of a two-site Fermi-Hubbard model on a trapped ion quantum computer_ , Phys. Rev. A 98, 052334 (2018).
  * [88] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch, C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F. Roos, _Probing Rényi entanglement entropy via randomized measurements_ , Science 364, 260 (2019).
  * [89] M. Fishman, S. R. White, and E. M. Stoudenmire, (2020), eprint 2007.14822.
  * [90] T. Holstein and H. Primakoff, _Field dependence of the intrinsic domain magnetization of a ferromagnet_ , Phys. Rev. 58, 1098 (1940).
  * [91] N. N. Bogolyubov, J. Phys. USSR 11, 23 (1947).
  * [92] Note that a different value of $h/J$ is used in the two figures and the values of $V_{\textrm{CE}}$ in Fig. 3(c) and $V_{\textrm{SE}}$ in Fig. 4(c) in the main paper are not related.

Supplemental Material for
Spreading of Correlations and Entanglement in the Long-Range Transverse Ising
Chain

In this Supplemental Material, we provide additional technical details of the
methodology and analysis underlying the central results discussed in the main
paper. In Sec. S1, we detail the fitting procedure to extract the scaling
exponent of the spin-correlation, local magnetization, and entanglement edges.
In Sec. S2, we discuss the spreading of spin correlations in the $x$ direction
and in Sec. S3 results for large quenches. Sections S4 and S5 outline the LSWT
theory for the spreading of the local magnetization $\langle
S^{z}_{R}(t)\rangle$ and entanglement entropies $\mathcal{S}_{n}(R,t)$,
respectively, after a local quench. Finally, Sec. S6 gives additional TDVP and
LSWT results for the spreading of entanglement considering various Rényi
entropies.

##  S1 Determination of the edge

For all the data reported in the main paper, we have determined the spin-
correlation, local magnetization, and entanglement edges by tracking the
ensemble of points in the $R$–$t$ plane where the signal reaches a fraction
$\epsilon$ of its maximal value. Since this threshold line depends on the
value of $\epsilon$, we have systematically scanned $\epsilon$, e.g. from
$0.01$ to $0.12$ for $G_{z}(R,t)$ and from $0.20$ to $0.80$ for
$\mathcal{S}_{n}(R,t)$. In all the cases considered in this work, we find that
the edge is well fitted by the algebraic law

$t=a\times R^{\beta}.$ (S1)

While the coefficient $a$ fundamentally depends on $\epsilon$, the scaling law
exponent $\beta$ is nearly independent of $\epsilon$. An example of such an
analysis, plotted in log-log scale, is shown in Fig. S1 where the various
lines correspond to different values of $\epsilon$. The fact that they are
parallel straight lines validates the scaling law $t\sim\times
R^{\beta_{\textrm{\tiny CE}}}$ with an exponent $\beta_{\textrm{\tiny CE}}$
nearly independent of $\epsilon$.

Figure S1:  Spreading of the $G_{z}$ spin correlation function in the quasi-
local regime ($\alpha=1.7<2$), for a global quench from
$(h/J)_{\mathrm{i}}=50$ to $(h/J)_{\mathrm{f}}=1$, both in the $z$ polarized
phase [same data as in the Fig. 2(a) of the main paper]. The squares, disks,
and diamonds indicate points where $G_{z}(R,t)$ reaches a fraction $\epsilon$
of its maximum value for various values of $\epsilon$. The corresponding lines
are linear fits to these points in log-log scale, consistently with Eq. (S1).
The brown squares correspond to $\epsilon=2.6\%$ for which we find
$\beta_{\textrm{\tiny CE}}\simeq 1.27$, the green disks to $\epsilon=6.5\%$
for which we find $\beta_{\textrm{\tiny CE}}\simeq 1.23$, and the white
diamonds to $\epsilon=11.2\%$ for which we find $\beta_{\textrm{\tiny
CE}}\simeq 1.22$.

##  S2 Spreading of the $\mathbf{G_{x}}$ spin correlation function

In this section, we discuss analytic and numerical results for the spreading
of spin correlations along the $x$ direction,

$G_{x}(R,t)=G_{x}^{0}(R,t)-G_{x}^{0}(R,0)\qquad\textrm{with}\qquad
G_{x}^{0}(R,t)=\langle\hat{S}^{x}_{R}(t)\hat{S}^{x}_{0}(t)\rangle-\langle\hat{S}^{x}_{R}(t)\rangle\langle\hat{S}^{x}_{0}(t)\rangle,$
(S2)

where $\hat{S}^{x}_{R}(t)$ is the spin operator along $x$ at position $R$ and
time $t$.

###  S2.1 Linear spin wave theory

In the $z$-polarized phase, the LRTI model can be diagonalized using a
Holstein-Primakoff transformation [90],

$\hat{S}^{x}_{R}\simeq\frac{1}{2}\left(\hat{a}^{{\dagger}}_{R}+\hat{a}_{R}\right)\qquad,\qquad\hat{S}^{y}_{R}\simeq\frac{-1}{2i}\left(\hat{a}^{{\dagger}}_{R}-\hat{a}_{R}\right)\qquad,\quad\textrm{and}\qquad\hat{S}^{z}_{R}\simeq\frac{1}{2}-\hat{a}^{{\dagger}}_{R}\hat{a}_{R},$
(S3)

where $\hat{a}_{R}(t)$ is the annihilation operator of a boson at site $R$ and
time $t$. Inserting these formulas into the LRTI Hamiltonian [Eq. (1) of the
main paper], we find a quadratic Bose Hamiltonian, see for instance Refs. [61,
65] for details. The latter is readily diagonalized using a standard
Bogoliubov transformation [91]. The elementary excitations are magnons (spin-
wave excitations), characterized by the gapped spectrum

$E_{k}=2\sqrt{h[h+JP_{\alpha}(k)]}=\sqrt{\mathcal{A}_{k}^{2}-\mathcal{B}_{k}^{2}}\qquad,\quad\mathrm{with}\qquad\mathcal{A}_{k}=JP_{\alpha}(k)+2h\qquad\mathrm{and}\qquad\mathcal{B}_{k}=JP_{\alpha}(k),$
(S4)

where $P_{\alpha}(k)=\sum_{R}e^{ikR}/R^{\alpha}$ denotes the Fourier transform
of the long-range term in the Hamiltonian.

For a quench on the exchange amplitude $J$ from $J_{\mathrm{i}}$ to
$J_{\mathrm{f}}$, the $G_{x}(R,t)$ correlation function can be cast in the
form [72]

$G_{x}(R,t)\simeq
g(R)-\int_{-\pi}^{\pi}\frac{\mathrm{d}k}{2\pi}\mathcal{F}(k)\left\\{\frac{e^{i(kR+2E_{k}^{\mathrm{f}}t)}+e^{i(kR-2E_{k}^{\mathrm{f}}t)}}{2}\right\\}\quad,\quad\mathrm{with}\quad\mathcal{F}(k)=\frac{h(J_{\mathrm{i}}-J_{\mathrm{f}})P_{\alpha}(k)}{8[h+J_{\mathrm{f}}P_{\alpha}(k)]\sqrt{h[h+J_{\mathrm{i}}P_{\alpha}(k)]}},$
(S5)

where the index "i" refers to the pre-quench (initial) Hamiltonian and the
index "f" to the post-quench (final) Hamiltonian. The function
$\mathcal{F}(k)$ gives the weight of each quasi-particle according to their
mode $k$. It depends on the observable and the quench.

To characterize the asymptotic behavior, $R,t\rightarrow+\infty$ along a
constant line $R/t$, one can then rely on the stationary-phase approximation,
applied to Eq. (S5). It yields

$G_{x}(R,t)\sim\frac{\mathcal{F}(k_{\mathrm{sp}})}{(|\partial^{2}_{k}E_{k_{\mathrm{sp}}}|t)^{1/2}}\cos(k_{\mathrm{sp}}R-2E_{k_{\mathrm{sp}}}t+\phi),$
(S6)

where $\phi$ is a constant phase irrelevant to our study, $k_{\mathrm{sp}}$ is
the stationary-phase quasi-momentum, given by the solution of the equation
$2V_{\mathrm{g}}(k_{\mathrm{sp}})=2\partial_{k}E_{k_{\mathrm{sp}}}=R/t>0$ with
$V_{\textrm{g}}$ the group velocity associated to the spin-wave excitations.

#### Quasi-local regime –

The quasi-local regime corresponds to the case where the quasi-particle energy
$E_{k}$ is finite over the whole first Brillouin zone, $k\in[-\pi,+\pi]$, but
the group velocity $V_{\textrm{g}}(k)$ presents a divergence. The space-time
behavior of $G_{x}$ in the vicinity of the CE (correlation edge) is dominated
by the quasi-particles propagating with the highest group velocities. Due to
the infrared divergence, only the limit behavior in $k\rightarrow 0$ is
relevant. There the long-range term reads as $P_{\alpha}(k)\approx
P_{\alpha}(0)+P_{\alpha}^{\prime}|k|^{\alpha-1}$, with $P_{\alpha}(0)>0$,
$P^{\prime}_{\alpha}<0$, and the excitation energy is
$E_{k}\simeq\Delta-c|k|^{z}$, where $\Delta=2\sqrt{h[h+JP_{\alpha}(0)]}>0$ is
the gap, $z=\alpha-1\geq 0$ is the dynamical exponent, and
$c=\sqrt{h/[h+JP_{\alpha}(0)]}J|P_{\alpha}^{\prime}|>0$ is a prefactor. Then,
computing $k_{\mathrm{sp}}$ and injecting it into Eq. (S6), we find

$G_{x}(R,t)\sim\frac{t^{\gamma}}{R^{\chi}}\cos[A_{z}\Big{(}\frac{t}{R^{z}}\Big{)}^{\frac{1}{1-z}}-2\Delta
t+\phi\Big{]},$ (S7)

with

$A_{z}=2c(2cz)^{z/(1-z)}(1-z)\qquad\mathrm{,}\qquad\gamma=\frac{\nu+1/2}{1-z}\qquad\mathrm{,}\quad\mathrm{and}\qquad\chi=\frac{\nu+(2-z)/2}{1-z},$
(S8)

where $\nu\geq 0$ is the scaling exponent of the amplitude function in the
infrared limit, $\mathcal{F}(k)\sim|k|^{\nu}$. It follows from Eq. (S5) and
the approximation of $P_{\alpha}$ in the infrared limit that $\nu=0$ for the
$G_{x}$ spin correlation function. On the one hand, the CE is found by
imposing the condition that the prefactor is constant. The latter leads to the
algebraic form $t\propto R^{\beta_{\mathrm{CE}}}$ with $\beta_{\textrm{\tiny
CE}}=\chi/\gamma=3-\alpha$. Since $1\leq\alpha<2$ in the quasi-local, the CE
is always sub-ballistic, i.e. $\beta_{\mathrm{CE}}>1$. On the other hand, the
spreading law of the local extrema is determined by the equation

$A_{z}\Big{(}\frac{t}{R^{z}}\Big{)}^{\frac{1}{1-z}}-2\Delta
t+\phi=\mathrm{cst}\qquad\mathrm{leading\leavevmode\nobreak\ to}\qquad
A_{z}\Big{(}\frac{t}{R}\Big{)}^{\frac{z}{1-z}}-2\Delta\rightarrow 0.$ (S9)

The maxima are thus ballistic, i.e. $t\sim R^{\beta}_{\textrm{\tiny m}}$ with
$\beta_{\textrm{\tiny m}}=1$.

#### Local regime –

In the local regime, corresponding to the case where both the quasi-particle
energy $E_{k}$ and the group velocity $V_{\textrm{g}}(k)$ are finite over the
whole first Brillouin zone, there exists a quasi-momentum $k^{*}$ such that
the group velocity is maximum,
$V_{\textrm{g}}(k^{*})=\mathrm{max}_{k}[V_{\textrm{g}}(k)]$. Hence, the
stationary-phase condition $2V_{\textrm{g}}(k_{\mathrm{sp}})=R/t$ has a
solution only for $R/t\leq 2V_{\textrm{g}}(k^{*})$. The CE is determined by
the spreading of the quasi-particles with a quasi-momentum
$k_{\mathrm{sp}}\simeq k^{*}$. It is thus ballistic, with the CE velocity
$V_{\textrm{\tiny CE}}=2V_{\textrm{g}}(k^{*})$. Moreover, in the vicinity of
the CE, the motion of the local maxima is determined by the phase factor in
Eq. (S6) with $k_{\mathrm{sp}}=k^{*}$. It follows that the local extrema
propagate ballistically at the velocity $V_{\textrm{m}}=2V_{\varphi}(k^{*})$
where $V_{\varphi}(k)=E_{k}/k$ is the quasi-particle phase velocity. Note that
here, $k^{*}<0$ and $V_{\textrm{m}}=2V_{\varphi}(k^{*})=2E_{k^{*}}/k^{*}<0$.

###  S2.2 Numerical TDVP results

We have performed numerical TDVP calculations for the spreading of the spin
correlation function along $x$ for the same parameters as for its counterpart
in the $z$ direction, discussed in the text.

The left panel of Fig. S2 shows a typical result, plotted in log-log scale, in
the quasi-local regime. The underlying checkerboard-like structure is
characteristic of the $G_{x}$ function and was similarly found using meanfield
calculations in Ref. [72]. In spite of this complex structure, we can identify
a CE, using the same $\epsilon$ method used for the $G_{z}$ correlation
function, as well as propagating local maxima. The green solid and dashed blue
lines correspond to alegbraic fits to the CE and the local extrema,
respectively. The associated scaling law exponents $\beta_{\textrm{\tiny CE}}$
and $\beta_{\textrm{\tiny m}}$ are in very good agreement with the theoretical
ones discussed in Sec. S2.1, see filled symbols in Fig. 2(b) of the main
paper.

The right panel of Fig. S2 shows a typical result in the local regime, now
plotted in linear scale. The solid green and dashed blue lines are linear fits
to the CE and local extrema, respectively. The corresponding numerical
velocities $V_{\textrm{\tiny CE}}$ and $V_{\textrm{m}}$ are in good agreement
with the theoretical ones, see filled symbols in Fig. 2(c) of the main paper.

Figure S2:  Spreading of the $G_{x}$ spin correlations. Left panel (log-log
scale): Quasi-local regime with $\alpha=1.7<2$ for a global quench from
$(h/J)_{\mathrm{i}}=50$ to $(h/J)_{\mathrm{f}}=1$, both in the $z$ polarized
phase of the LRTI model. Right panel (linear scale): Local regime with
$\alpha=3<2$ for a global quench starting from $(h/J)_{\mathrm{i}}=0.9$ to
$(h/J)_{\mathrm{f}}=1$. The solid green and dashed blue lines are algebraic
(for $\alpha<2$) or linear (for $\alpha>2$) fits to the CE and the extrema,
respectively.

##  S3 Large quenches

In this section, we present numerical results for large quenches using the
same TDVP calculations as for the results of the main paper.

Figure S3 shows results for both the $G_{x}(R,t)$ (upper row) and the
$G_{z}(R,t)$ (lower panel) correlation functions in the quasi-local regime,
$\alpha=1.7$. Form left to right, we quench from decreasing values of the
initial parameter $(h/J)_{\textrm{i}}$, ranging from the $z$-polarized phase
to the $x$-Néel phase. The post-quench parameter $(h/J)_{\textrm{f}}$ is fixed
to a value corresponding to the $z$-polarized phase. When changing the initial
value $(h/J)_{\textrm{i}}$ towards the Néel phase, we find that an additional
staggered structure, reminiscent of the Néel ordering, appears. For the
$G_{x}(R,t)$ correlation function, it rapidly dominates the signal and it
becomes difficult to identify propagating local maxima and a correlation edge.
In contrast, the signal is clearer for the $G_{z}(R,t)$ correlation function.
While a checkerboard-like structure appears, it is still possible to identify
the local maxima and the correlation edge. We find that both are algebraic
with dynamical exponents $\beta_{\textrm{\tiny m}}$ and $\beta_{\textrm{\tiny
CE}}$ that do not significantly change with the strength of the quench.

Figure S3:  Spreading of the $G_{x}$ (upper row) and $G_{z}$ (lower row) spin
correlations in the quasi-local regime, $\alpha=1.7$ (log-log scale). The
quench is performed from an increasing value of the initial parameter
$(h/J)_{\mathrm{i}}$, ranging from the $z$ polarized phase to the $x$-Néel
phase, to the fixed value $(h/J)_{\mathrm{f}}=1$, corresponding to the $z$
polarized phase. (a) $\hat{H}_{\mathrm{i}}$ deep in the $z$ polarized phase,
$(h/J)_{\mathrm{i}}=0.8$; (b) $\hat{H}_{\mathrm{i}}$ in the $z$ polarized
phase close to the critical point, $(h/J)_{\mathrm{i}}=0.4$; (c)
$\hat{H}_{\mathrm{i}}$ in the $x$ Néel phase, $(h/J)_{\mathrm{i}}=0.25$; (d)
$\hat{H}_{\mathrm{i}}$ deep in the $x$ Néel phase, $(h/J)_{\mathrm{i}}=0.1$.
The dashed blue and solid green lines are algebraic fits to propagating local
extrema and to the CE respectively.

Figure S4 shows the same calculations, now performed in the local regime,
$\alpha=3$. Similar conclusions as for the quasi-local regime hold. The
$G_{x}(R,t)$ correlation function is rapidly blurred. In contrast, the
$G_{z}(R,t)$ is more robust and allows us to fit the dynamical exponents
$\beta_{\textrm{\tiny m}}$ and $\beta_{\textrm{\tiny CE}}$, yielding values
nearly independent of $(h/J)_{\textrm{i}}$. In fact, for the $G_{z}(R,t)$
correlation function, we find no significant impact of increasing the strength
of the quench.

---
Figure S4: Spreading of the $G_{x}$ (upper row) and $G_{z}$ (lower row) spin
correlations in the local regime, $\alpha=3$ (linear scale). The quench is
performed from an increasing value of the initial parameter
$(h/J)_{\mathrm{i}}$, ranging from the $z$ polarized phase to the $x$-Néel
phase, to the fixed value $(h/J)_{\mathrm{f}}=1$, corresponding to the
$z$-polarized phase. (a) $\hat{H}_{\mathrm{i}}$ deep in the $z$-polarized
phase, $(h/J)_{\mathrm{i}}=0.9$; (b) $\hat{H}_{\mathrm{i}}$ in the $z$
polarized phase close to the critical point, $(h/J)_{\mathrm{i}}=0.59$; (c)
$\hat{H}_{\mathrm{i}}$ in the $x$ Néel phase, $(h/J)_{\mathrm{i}}=0.25$; (d)
$\hat{H}_{\mathrm{i}}$ deep in the $x$ Néel phase, $(h/J)_{\mathrm{i}}=0.1$.
The dashed blue and solid green lines are linear fits to propagating local
extrema and to the CE respectively.

##  S4 Scaling laws for $\mathbf{\langle S^{z}_{R}(t)\rangle}$ after a local
quench

The local magnetization $\langle\hat{S}^{z}_{R}(t)\rangle$ can be computed
analytically for a local quench in the $z$-polarized phase where $h\gg J$,
using an approach similar to that presented in Sec. S2.1. Starting from the
initial state $\ket{\Psi_{0}}\propto\hat{S}^{-}_{0}\ket{\Psi_{\textrm{\tiny
GS}}}$ obtained by lowering the central ($R=0$) spin of the ground state
(which is close to the product state
$\ket{\uparrow...\uparrow\downarrow\uparrow...\uparrow}$), and applying both
the Holstein-Primakoff transformation and the bosonic Bogoliubov approach, we
find

$\displaystyle
1/2-\langle\hat{S}^{z}_{R}(t)\rangle\simeq\left|\int_{-\pi}^{\pi}\frac{\mathrm{d}k}{2\pi}\mathcal{F}_{1}(k)\frac{e^{i(kR+E_{k}t)}+e^{-i(kR-
E_{k}t)}}{2}\right|^{2}+\left|\int_{-\pi}^{\pi}\frac{\mathrm{d}k}{2\pi}\mathcal{F}_{2}(k)\frac{e^{i(kR+E_{k}t)}+e^{-i(kR-
E_{k}t)}}{2}\right|^{2}.$ (S10)

The quantities $\mathcal{F}_{1}(k)$ and $\mathcal{F}_{2}(k)$ are two different
quasi-momentum-dependent amplitude functions, which read as

$\mathcal{F}_{1}(k)=\frac{1}{2}\left(\frac{\mathcal{A}_{k}}{E_{k}}+1\right)=\frac{1}{2}\left(\frac{2h+JP_{\alpha}(k)}{2\sqrt{h[h+JP_{\alpha}(k)]}}+1\right)\qquad\textrm{and}\qquad\mathcal{F}_{2}(k)=-\frac{\mathcal{B}_{k}}{2E_{k}}=-\frac{JP_{\alpha}(k)}{4\sqrt{h[h+JP_{\alpha}(k)]}}.$
(S11)

For large transverse fields, $h\gg J$, we find $\mathcal{F}_{1}(k)\simeq
1\gg\mathcal{F}_{2}(k)\sim J/h$ and the second term in Eq. (S10) may be
neglected.

#### Quasi-local regime –

In the quasi-local regime, owing to the infrared divergence of the quasi-
particle velocities, the spreading of the local magnetization (as well as that
of correlations) is dominated by the Fourier components with $k\rightarrow 0$,
see Eq. (S10). There, we may write $E_{k}\simeq\Delta-c|k|^{z}$, where
$\Delta=2\sqrt{h[h+JP_{\alpha}(0)]}$ is the gap and $z=\alpha-1$ is the
dynamical exponent, see also Sec. S2.1. Inserting this expression into Eq.
(S10), we find that the gap $\Delta$ may be factorized, taken out of the
integrals, and then vanishes owing to the square moduli. It explains the
irrelevance of the gap for the spreading of the local magnetization in the
quasi-local regime. Applying the stationary phase approximation to evaluate
the integrals along the lines of Ref. [72], we then find that the spreading of
$\langle\hat{S}^{z}_{R}(t)\rangle$ is super-ballistic, with
$\beta_{\textrm{\tiny m}}=1-z$. This is consistent with the TDVP results shown
in the main paper, see Fig. 3(b) of the main paper.

#### Local regime –

In the local regime, both the quasi-particle energy $E_{k}$ and group velocity
$V_{\textrm{g}}(k)$ are finite over the whole Brillouin zone. The spin edge is
thus determined by the maximum group velocity, $Vg(k^{*})$. The calculation is
similar to that outlined for the spreading of the spin correlation in Sec.
S2.1. The only relevant difference is that the phase factor for the spreading
of the spin correlations has an extra factor of $2$, which is absent for the
local magnetization, see Eq. (S5) versus Eq. (S10). As a result, the spin edge
propagates at the velocity $V_{\textrm{\tiny SE}}=V_{\textrm{g}}(k^{*})$,
while the correlation edge propagates at the velocity $V_{\textrm{\tiny
CE}}=2V_{\textrm{g}}(k^{*})$. This is consistent with the TDVP results, see
Fig. 3(c) of the main paper 222Note that a different value of $h/J$ is used in
the two figures and the values of $V_{\textrm{\tiny CE}}$ in Fig. 3(c) and
$V_{\textrm{\tiny SE}}$ in Fig. 4(c) in the main paper are not related..

More precisely, applying the stationary phase approximation, we find

$\displaystyle
1/2-\langle\hat{S}^{z}_{R}(t)\rangle\sim\left|\frac{\mathcal{F}_{1}(k^{*})}{(|\partial^{2}_{k}E_{k^{*}}|t)^{1/2}}[\cos(k^{*}R-E_{k^{*}}t+\phi)-i\sin(k^{*}R-E_{k^{*}}t+\phi)]\right|^{2}$
(S12)

where $\phi=-(\pi/4)\mathrm{sgn}(\partial_{k}^{2}E_{k_{\mathrm{sp}}}t)$ is a
constant phase term. The local magnetization is thus the sum of two terms,
associated to the $\cos$ and $\sin$ terms. Each has local maxima that
propagate at the velocity $V_{\varphi}(k^{*})$. The two terms are, however, in
phase opposition, which mutually cancel the local maxima. We hence expect a
single structure in the vicinity of the spin edge, consistently with the TDVP
results. A similar effect has been found for density correlations deep in the
Mott insulator phase of the Bose-Hubbard model [58]. This effect is also
expected in the quasi-local regime but the extinction of the local maxima
happens on very large distances and in the TDVP we still observe them up to
the largest accessible times and distances, see Fig. 2 of the main paper.

##  S5 Density matrix and entanglement entropies after a local quench

In this section we discuss analytic and numerical results for the entanglement
entropy following a local quench, using linear spin wave theory (LSWT).

###  S5.1 Spin wave analysis

The initial state prepared by the local quench is obtained from the ground
state (GS) by an up-down spin flip at the central lattice site. The time
dependence of the system’s state is given by

$\ket{\Psi_{0}(t)}=\tilde{\mathcal{N}}\hat{S}^{-}_{0}(t)\ket{\Psi_{\mathrm{GS}}}=\tilde{\mathcal{N}}\hat{a}^{\dagger}_{0}(t)\ket{0_{b}}=\tilde{\mathcal{N}}\sum_{k=-\pi}^{\pi}u_{k}\mathrm{e}^{\mathrm{i}E_{k}t}\,\hat{b}_{k}^{\dagger}(t=0)\ket{0_{b}}\,,$
(S13)

where $\hat{S}^{-}_{0}=S^{x}_{0}-iS^{y}_{0}$ is the spin lowering operator at
the central site (labelled by $j=0$), $\hat{a}_{0}^{\dagger}$ is the bosonic
Holstein-Primakoff creation operator [see Eq. (S3)], $\hat{b}_{k}(0)$ is the
Bogoliubov quasi-particle annihilation operator [91],
$u_{k}=\mathrm{sign}(\mathcal{A}_{k})\sqrt{\frac{1}{2}\left(\frac{\absolutevalue{\mathcal{A}_{k}}}{E_{k}}+1\right)}$,
$\ket{0_{b}}$ is the quasi-particle vacuum, and
$\tilde{\mathcal{N}}=\left(\sum_{k}u_{k}^{2}\right)^{-1/2}$ stands for the
normalization of the state. Crucially, for $h/J\gg 1$, the initial state is a
superposition of single-particle excitations, which allows us to proceed
analytically. Defining $S$ as the entire system and $A$ and $B$ as
complementary subsystems, we may decompose the initial state into the cases
where the excitation is in either subsystem $A$ or $B$ as
$\sum_{m\in{S}}\ket{1_{m}}=\sum_{m\in
A}\ket{1_{m}}\otimes\ket{0_{B}}+\sum_{m\in B}\ket{0_{A}}\otimes\ket{1_{m}}$,
where $\ket{0_{{A},{B}}}$ is the unique vacuum of the subsystems $A$ and $B$,
respectively. The state of the entire system thus reads as

$\displaystyle\ket{\Psi_{0}(t)}$ $\displaystyle=\tilde{\mathcal{N}}\sum_{m\in
B}\sum_{k=-\pi}^{\pi}u_{k}\mathrm{e}^{\mathrm{i}\left(E_{k}t+mk\right)}\ket{0_{A}}\otimes\ket{1_{m}}+\tilde{\mathcal{N}}\sum_{m\in
A}\sum_{k=-\pi}^{\pi}u_{k}\mathrm{e}^{\mathrm{i}\left(E_{k}t+mk\right)}\ket{1_{m}}\otimes\ket{0_{B}}$
(S14)
$\displaystyle\coloneqq\sqrt{\lambda_{1}}\mathrm{e}^{\mathrm{i}\theta_{1}}\ket{0_{A}}\otimes\ket{\chi_{B}}+\sqrt{\lambda_{2}}\mathrm{e}^{\mathrm{i}\theta_{2}}\ket{\chi_{A}}\otimes\ket{0_{B}}\,,$
(S15)

where $\lambda_{1}$ and $\lambda_{2}$ are real-valued numbers such that
$\lambda_{1}+\lambda_{2}=1$, and $\theta_{1}$ and $\theta_{2}$ are some
phases. The density matrix of the entire system thus contains four terms,

$\displaystyle\rho(t)$
$\displaystyle=\lambda_{1}\outerproduct{0_{A}}{0_{A}}\otimes\outerproduct{\chi_{B}}{\chi_{B}}+\lambda_{2}\outerproduct{\chi_{A}}{\chi_{A}}\otimes\outerproduct{0_{B}}{0_{B}}$
$\displaystyle\quad+\sqrt{\lambda_{1}\lambda_{2}}\mathrm{e}^{i(\theta_{1}-\theta_{2})}\outerproduct{0_{A}}{\chi_{A}}\otimes\outerproduct{\chi_{B}}{0_{B}}+\sqrt{\lambda_{1}\lambda_{2}}\mathrm{e}^{i(\theta_{2}-\theta_{1})}\outerproduct{\chi_{A}}{0_{A}}\otimes\outerproduct{0_{B}}{\chi_{B}}\,,$
(S16)

while the reduced density matrix of subsystem $A$ contains two terms,

$\displaystyle\rho_{A}(t)$
$\displaystyle=\tr_{B}(\rho(t))=\lambda_{1}(t)\outerproduct{0_{A}}{0_{A}}+\lambda_{2}(t)\outerproduct{\chi_{A}}{\chi_{A}}\,.$
(S17)

In the limit of $h/J\gg 1$ the Bogoliubov coefficient approaches unity,
$u_{k}\simeq 1$. Since the two eigenvalues are constrained to add to unity and
have the same functional form, we may without loss of generality restrict our
analysis to one of them. The $t\rightarrow\infty\,,R\rightarrow\infty$
behavior of the second eigenvalue of the reduced density matrix in the
infinite volume limit using the stationary phase approximation is given by

$\displaystyle\lambda_{2}$
$\displaystyle\simeq{\tilde{\mathcal{N}}^{2}}\sum_{m\in
A}\absolutevalue{\sum_{k=-\pi}^{\pi}\mathrm{e}^{\mathrm{i}\left(E_{k}t+mk\right)}}^{2}\overset{N\rightarrow\infty}{\propto}\tilde{\mathcal{N}}^{2}\,\sum_{m\in
A}\frac{1}{t\absolutevalue{\partial^{2}_{k}E_{k_{\mathrm{sp}}}}}\cos^{2}\\!\big{[}E(k_{\mathrm{sp}})t-k_{\mathrm{sp}}m+\phi\big{]}$
(S18)

where $k_{sp}$ is defined by the stationary phase condition
$V_{\textrm{g}}(k_{\mathrm{sp}})=m/t$ with $V_{\textrm{g}}(k)=\partial
E_{k}/\partial k$.

###  S5.2 Entanglement Entropies

In order to compute the entanglement entropies from the reduced density
matrix, we now consider the partitions $A=\left\\{R,R+1,\ldots,N/2\right\\}$
and $B=\left\\{-N/2,-N/2+1,\ldots,R-1\right\\}$ and look for an analytic
expression for its eigenvalues.

#### Quasi-local regime –

The spectrum in the quasi-local regime is bounded but its first derivative
$V_{g}(k)$ has an infrared divergence ($k_{sp}\to 0$), as $V_{g}(k)\sim
1/k^{2-\alpha}$, which dominate the eigenvalues of the reduced density matrix.
Inserting this into Eq. S18, we find that the second eigenvalue can be
approximated in the infinite volume limit by,

$\displaystyle\lambda_{2}$
$\displaystyle\propto\tilde{\mathcal{N}}^{2}\,t^{\frac{1}{1-z}}\sum_{m=0}^{\infty}\frac{1}{(m+R)^{\frac{z-2}{z-1}}}\cos^{2}\left[\xi_{z}\left(\frac{t}{(m+R)^{z}}\right)^{\frac{1}{1-z}}+\phi\right]$
(S19)

with $\xi_{z}=(cz)^{\frac{1}{1-z}}\left(1+c^{z+1}z^{z}\right)$. In contrast
with the earlier study of spin wave excitations for magnetic properties (see
for instance Secs. S2.1 and S4), here the site $R$ is contained in every term
of the sum over subsystem $A$, and we cannot extract the entanglement edge
without considering every term in the sum. Moreover, this means that we cannot
_a priori_ simply look at the scaling of the coefficient without taking into
account the oscillatory cosine term. However, in both our TDVP and LSWT
numerical results, we find no significant internal oscillations of the Rényi
entanglement entropies, implying that the oscillatory term does not contribute
significantly to the observed behavior of the entanglement edge. Additionally,
we find the frequency $\xi_{z}$ to be very large such that the cosine term
varies rapidly compared to the algebraically decaying amplitude. We can then
safely substitute the square cosine function by its average value, $1/2$,
allowing us to analytically perform the sum in Eq. (S19). We find

$\displaystyle\lambda_{2}$
$\displaystyle\propto\tilde{\mathcal{N}}^{2}\frac{2\pi(cz)^{\kappa}}{\absolutevalue{cz(z-1)}}t^{\frac{\kappa}{2-z}}\zeta\left(\kappa,R\right)\,,$
(S20)

where $\zeta(s,q)=\sum_{m=0}^{\infty}{(m+q)^{-s}}$ is the Hurwitz zeta
function [84, 85] and $\kappa=\frac{z-2}{z-1}=\frac{\alpha-3}{\alpha-2}$.
Using the series representation of the Hurwitz zeta function in the limit of
large $R$ and positive $\kappa$, this may be written as

$\displaystyle\lambda_{2}$
$\displaystyle=\tilde{\mathcal{N}}^{2}\frac{2\pi(cz)^{\kappa}}{\absolutevalue{cz(z-1)}}t^{\frac{\kappa}{2-z}}\left[\frac{R^{1-\kappa}}{\kappa-1}+\frac{R^{-\kappa}}{2}+\sum_{j=1}^{\infty}\frac{B_{2j}}{(2j)!}\frac{\Gamma(\kappa+2j-1)}{\Gamma(\kappa)}R^{1-\kappa-2j}\right]\,,$
(S21)

where $B_{n}$ are the Bernoulli numbers, and
$\Gamma(s)=\int_{0}^{\infty}\mathrm{e}^{-u}u^{s-1}\differential{u}$ is the
gamma function [85]. In the quasi-local regime, $\kappa>2$ and so in the limit
of large $R$ we may neglect the second and third term which contain
contributions of the order equal to $R^{-\kappa}$ and higher, and only
consider the leading term of order $R^{1-\kappa}$. It yields

$\displaystyle\lambda_{2}$
$\displaystyle\approx\tilde{\mathcal{N}}^{2}\frac{2\pi(cz)^{\kappa}}{\absolutevalue{cz(z-1)}}\frac{1}{\kappa-1}\left(\frac{t}{R}\right)^{\frac{1}{1-z}}\,$
(S22)

which scales in $R$ and $t$ with the same power. Note that the stationary
phase approximation is not valid when taking the limit $t\to\infty$ alone, and
the above equation does not hold in this limit. Consequently, any analytic
function of $\lambda_{2}$, particularly Rényi entropies, will scale in $R$ and
$t$ with the same power, so that

$\displaystyle\mathcal{S}_{n}(R,t)\simeq\frac{1}{1-n}\log\\!\big{\\{}\lambda_{2}\\!\left({R}/{t}\right)^{n}+\big{[}1-\lambda_{2}\\!\left({R}/{t}\right)\big{]}^{n}\big{\\}}\,.$
(S23)

This yields a dynamical exponent of the EE of unity, $\beta^{n}_{EE}=1\
\forall\,n$, in close agreement with our numerical results using both TDVP and
LSWT. The additional finite-size corrections neglected in moving from Eq.
(S21) to Eq. (S22) are likely responsible for the deviation from ballistic
behavior seen in our simulations on small system sizes.

#### Local regime –

The calculation in the local regime $\alpha>2$ proceeds similarly. However,
both the quasi-particle energy $E_{k}$ and group velocity $V_{\textrm{g}}(k)$
are finite over the whole Brillouin zone such that there exists a quasi-
momentum $k^{\star}$ where the group velocity is maximum,
$\max_{k}V_{g}(k)=V_{g}(k^{\star})$. The stationary phase condition in the
local regime thus has only a solution for $m/t\leq V_{g}(k^{\star})$. The
eigenvalues of the reduced density matrix are then dominated by the momenta
$k\sim k_{\mathrm{sp}}=k^{\star}$, which is independent of $R$ and $t$, and
thus scale ballistically. Finally, any analytic function of the eigenvalues
will also scale ballistically.

##  S6 Rényi entropies after a local quench

In this section we present complementary numerical results for different Rényi
entropies in the local and quasi-local regimes obtained via TDVP simulations
and LSWT calculations.

In Fig. S5(a) and (b), the counterparts of Fig. 4(b) and (c) in the main text,
dynamical exponents $\beta_{EE}^{n}$ are shown for the different values of the
Rényi parameter considered ($n=1/2$, $n=1$, $n=2$) obtained via TDVP and LSWT
with system size $N=96$ and $N=512$, respectively. The dynamical exponent of
the EE remains close to unity almost independently of the order of the Rényi
entropy, both in TDVP and LSWT. Each of the $\mathcal{S}_{n}(R,t)$ considered
displays a linear causal structure fully characterized by an entanglement edge
(EE) and its behavior is very similar to Fig. 4(a).

We have also computed the entanglement entropy on much larger system sizes
using LSWT. This allows us to obtain a better estimate of the spreading
exponent $\beta_{\textrm{\tiny EE}}^{n}$ in the thermodynamic limit. The
results are shown in Fig. S5(c)–(h) for three values of the Rényi order
$n=1/2$, $1$, and $2$ and two values of $\alpha$ for a chain of length
$N=512$. The results of the fits to the EE are shown in Fig. S5(b), where we
see that the error bars of the von Neumann entanglement entropy are much
smaller and give a more precise indication of ballistic spreading of the EE.

The determination of the EE emerges as a challenge, due to the absence of a
sharply defined edge in the quasi-local regime. In practice we chose a large
range and average over it to find a value for $\beta_{EE,n}$ in Fig. 4(b) and
(c). The threshold ranges for all Rényi entropies $\mathcal{S}_{n}$ obtained
from TDVP and LSWT in all system sizes is $\epsilon\in[20\%,80\%]$. The
associated error is taken to be the difference between the maximal and minimal
value for the dynamical exponent within this range of thresholds.

Figure S5:  EE after a local quench in the $z$-polarized phase ($h/J=50$). (a)
TDVP ($N=96$): dynamical exponent $\beta_{EE}^{n}$ of the three Rényi
entropies over different values of $\alpha$. (b) LSWT ($N=512$): dynamical
exponent $\beta_{EE}^{n}$ of the three Rényi entropies over different values
of $\alpha$. (c)–(h) LSWT ($N=512$): Space-time behavior of the three Rényi
entropies $\mathcal{S}_{n}(R,t)$ considered in the local regime ($\alpha=2.5$)
with $n=1/2$ (c), $n=1$ (d), and $n=2$ (e) as well as in the quasi-local
regime ($\alpha=1.5$) with $n=1/2$ (f), $n=1$ (g), and $n=2$ (h). The solid
green lines represent power law fits to the entropy edges with dynamical
exponents, (c) $\beta_{EE}^{n=1/2}=1.006\pm 0.001$, (d)
$\beta_{EE}^{n=1}=0.9907\pm 0.0002$, (e) $\beta_{EE}^{n=2}=1.014\pm 0.002$,
and (f) $\beta_{EE}^{n=1/2}=1.018\pm 0.003$, (g) $\beta_{EE}^{n=1}=1.018\pm
0.001$, (g) $\beta_{EE}^{n=2}=1.031\pm 0.002$.