# Universal Scaling Laws for Correlation Spreading in Quantum Systems with Short- and Long-Range Interactions

Lorenzo Cevolani, Julien Despres, Giuseppe Carleo, Luca Tagliacozzo, Laurent Sanchez-Palencia

► To cite this version:

Lorenzo Cevolani, Julien Despres, Giuseppe Carleo, Luca Tagliacozzo, Laurent Sanchez-Palencia. Universal Scaling Laws for Correlation Spreading in Quantum Systems with Short- and Long-Range Interactions. Physical Review B, 2018, 98 (2), 10.1103/PhysRevB.98.024302. [https://arxiv.org/abs/1708.03047](https://arxiv.org/abs/1708.03047)

HAL Id: hal-01532720

https://hal.science/hal-01532720v4

Submitted on 9 Jul 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Universal scaling laws for correlation spreading in quantum systems with short- and long-range interactions

Lorenzo Cevolani,¹ Julien Despres,² Giuseppe Carleo,³,⁴ Luca Tagliacozzo,⁵ and Laurent Sanchez-Palencia²

¹Institut für Theoretische Physik, Georg-August-Universität Göttingen, 37077 Göttingen, Germany

²Centre de Physique Théorique, Ecole Polytechnique,
CNRS, Univ Paris-Saclay, F-91128 Palaiseau, France

³Institute for Theoretical Physics, ETH Zurich, Wolfgang-Pauli-Str. 27, 8093 Zurich, Switzerland

⁴Center for Computational Quantum Physics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA

⁵Department of Physics and SUPA, University of Strathclyde, Glasgow G4 0NG, United Kingdom
(Dated: July 9, 2018)

The spreading of correlations after a quantum quench is studied in a wide class of lattice systems, with short- and long-range interactions. Using a unifying quasi-particle framework, we unveil a rich structure of the correlation cone, which encodes the footprints of several microscopic properties of the system. When the quasi-particle excitations propagate with a bounded group velocity, we show that the correlation edge and correlation maxima move with different velocities that we derive. For systems with a divergent group velocity, especially relevant for long-range interacting systems, the correlation edge propagates slower than ballistic. In contrast, the correlation maxima propagate faster than ballistic in gapless systems but ballistic in gapped systems. Our results shed new light on existing experimental and numerical observations, and pave the way to the next generation of experiments. For instance, we argue that the dynamics of correlation maxima can be used as a witness of the elementary excitations of the system.

# I. INTRODUCTION

The ability of a quantum system to establish long-distance correlations and entanglement, and possibly equilibrate, is determined by the speed at which information can propagate within the system. For lattice models with short-range interactions, Lieb and Robinson (LR) have unveiled a bound that forms a linear causality cone beyond which information decays exponentially¹. This bound implies ballistic propagation of equal time-correlation functions² that has been observed experimentally³,⁴ and characterized numerically⁵⁻⁹. Generalized LR bounds have been derived for long-range systems where the interactions decay algebraically, $1/R^\alpha$, with the distance $R$, see Ref. 10 and 11. The related experiments and numerical investigations have, however, lead to conflicting pictures¹²⁻¹⁸. For instance, experiments on ion chains¹⁵ and numerical simulations within truncated Wigner approximation¹⁹ for the one-dimensional (1D) long-range XY (LRXY) model point towards bounded, super-ballistic, propagation for all values of $\alpha$. In contrast, experiments on the long-range transverse Ising (LRTI) model reported ballistic propagation of correlation maxima with, however, observable leaks that increase when $\alpha$ decreases¹⁴. Moreover, time-dependent density matrix renormalization group (t-DMRG) and variational Monte-Carlo (t-VMC) numerical simulations indicate the existence of three distinct regimes, namely instantaneous, sub-ballistic, and ballistic, for increasing values of the exponent $\alpha$, see Ref. 12, 13, 16–18, and 20.

In this paper we shed light on these apparent contradictions. We focus on equal-time correlation functions that are relevant experimentally. Implications of bounded correlation spreading on universal LR bounds are not yet completely understood (see however Ref. 21) so we do

not draw explicit conclusions of the latter.

Using a universal picture based on quasi-particles that can be applied to both short- and long-range models, we unveil a double causality structure for correlation spreading. The outer structure determines the correlation edge (CE), while the inner structure determines the propagation of local extrema. For short-range interactions, the two structures are determined by the dispersion relation and can be associated to, respectively, the group and phase velocities of the quasi-particles. For long-range interactions, the inner structure is still determined by the dispersion relation. It is super-ballistic for gapless models and ballistic for gapped models. It implies that quantum quenches can be used experimentally as a witness to detect the presence of the gap and the value of the dynamical exponent of the underlying model, something that as far as we know was not realized previously. The outer structure depends both on the dispersion relation and on the considered observable, and is thus less universal. Besides pathologic cases, it is always sub-ballistic.

The identification of this double structure (i.e. edge versus local maxima of correlations), and the lack of universality of the outer-edge in long-range systems permits to accomodate and explain previous observations in a unified picture. In particular it constitutes an important result to predict the spreading of specific observables and design the next generation of experiments within a large class of long-range systems, e.g. Rydberg gases²²⁻²⁵, nonlinear optical media²⁶, polar molecules²⁷⁻²⁹, magnetic atoms³⁰⁻³⁴, superconductors³⁵, ion chains³⁶⁻⁴⁰, and solid-state defects⁴¹⁻⁴³.
## II. TIME EVOLUTION OF LOCAL CORRELATIONS

Consider a quantum system defined on a hypercubic lattice of dimension $D$ and governed by a translation-invariant Hamiltonian of the form

$$ \hat{H} = \sum_{\mathbf{R}} h(\mathbf{R}) \hat{K}_1(\mathbf{R}) + \sum_{\mathbf{R},\mathbf{R}'} J(\mathbf{R}, \mathbf{R}') \hat{K}_2(\mathbf{R}, \mathbf{R}'), \quad (1) $$

where $\mathbf{R}$ and $\mathbf{R}'$ span the lattice sites. The first term accounts for local interactions and the second term for two-site couplings. It applies to a variety of models including the Bose-Hubbard (BH), see App. A, the LRXY, and the LRTI models, see App. B, that we consider in the following. We start from the ground-state of $\hat{H}$ and quench the system out of equilibrium by changing the couplings at time $t=0$. We characterize the evolution by computing equal-time connected correlation functions with respect to the pre-quench equilibrium value. They read $G(\mathbf{R}, t) \equiv G_0(\mathbf{R}, t) - G_0(\mathbf{R}, 0)$ with $G_0(\mathbf{R}, t) \equiv \langle \hat{A}_X(t) \hat{B}_Y(t) \rangle - \langle \hat{A}_X(t) \rangle \langle \hat{B}_Y(t) \rangle$, where $\hat{A}_X$ and $\hat{B}_Y$ are local operators with support on regions $X$ and $Y$ separated by $\mathbf{R}$. Such correlations can be measured in state-of-the-art experiments$^{3,4,14,15}$. We describe quenches where the dynamics is driven by the low-energy sector of $\hat{H}$ that may be assumed to consist of quasi-particle excitations. Due to translation invariance, they are characterized by well-defined quasi-momentum $\mathbf{k}$ and energy $E_{\mathbf{k}}$. The correlation functions may be written

$$ G(\mathbf{R}, t) = g(\mathbf{R}) - \int_B \frac{d\mathbf{k}}{(2\pi)^D} \mathcal{F}(\mathbf{k}) \frac{e^{i(\mathbf{k}\cdot\mathbf{R} + 2E_{\mathbf{k}}t)} + e^{i(\mathbf{k}\cdot\mathbf{R} - 2E_{\mathbf{k}}t)}}{2}, \quad (2) $$

where the integral spans the first Brillouin zone $\mathcal{B}$. The quantity $g(\mathbf{R})$ can be dropped, since it does not depend on time. Equation (2) represents the motion of counter-propagating quasi-particle pairs, with velocities determined by $E_{\mathbf{k}}$, and where the amplitude $\mathcal{F}(\mathbf{k})$ encodes the overlap of the initial state with the quasi-particle wave functions and the matrix elements of $\hat{A}$ and $\hat{B}$. It can be derived explicitly in exactly-solvable models and quadratic systems, which can be diagonalised by means of canonical transformations. Many models, in various regimes, can be mapped into this form (see for instance Refs. 6, 12, 16-18, 44, and 45 in the context of out-of-equilibrium dynamics). The concept of quasi-particles also applies to non exactly-solvable models, where they can be determined using tensor-network techniques$^{46,47}$ for instance, and we expect that our results also hold for such systems.

## III. SHORT-RANGE COUPLINGS

Consider first the case of nearest-neighbour interactions for which the quasi-particle group velocity is

bounded. In the infinite time and distance limit along the line $R/t = \text{const}$, the integral in Eq. (2) is dominated by the momentum contributions with a stationary phase (sp), i.e. $\nabla_k(kR \mp 2E_k t) = 0$ or, equivalently,

$$ 2V_g(k_{sp}) = \pm R/t, \quad (3) $$

where $V_g = \nabla_k E_k$ is the group velocity. Since the latter is upper bounded by some value $V_g^*$, Eq. (3) has a solution only for $R/t < 2V_g^*$. The correlation function then reads$^{48}$

$$ G(R, t) \propto \frac{\mathcal{F}(k_{sp})}{(|\nabla_k^2 E_{k_{sp}}|t)^{\frac{D}{2}}} \cos \left( k_{sp}R - 2E_{k_{sp}}t + \frac{\pi}{4} \right). \quad (4) $$

For $R/t > 2V_g^*$, Eq. (3) has no solution and $G(R, t)$ is vanishingly small. The correlations are thus activated ballistically at the time $t = R/2V_g^*$. It defines a linear correlation edge (CE) with velocity $V_{CE} = 2V_g^*$, consistently with the Calabrese-Cardy picture$^{49}$.

Yet, Eq. (4) does not only yield the CE but also a series of local maxima. In the vicinity of the CE cone, only the quasi-particles with momenta $k \simeq k^*$, which move at $V_g^*$, contribute to the correlations. There the maxima (m), defined by the equation $k^*R - 2E_{k^*}t = \text{const}$, propagate at the velocity $V_m = 2V_\varphi^* \equiv 2E_{k^*}/k^*$, i.e. twice the phase velocity at the maximum of the group velocity, $k^*$. Since the phase and group velocities are generally different, the CE is expected to feature a double structure characterized by these two velocities. This observation and its counterpart for long-range systems (see below) have fundamental consequences on correlation spreading and is the pivotal result of this work.

To illustrate it, let us consider the BH model. In the superfluid regime, the dispersion relation is bounded and the group velocity has a local maximum at some momentum $0 < k^* < \pi$, see inset of Fig. 1(a). The main panel of Fig. 1(a) shows the connected one-body correlation function versus distance and time in this regime. Its value is determined from numerical integration of Eq. (2) with the coefficients calculated using Bogoliubov theory. The latter holds for weak interactions, $Jn \gg U$, with $J$ the hopping and $U$ the interaction strength, see App. A for details. As expected, the correlation cone is determined by the velocity $V_{CE} \simeq 2V_g^*$ (solid green line). Moreover, the correlations show a series of local maxima, all propagating at the same speed, approximately twice the phase velocity at the momentum $k^*$, $V_m \simeq 2V_\varphi^*$ (dashed blue lines). These observations are confirmed quantitatively in Fig. 1(c) where we compare the values of the velocities found from fits to the correlation edge ($V_{CE}$) and local maxima ($V_m$) on the one hand, to twice the group ($V_g^*$) and phase ($V_\varphi^*$) velocities at $k^*$ on the other hand.

The distinction between the edge propagating at $V_g^*$ and the maxima emanating from it and propagating at $V_\varphi^*$ permits to understand previously unexplained observations. The propagation velocity extracted from t-VMC calculations in Ref. 7 quantitatively agrees with the value $2V_\varphi^*$ calculated here. Our analysis shows that it should
<figure><img src="image_2.png" /></figure>

Figure 1. Upper panel: Spreading of the connected one-body correlation function $G(R,t) = \langle \hat{a}_R^\dagger(t) \hat{a}_0(t) \rangle - \langle \hat{a}_R^\dagger(0) \hat{a}_0(0) \rangle$ for the 1D Bose-Hubbard model. (a) Superfluid phase for a quench from the initial value $U_i n = J$ to the final value $U_f n = 0.5 J$. (b) Mott-insulator phase with $n = 1$ for a quench from $U_i = \infty$ to $U_f = 18 J$. The solid green and dashed blue lines indicate ballistic spreading at twice the maximum group velocity, $2V_g^*$, and twice the corresponding phase velocity, $2V_\varphi^*$, respectively. Lower panel: Comparison between the maximum group velocity ($V_g^*$, solid green line), the corresponding phase velocity ($V_\varphi^*$, dashed blue line), the sound velocity (c, dotted purple line), and fits to the LR cone velocity ($V_{CE}$, green diamonds) and to the velocity of the maxima ($V_m$, blue disks) for the (c) superfluid and (d) Mott insulator phases with the same initial values as for (a) and (b).

thus be assimilated to the propagation of the local maxima, i.e. the inner structure of the causal region$^{50}$. In contrast, the CE is determined by the raise of the envelop of these maxima, and moves at the velocity $2V_g^*$. The same analysis applies to the results of t-DMRG calculations for a quench in the superfluid regime of the Fermi-Hubbard model$^5$.

So far experimental characterization of correlation spreading for a quench in the Mott-insulator regime was performed close to the critical point where there is a single characteristic velocity and as a consequence no inner structure was observed$^3$. A richer behaviour occurs deeper in the Mott regime, where the presence of a gap permits to find a regime, $U > \pi(2n + 1)J$, where $V_\varphi^* > V_g^*$, see inset of Fig. 1(b). Here the excitation spectrum is found using strong-coupling perturbation theory, which holds for $n \in \mathbb{N}^*$ and $U \gg Jn$, see App. A. The connected one-body correlation function plotted in Fig. 1(b) is found from numerical integration of Eq. (2) with the corresponding coefficients. In this case, the local maxima propagate (still at $V_\varphi^*$) faster than the correlation cone (still at $V_g^*$), and vanish when reaching it, see Figs. 1(b) and (d).

<figure><img src="image_3.png" /></figure>

Figure 2. Spreading of the connected spin correlation function $G(R,t) = G_0(R,t) - G_0(R,0)$ for the following 1D models: (a) LRXY model with $\alpha = 2.3$, $G_0(R,t) = \langle S_R^z(t) S_0^z(t) \rangle - \langle S_R^z(t) \rangle \langle S_0^z(t) \rangle$ for a quench from the ground state of the XXZ model ($\epsilon = 0.2$, see App. B); (b) LRTI model with $\alpha = 1.7$, $G_0(R,t) = \langle S_R^x(t) S_0^x(t) \rangle - \langle S_R^x(t) \rangle \langle S_0^x(t) \rangle$ for the quench in the polarized phase, from $J_i/h = 0.02$ to $J_f/h = 1$, see Ref. 52. They feature a double algebraic structure (straight lines in log-log scale): A sub-ballistic correlation edge (solid green line) and super-ballistic or ballistic spreading of local maxima (dashed blue lines). The white dotted line indicates ballistic spreading for reference. The light blue, dashed lines are guides to the eye.

## IV. LONG-RANGE COUPLINGS

We now turn to long-range systems with power-law couplings, $J_{R,R'} \sim J/|\mathbf{R} - \mathbf{R}'|^\alpha$. We assume that the spectrum is regular in the whole Brillouin zone, except for a cusp at $k = 0$. There, the dispersion relation may be written $E_k \simeq \Delta + ck^z$, with $z$ the dynamical exponent and $\Delta$ the (possibly vanishing) gap. For $0 < z < 1$ the quasi-particle energy $E_k$ is bounded but the group velocity $V_g(k)$ diverges. In the following, we consider connected spin correlation functions for two spin models, as found from Eq. (2) and linear spin wave theory, see App. B and references therein. All quenches are performed in a single polarized phase, without crossing any critical line.

Figure 2(a) corresponds to the LRXY model. Owing to continuous spin rotation symmetry, it is gapless, $\Delta = 0$, and $z = (\alpha - D)/2$ for $D < \alpha < D + 2$, see Ref. 51. Figure 2(b) corresponds to the LRTI model, where the transverse magnetic field opens a gap, $\Delta > 0$, and $z = \alpha - D$ for $D < \alpha < D + 1$, see Ref. 12, 16, and 17. For both models, we find a double structure reminiscent of the one of short-range models, although with crucial differences. First, the CE is not linear but algebraic (note the log-log scales in Fig. 2). While the known extended LR bounds$^{10,11}$ are all super-ballistic, we find a sub-ballistic CE, $t \sim R^{\beta_{CE}}$ with $\beta_{CE} > 1$ (the edges are marked by solid green lines and, for reference, ballistic spreading by white dotted lines). Second, the inner structure shows a strongly model-dependent behaviour: For the LRXY model [Fig. 2(a)], the correlation maxima (dashed blue lines) are super-ballistic, $t \sim R^{\beta_m}$ with $\beta_m < 1$ while for the LRTI model [Fig. 2(b)] they are ballistic, $t \sim R$.

To understand these behaviours, let us use again the stationary-phase approximation. Equations (3) and
(4) still hold. However, the group velocity, $V_g(k) = |c|z/k^{1-z}$, now diverges at $k \to 0$. Hence, for any combination of $t$ and $R$, there is a quasi-particle with the corresponding group velocity, at the momentum $k_{sp} = (2|c|zt/R)^{1/(1-z)}$. The CE is thus dominated by the infrared divergence, where we now need to analyse the amplitude function $\mathcal{F}$. Inserting the assumed scaling $\mathcal{F}(k) \sim k^\nu$, with $\nu \ge 0$, into Eqs. (3) and (4), we find

$$G_c(R, t) \propto \frac{t^\gamma}{R^\chi} \cos \left[ A_z \left( \frac{t}{R^z} \right)^{\frac{1}{1-z}} - 2\Delta t + \frac{\pi}{4} \right], \quad (5)$$

with $\gamma = \frac{\nu+D/2}{1-z}$, $\chi = \frac{\nu+D(2-z)/2}{1-z}$, and $A_z = 2|c|(1-z)(2|c|z)^{\frac{1-z}{1-z}}$. The CE is found by imposing that the amplitude of the correlation function becomes of order one. It yields to the algebraic form

$$t^* \propto R^{\beta_{CE}}, \quad \beta_{CE} = \chi/\gamma. \qquad (6)$$

Hence, the scaling of the CE does not depend only on the dynamical exponent $z$ but also on the specific correlation function, via the exponent $\nu$, and on the dimension $D$. This contrasts with the short-range case, where a ballistic propagation independent of the dimension and of the observable is found$^{53}$. Since $\chi = \gamma + D/2$ the CE is always sub-ballistic, $\beta_{CE} > 1$. For the LRXY model and spin-spin correlations, we have $\nu = z = (\alpha - D)/2$, which yields $\beta_{CE} = 1 + D/(2\alpha)(2 + D - \alpha)$. In the numerical calculations of Fig. 2(a), the CE is found by tracing the points in the $R-t$ plane where the correlations reach $\epsilon = 2\%$ of the maximal value. The activation time $t^*$ as a function of the distance $R$ is then fitted by a power law, $t^* \sim R^{\beta_{CE}^*}$, see details in App. D. For $\alpha = 2.3$ and $D = 1$ [Fig. 2(a)], we find $\beta_{CE}^* \approx 1.083 \pm 0.013$, in good agreement with the theoretical value $\beta_{CE} \approx 1.15$.

For the LRTI model, we have $z = \alpha - D$ and $\nu = 0$. It yields the exponent $\beta_{CE} = 2 - z = 2 + D - \alpha$ completely determined by the dynamical exponent $z$. For $\alpha = 1.7$ and $D = 1$ [Fig. 2(b)], analyzing the numerical results as before we find the CE exponent $\beta_{CE}^{\text{fit}} = 1.28 \pm 0.02$, in excellent agreement with the theoretical value $\beta_{CE} = 1.3$. Note that the general formula for $\beta_{CE}$ matches the exact result of Ref. 16 for $D = 1$ and $\alpha = 3/2$ (i.e. $z = 1/2$), also confirmed by t-VMC calculations, and it is in fair agreement with the analysis of Ref. 17 for the 1D and 2D LRTI models.

On the other hand, the inner structure of the causal region is determined by the local maxima of the cosine function in Eq. (5). It does not depend on the observable but on the presence or absence of a gap. For a gapless system ($\Delta = 0$), we find

$$t_m \propto R^{\beta_m}, \quad \beta_m = z. \qquad (7)$$

The correlation maxima are thus always super-ballistic, $\beta_m < 1$. For the LRXY model and $\alpha = 2.3$ [Fig. 2(a)], we find the theoretical value $\beta_m = 0.65$. In the numerics, we study the internal structure of the correlation function

by tracking the position of the first local maximum as a function of time. We then fit the corresponding function by $t_m = aR^{\beta_{m}^{\text{fit}}} + b$. For the parameters of Fig. 2(a), it yields $\beta_{m}^{\text{fit}} \approx 0.634 \pm 0.014$, in excellent agreement with the theoretical value. It is also consistent with the experimental observation of super-ballistic dynamics in the 1D LRXY model realized with trapped ion chains for $\alpha > 1$, see Ref. 15 and in rough agreement with the analysis of numerical calculations performed within the truncated Wigner approximation for 1D and 2D LRXY models$^{19}$. The same result as Eq. (7) was found in Ref. 45, which appeared while completing the present work. Our analysis shows that this super-ballistic behaviour characterizes the inner structure but not the CE.

For a gapped system ($\Delta > 0$), the momentum dependence of the dispersion relation becomes irrelevant in the infrared limit and the argument of the cosine function in Eq. (5) is constant in the large $t$ and $R$ limit for $t \propto R$. It follows that the local maxima are here always ballistic, $\beta_m = 1$. This case applies to the LRTI model. It is confirmed in Fig. 2(b), where we observe that the local maxima converge to a ballistic propagation for sufficiently long times. Performing the analysis as above, we find $\beta_m^{\text{fit}} \approx 1.0045 \pm 0.0003$, in excellent agreement with the theoretical prediction. This result is consistent with the observation of ballistic motion of local maxima for the 1D LRTI model realized with trapped ion chains$^{14}$.

## V. CONCLUSIONS

In this work we have shown that the spreading of equal-time correlations has a double structure whose scaling laws can be related to different characteristic spectral properties. For short-range systems, they are readily associated to the group and phase velocities, which generally differ. For long-range systems with a diverging group velocity, the CE is observable-dependent and sub-ballistic. Close to the CE, the local maxima propagate ballistically in gapped systems and super-ballistically in gapless systems. Their observation can thus be used as an experimental footprint for the presence of a spectral gap.

This double structure can be observed experimentally. Our analysis provides just the first step of an important research problem that aims at enveiling the physical information encoded in correlation spreading and how this can be extracted in the next generation of experiments (see also$^9$ for recent results in this direction). In practice, the dynamics of the local maxima is easier to observe, and, as discussed above, our predictions are consistent with the existing observations. Our analysis shows, however, that in generic experiments characterizing the spreading of correlations the data need be interpreted very carefully. The propagation of local extrema does not characterize the correlation edge. Identifying the latter requires an accurate scaling analysis of the leaks. Existing experimental data have been collected either in a
regime of parameters where the two structures coincide³, or on small systems where quantitative analysis is obfuscated by strong finite-size effects. However the next generation of experiments based on Rydberg atoms, tampered wave-guides, and larger trapped-ion systems provide the natural setup to discern between the CE and the local features as our calculations suggest.

## ACKNOWLEDGMENTS

We thank Anton Buyskikh and Jad C. Halimeh for fruitful discussions. This research was supported by the European Commission FET-Proactive QUIC (H2020 grant No. 641122) and the Deutsche Forschungsgemeinschaft (DFG) through SFB/CRC1073 (Projects B03 and C03).

## Appendix A: Bose-Hubbard model

The Bose-Hubbard (BH) model,

$$ \hat{H} = -J \sum_{\langle \mathbf{R}, \mathbf{R}' \rangle} \left( \hat{a}_{\mathbf{R}}^{\dagger} \hat{a}_{\mathbf{R}'} + \text{H.c.} \right) + \frac{U}{2} \sum_{\mathbf{R}} \hat{n}_{\mathbf{R}} (\hat{n}_{\mathbf{R}} - 1), \quad (A1) $$

is constructed using the particle operators $\hat{K}_1(\mathbf{R}) \equiv \hat{n}_{\mathbf{R}} (\hat{n}_{\mathbf{R}} - 1)$ and $\hat{K}_2(\mathbf{R}, \mathbf{R}') \equiv -\hat{a}_{\mathbf{R}}^{\dagger} \hat{a}_{\mathbf{R}'} - \hat{a}_{\mathbf{R}'}^{\dagger} \hat{a}_{\mathbf{R}}$, where $\hat{a}_{\mathbf{R}}$ and $\hat{n}_{\mathbf{R}} = \hat{a}_{\mathbf{R}}^{\dagger} \hat{a}_{\mathbf{R}}$ are, respectively, the annihilation and number operators on the lattice site $\mathbf{R}$. The amplitudes are, respectively, the two-body interaction strength, $h(\mathbf{R}) = U/2$, and the tunnel amplitude $J(\mathbf{R}, \mathbf{R}') = J$. The Bose-Hubbard model has two phases, namely the superfluid (SF) phase for $J \gg U$ and the Mott-insulator (MI) phase for $J \ll U$. The precise critical point depends on the dimension and on the average number of particles per site. For a review, see for instance Ref. 54.

### 1. Superfluid phase

In the superfluid phase and for high-enough average particle density in 1D, $n \gg U/2J$, we may rely on Bogoliubov meanfield approximation. Assuming small density fluctuations, $\Delta n \ll n$, one develops the interaction term in Eq. (A1) up to quadratic order. The resulting quadratic form is then diagonalized using standard Bogoliubov transformation (see for instance Refs. 16 and 44). It yields the gapless dispersion relation

$$ E_k \approx 2\sqrt{2J \sin^2(k/2) [2J \sin^2(k/2) + nU]} \quad (A2) $$

It is phononic in the low-energy limit, $E_k \approx ck$, with the sound velocity $c = \sqrt{2nJU}$. At higher energy, it shows an inflection point at some finite momentum $0 < k^* < \pi$, corresponding to the maximum group velocity $V_g^* = V_g(k^*)$.

After the quench, the connected one-body correlation function, $G(R, t) = \langle \hat{a}_R^\dagger(t) \hat{a}_0(t) \rangle - \langle \hat{a}_R^\dagger(0) \hat{a}_0(0) \rangle$ considered in the main paper, is then cast into the form of Eq. (2) by mapping the particle operators onto Bogoliubov quasi-particles operators. It yields the amplitude function

$$ \mathcal{F}(k) = \frac{2n^2 J U_f (U_f - U_i)}{E_{k,i} E_{k,f}^2} \sin^2(k/2), \quad (A3) $$

where the indices 'i' and 'f' refer to the pre-quench and post-quench values, respectively.

### 2. Mott insulator phase

In the Mott insulator phase, the model develops a finite gap. The energy excitations may be found using strong-coupling expansions, see for instance Refs. 6, 55, and 56. For $n \in \mathbb{N}^*$ and $U \gg Jn$, it yields low-energy excitations made of doublon-holon pairs of energy

$$ 2E_k \approx U - 2(2n+1)J \cos(k). \quad (A4) $$

The maximum of the group velocity is at the center of the band, $k^* = \pi/2$, where the group and phase velocities are $V_g^* = (2n+1)J$ and $V_\varphi^* = U/\pi$, respectively. Hence, for $U > \pi(2n+1)J$, the phase velocity exceeds the group velocity, $V_\varphi^* > V_g^*$. Note that this regime is well inside the Mott regime where Eq. (A4) is accurate.

Similarly as for the superfluid phase, the connected one-body correlation function can be cast into the form of Eq. (2) with the amplitude function

$$ \mathcal{F}(k) = \frac{4Jn(n+1)}{iU} \sin(k), \quad (A5) $$

for the quench from $U_i = \infty$ to $U_f = U$ as considered in the paper.

## Appendix B: Long-range XY and XXZ models

For spin models, the operators $\hat{K}_j$ represent spin operators, the parameter $J(\mathbf{R}, \mathbf{R}')$ the exchange term, and $h(\mathbf{R})$ a magnetic field. For the long-range XY (LRXY) model, we use $\hat{K}_2(\mathbf{R}, \mathbf{R}') \equiv \hat{S}_\mathbf{R}^x \cdot \hat{S}_\mathbf{R}'^x + \hat{S}_\mathbf{R}^y \cdot \hat{S}_\mathbf{R}'^y$, $J(\mathbf{R}, \mathbf{R}') = -J/2 |\mathbf{R} - \mathbf{R}'|^\alpha$, and $h(\mathbf{R}) = 0$. For the initial state, it is generalized to the XXZ model by including an antiferromagnetic exchange coupling in the $z$ direction, which yields the Hamiltonian

$$ \hat{H} = \sum_{\mathbf{R} \neq \mathbf{R}'} \frac{J/2}{|\mathbf{R} - \mathbf{R}'|^\alpha} \left[ - (\hat{S}_\mathbf{R}^x \hat{S}_\mathbf{R}'^x + \hat{S}_\mathbf{R}^y \hat{S}_\mathbf{R}'^y) + \epsilon \hat{S}_\mathbf{R}^z \hat{S}_\mathbf{R}'^z \right]. \quad (B1) $$

For the LRXY case considered in the main paper, the quench is performed from the ground state of the XXZ model ($\epsilon \neq 0$) to the XY model ($\epsilon = 0$).
We study the phase where the rotational symmetry around the z axis is spontaneously broken and the spins are polarized along the x axis. There, the Hamiltonian can be diagonalized using standard Holstein-Primakoff transformation$^{57,58}$,

$$ 
\begin{aligned}
S_{\mathbf{R}}^{x} &= \frac{1}{2} - \hat{a}_{\mathbf{R}}^{\dagger} a_{\mathbf{R}}, \\
S_{\mathbf{R}}^{y} &\approx -\frac{\hat{a}_{\mathbf{R}}^{\dagger} - \hat{a}_{\mathbf{R}}}{2i}, \\
S_{\mathbf{R}}^{z} &\approx -\frac{\hat{a}_{\mathbf{R}} + \hat{a}_{\mathbf{R}}^{\dagger}}{2},
\end{aligned}
$$

where terms beyond second order in the boson operators $\hat{a}_{\mathbf{R}}$ and $\hat{a}_{\mathbf{R}}^{\dagger}$ are neglected. Inserting these transformations into Eq. (B1) yields a quadratic Bose Hamiltonian, which can be diagonalized using canonical Bogoliubov transformations, see for instance Ref. 51. For $\epsilon = 0$ (LRXY model), it yields the dispersion relation for $D = 1$

$$ E_k = \frac{JP_\alpha(0)}{2} \sqrt{1 - \frac{P_\alpha(k)}{P_\alpha(0)}}, \qquad (B2) $$

where $P_{\alpha}(k) = \int dR e^{-ik \cdot R}/|R|^{\alpha}$ is the Fourier transform of the long-range term. In the infrared limit, it can be written

$$ P_{\alpha}(k) \approx P_{\alpha}(0) + P'_{\alpha} k^{\alpha-D}, \qquad (B3) $$

where $P_{\alpha}(0)$ and $P'_{\alpha}$ are finite constants. Hence, we find $E_k \propto |k|^z$ with $z = (\alpha - D)/2$. For $D < \alpha < D + 2$, the quasi-particle energy is finite but the group velocity $V_g$ diverges in the infrared limit $k \to 0$.

The connected spin-spin correlation function along the z direction for a quench from $\epsilon_i \neq 0$ to $\epsilon_f = 0$, $G_0(R, t) = \langle S_R^z(t)S_0^z(t) \rangle - \langle S_R^z(t) \rangle \langle S_0^z(t) \rangle$, used in the paper is cast into the form of Eq. (2) using the quasi-particle amplitudes, which yields

$$ \mathcal{F}(k) = \frac{\epsilon_i P_{\alpha}(k)}{8 P_{\alpha}(0)} \sqrt{\frac{P_{\alpha}(0) - P_{\alpha}(k)}{P_{\alpha}(0) + \epsilon_i P_{\alpha}(k)}}. \qquad (B4) $$

In the infrared limit, it scales as $\mathcal{F}(k) \sim k^{\nu}$ with $\nu = (\alpha - D)/2 = z$.

The linearization of the Holstein-Primakoff transformation holds for $|1/2 - S_R^x| \ll 1$, see Ref. 58. For the calculations corresponding to Fig. 2 (a), we find $\max\{|1/2 - S_R^x|\} \simeq 0.12$. It validates the spin-wave approximation used in the paper. This result agrees with the predictions for the same model made in Ref. 51 where the validity of the spin wave approach for that model is extensively studied.

## Appendix C: Long-range transverse Ising model

The long-range transverse Ising (LRTI) model corresponds to the spin operators $\hat{K}_1(\mathbf{R}) \equiv \hat{S}_{\mathbf{R}}^z$ and

$\hat{K}_2(\mathbf{R}, \mathbf{R}') \equiv \hat{S}_{\mathbf{R}}^x \cdot \hat{S}_{\mathbf{R}'}^x$ with a uniform magnetic field $h(\mathbf{R}) = -2h$ and the algebraically decaying exchange amplitude $J(\mathbf{R}, \mathbf{R}') = 2J/|\mathbf{R} - \mathbf{R}'|^\alpha$, which yields

$$ \hat{H} = \sum_{\mathbf{R} \neq \mathbf{R}'} \frac{2J}{|\mathbf{R} - \mathbf{R}'|^\alpha} \hat{S}_{\mathbf{R}}^x \hat{S}_{\mathbf{R}'}^x - 2h \sum_{\mathbf{R}} \hat{S}_{\mathbf{R}}^z. \qquad (C1) $$

The LRTI has two phases$^{59}$. In the z-polarized phase, the dispersion relation can be found using again the Holstein-Primakoff transformation

$$ 
\begin{aligned}
S_{\mathbf{R}}^{x} &\approx \frac{\hat{a}_{\mathbf{R}} + \hat{a}_{\mathbf{R}}^{\dagger}}{2}, \\
S_{\mathbf{R}}^{y} &\approx -\frac{\hat{a}_{\mathbf{R}}^{\dagger} - \hat{a}_{\mathbf{R}}}{2i}, \\
S_{\mathbf{R}}^{z} &= \frac{1}{2} - \hat{a}_{\mathbf{R}}^{\dagger} a_{\mathbf{R}},
\end{aligned}
$$

One then finds the dispersion relation for $D = 1$

$$ E_k = 2\sqrt{h[h + JP_{\alpha}(k)]}. \qquad (C2) $$

In the infrared limit, where Eq. (B3) holds, it can be expressed as

$$ E_k = \Delta + c|k|^z, \qquad (C3) $$

where the gap $\Delta = 2\sqrt{h[h + JP_{\alpha}(0)]}$ is finite, $c = \sqrt{\frac{h}{h + JP_{\alpha}(0)}}JP'_{\alpha}$, and $z = \alpha - D$, see Ref. 17. Hence, the quasi-particle energy is finite and the group velocity diverges for $D < \alpha < D + 1$.

In the main paper, we consider the connected spin-spin correlation function along the x direction, $G_0(R, t) = \langle S_R^x(t)S_0^x(t) \rangle - \langle S_R^x(t) \rangle \langle S_0^x(t) \rangle$, which can be written in the form of Eq. (2) with

$$ \mathcal{F}(k) = \frac{h(J_i - J_f) P_{\alpha}(k)}{8 [h + J_f P_{\alpha}(k)] \sqrt{h [h + J_i P_{\alpha}(k)]}}. \qquad (C4) $$

In the infrared limit, it converges to a finite value. Hence $\mathcal{F}(k) \sim k^{\nu}$ with $\nu = 0$.

For the calculations corresponding to the Ising model in Fig. 2 (b), we find $\max\{|1/2 - S_R^z|\} = 0.11$, which validates the spin-wave approximation also in this case.

## Appendix D: Numerical analysis of the local extrema and the correlation edge for the LRXY and LRTI models

In the main paper, it is shown that the causality cone features a double structure: an outer structure, which determines the correlation edge (CE), and an inner structure where local extrema propagate. Here we provide details on the numerical analysis of the CE and of the trajectory for the first local extrema, for both the LRXY and LRTI models considered in the main manuscript.
<figure><img src="image_4.png" /></figure>

Figure 3. Spreading of the connected spin-spin correlation function for the LRXY model with $\alpha = 2.3$ (same data as in Fig. 2 (a) of the main paper; log-log scale). The filled blue points correspond for each distance $R$ to the first time where the correlation reaches 2% of its maximum value. The solid green line shows the power law predicted theoretically with a fitted multiplicative factor.

## 1. LRXY model

For the LRXY model, we consider the time evolution of the connected spin-spin correlation function $G(R, t)$ along the $z$ axis. Figure 3 shows the same data as Fig. 2 (a) in the main paper. To find the CE, we proceed as follows. For each distance $R$, we trace the activation time $t^*(R)$ corresponding to the first time when a fraction (2%) of the absolute maximum of the correlation function is reached. It yields the filled blue points on Fig. 3. The latter feature a linear trajectory in the log-log scale of the figure, that is a power law behavior in lin-lin scale. The latter is in excellent agreement with the theoretical prediction $\beta_{CE} \simeq 1.15$ shown as a solid green line on the figure. We have also fitted a power law function, $t^* \propto R^{\beta_{CE}^{\text{fit}}}$, to the blue points for $20 < R < 175$ (not shown on the figure). It yields $\beta_{CE}^{\text{fit}} = 1.083 \pm 0.013$ in good agreement with the prediction.

A similar result is obtained on a length scale closer to the one accessible in state-of-the-art experiments$^{60}$. Fitting the same algebraic function to the correlation edge in the range $10 < R < 30$ and $t < 10/J$ yields $\beta_{CE}^{\text{fit}} = 1.121 \pm 0.012$. It is already in good agreement with the fit in the larger range and with the theoretical prediction.

To analyze the behaviour of the local extrema, we trace them from the data of Fig. 3 where they are clearly visible. The result for the first one is plotted on Fig. 4 (solid red line) together with a fitted power law, $t_m = aR^{\beta_m^{\text{fit}}} + b$ (dashed blue line). The fit yields $\beta_m^{\text{fit}} = 0.634 \pm 0.014$ in excellent agreement with the theoretical exponent $\beta_m = 0.65$.

<figure><img src="image_5.png" /></figure>

Figure 4. Trajectory of the first extremum of the connected spin-spin correlation function for the LRXY model (lin-lin scale). The figure shows the numerical results found from Fig. 3 (solid red line) together with a fitted power law (dashed blue line).

<figure><img src="image_6.png" /></figure>

Figure 5. Spreading of the connected spin-spin correlation function for the LRTI model with $\alpha = 1.7$ (same data as in Fig. 2(b) of the main paper; log-log scale). The filled blue points correspond for each distance $R$ to the first time where the correlation reaches 2.8% of its maximum value. The solid green line shows the power law predicted theoretically with a fitted multiplicative factor.

## 2. LRTI model

We now turn to the LRTI model and perform the same analysis as above, up to details that we discuss below. Figure 5 shows the same data as Fig. 2(b) in the main paper but on a smaller time scale. The CE, corresponding to the trajectory of the first points where the correlation function reaches a fraction (2.8%) of its maximum (filled blue points) matches very well the theoretical prediction $t^* \propto R^{\beta_{CE}}$ with $\beta_{CE} = 1.3$ (solid green line with
<figure><img src="image_7.png" /></figure>

Figure 6. Trajectory of an extremum of the connected spin-spin correlation function for the LRTI model (lin-lin scale). It corresponds to the dashed, dark blue line in Fig. 2(b) of the main paper. The figure shows the numerical results found from Fig. 5 (solid red line) together with a fitted power law (dashed blue line).

a fitted multiplicative factor). Moreover, fitting a power

law function to the blue points for $4 < R < 30$ yields $\beta_{CE}^{\text{fit}} = 1.28 \pm 0.02$ (not shown on the figure), in perfect agreement with our prediction.

We analyse the trajectory of the maxima in the same way as before. For the LRTI model, ballistic spreading is expected sufficiently large values of $t$ and $R$. This is confirmed by the behaviour of the local maxima in Fig. 2(b) of the main paper. The first local maxima are marked by dashed light blue curves. In the figure, they are curved owing to the competition between the different terms of the phase of the cosine in Eq. (5) of the main manuscript. They, however, clearly converge to a ballistic behaviour when $t$ and $R$ increase. For the range of $t$ and $R$ presented in Fig. 2 (b) the dark blue line is purely ballistic even for small values of $t$ and $R$. Its trajectory is shown in lin-lin scale on Fig. 6. To confirm the ballistic behaviour, we fitted the power-law function $t_m = aR^{\beta_m^{\text{fit}}} + b$ to the date (dashed blue line on the figure). The fit gives $\beta_m^{\text{fit}} = 1.0045 \pm 0.0003$ in excellent agreement with the theoretical exponent $\beta_m = 1$.

We have also performed similar fits on the other local maxima, which confirms the asymptotic ballistic behaviour. For instance, for the first local maximum (lowest dashed, light-blue lines in Fig. 2(b) of the main paper) and $10 < R < 20$, we find $\beta_m^{\text{fit}} = 1.0039 \pm 0.0002$, which is also in excellent agreement with the theoretical prediction.

<footnote n="1">1 E. H. Lieb and D. W. Robinson, *The finite group velocity of quantum spin systems*, Comm. Math. Phys. **28**, 251 (1972).</footnote>

<footnote n="2">2 S. Bravyi, M. B. Hastings, and F. Verstraete, *Lieb-Robinson bounds and the generation of correlations and topological quantum order*, Phys. Rev. Lett. **97**, 050401 (2006).</footnote>

<footnote n="3">3 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).</footnote>

<footnote n="4">4 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).</footnote>

<footnote n="5">5 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).</footnote>

<footnote n="6">6 P. Barmettler, D. Poletti, M. Cheneau, and C. Kollath, *Propagation front of correlations in an interacting Bose gas*, Phys. Rev. A **85**, 053625 (2012).</footnote>

<footnote n="7">7 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).</footnote>

<footnote n="8">8 L. Dell'Anna, O. Salberger, L. Barbiero, A. Trombettoni, and V. E. Korepin, *Violation of cluster decomposition and absence of light cones in local integer and half-integer spin chains*, Phys. Rev. B **94**, 155140 (2016).</footnote>

<footnote n="9">9 M. Kormos, M. Collura, G. Takacs, and P. Calabrese, *Real-time confinement following a quantum quench to a non-integrable model*, Nat. Phys. **13**, 246 (2017).</footnote>

<footnote n="10">10 M. B. Hastings and T. Koma, *Spectral gap and exponential decay of correlations*, Comm. Math. Phys. **265**, 781 (2006).</footnote>

<footnote n="11">11 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).</footnote>

<footnote n="12">12 P. Hauke and L. Tagliacozzo, *Spread of correlations in long-range interacting quantum systems*, Phys. Rev. Lett. **111**, 207202 (2013).</footnote>

<footnote n="13">13 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).</footnote>

<footnote n="14">14 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).</footnote>

<footnote n="15">15 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).</footnote>

<footnote n="16">16 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).</footnote>

<footnote n="17">17 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).</footnote>
18 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).

19 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).

20 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).

21 M. Kastner, *Entanglement-enhanced spreading of correlations*, New J. Phys. **17**, 123024 (2015).

22 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).

23 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).

24 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).

25 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).

26 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).

27 A. Micheli, G. K. Brennen, and P. Zoller, *A toolbox for lattice-spin models with polar molecules*, Nat. Phys. **2**(5), 341 (May 2006).

28 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).

29 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).

30 A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and T. Pfau, *Bose-Einstein condensation of chromium*, Phys. Rev. Lett. **94**, 160401 (Apr 2005).

31 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).

32 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).

33 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).

34 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).

35 A. A. Houck, H. E. Tureci, and J. Koch, *On-chip quantum simulation with superconducting circuits*, Nat. Phys. **8**, 292 (2012).

36 X.-L. Deng, D. Porras, and J. I. Cirac, *Effective spin quantum phases in systems of trapped ions*, Phys. Rev. A **72**, 063407 (2005).

37 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).

38 B. P. Lanyon, C. Hempel, D. Nigg, M. Müller, R. Geritsma, 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).

39 C. Schneider, D. Porras, and T. Schaetz, *Experimental quantum simulations of many-body physics with trapped ions*, Rep. Prog. Phys. **75**(2), 024401 (2012).

40 R. Blatt and C. F. Roos, *Quantum simulations with trapped ions*, Nat. Phys. **8**, 277 (2012).

41 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).

42 G. Balasubramanian, P. Neumann, D. Twitchen, M. Markham, R. Kolesov, N. Mizuuchi, J. Isoya, J. Achard, J. Beck, J. Tissler, et al., *Ultralong spin coherence time in isotopically engineered diamond*, Nat. Mater. **8**, 383 (2009).

43 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).

44 S. S. Natu and E. J. Mueller, *Dynamics of correlations in a dilute Bose gas following an interaction quench*, Phys. Rev. A **87**, 053607 (2013).

45 I. Frérot, P. Naldesi, and T. Roscilde, *Multispeed prethermalization in quantum spin models with power-law decaying interactions* (2018).

46 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).

47 N. Nakatani, S. Wouters, D. V. Neck, and G. K.-L. Chan, *Linear response theory for the density matrix renormalization group: Efficient algorithms for strongly correlated excited states*, J. Chem. Phys. **140**(2), 024108 (2014).

48 For simplicity, we assume that, if it exists, there is a single solution $k_{\text{sp}}$. Should Eq. (3) have several solutions, Eq. (4) is replaced by the sum of the corresponding contributions$^{16}$.

49 P. Calabrese and J. Cardy, *Time dependence of correlation functions following a quantum quench*, Phys. Rev. Lett. **96**, 136801 (2006).

50 No clear double structure is visible on Fig. 1(a) of Ref. 7 because it corresponds to $Un = 4$ where the group and phase velocities are almost equal.

51 I. Frérot, P. Naldesi, and T. Roscilde, *Entanglement and fluctuations in the XXZ model with power-law interactions*, Phys. Rev. B **95**, 245111 (2017).

52 For the LRTI model with $\alpha = 1.7$, the phase transition between the polarized and Néel phases is located around $J/h \simeq 3$, see Ref. 59.

53 The dependence on $D$ and $\nu$ completely disappears in the limit of a spectrum with a non-divergent velocity, $z \to 1^{-}$, in agreement with the general LR picture$^{1}$.

54 A. Georges and T. Giamarchi, *Many-Body Physics with Ultracold Gases*, C. Salomon, G. V. Shlyapnikov, and L. F. Cugliandolo eds., *Proceedings of the Les Houches Summer*
*School of Theoretical Physics, Vol. XCIV* (Oxford Univ. presse, Oxford, 2012), chap. Strongly correlated bosons and fermions in optical lattices.

55 E. Altman and A. Auerbach, *Oscillating superfluidity of bosons in optical lattices*, Phys. Rev. Lett. **89**, 250404 (2002).

56 S. D. Huber, E. Altman, H. P. Büchler, and G. Blatter, *Dynamical properties of ultracold bosons in an optical lattice*, Phys. Rev. B **75**, 085106 (2007).

57 T. Holstein and H. Primakoff, *Field dependence of the intrinsic domain magnetization of a ferromagnet*, Phys. Rev. **58**, 1098 (1940).

58 A. Auerbach, *Interacting Electrons and Quantum Magnetism* (Springer, New York, 1994).

59 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).

60 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).