Multi-speed prethermalization in spin models with power-law decaying interactions

TL;DR

This paper explores how long-range power-law interactions in spin models cause super-ballistic entanglement and correlation propagation, significantly altering relaxation dynamics compared to short-range systems.

Contribution

It reveals that in certain regimes, long-range interactions induce non-linear dispersion and super-ballistic propagation, a novel insight into quantum relaxation dynamics.

Findings

Super-linear growth of correlation fronts over time.

Sub-linear relaxation times for subsystem observables.

Extreme wavelength dependence of relaxation times.

Abstract

The relaxation of uniform quantum systems with finite-range interactions after a quench is generically driven by the ballistic propagation of long-lived quasi-particle excitations triggered by a sufficiently small quench. Here we investigate the case of long-range ($1/r^{\alpha}$) interactions for $d$-dimensional lattice spin models with uniaxial symmetry, and show that, in the regime $d < \alpha < d+2$, the entanglement and correlation buildup is radically altered by the existence of a non-linear dispersion relation of quasi-particles, $\omega\sim k^{z<1}$, at small wave vectors, leading to a divergence of the quasiparticle group velocity and super-ballistic propagation. This translates in a super-linear growth of correlation fronts with time, and sub-linear growth of relaxation times of subsystem observables with size, when focusing on $k=0$ fluctuations. Yet the large dispersion in group velocities leads to an extreme wavelength dependence of relaxation times of finite-$k$ fluctuations, with entanglement being susceptible to the longest of them. Our predictions are directly relevant to current experiments probing the nonequilibrium dynamics of trapped ions, or ultracold magnetic and Rydberg atoms in optical lattices.

Equations

$${\mathcal{H}}_{\alpha \mathrm{XX}}=-J\sum _{i\ne j}|{\mathbf{r}}_i-{\mathbf{r}}_j{|}^{-\alpha }\left(S_i^x{S}_j^x+S_i^y{S}_j^y\right).$$

$${\omega }_{\mathbf{k}}=2sJ{\gamma }_0\sqrt{1-{\gamma }_k/{\gamma }_0}$$

$$S^{zz}(\mathbf{k},t)=g_{\mathbf{k}}^{(0)}+g_{\mathbf{k}}^{(1)}\cos (2{\omega }_{\mathbf{k}}t)$$

$$S_A^{zz}(\mathbf{k},t)=l_A^{-d}\sum _{\mathbf{q}}{f}_A(\mathbf{q})S^{zz}(\mathbf{k}+\mathbf{q},t).$$

$$S_A^{zz}(\mathbf{k},t)\approx g_{\mathbf{k}}^{(0)}+g_{\mathbf{k}}^{(1)}\cos (2{\omega }_{\mathbf{k}}t)\prod _{a=1}^d\mathrm{sinc}\left(\pi t/t_{\mathbf{k},a}^{\ast }\right)$$

$$E_A(t)=\sum _{\mathbf{k}}f(t/t_{\mathbf{k}}^{\ast })s_{\mathrm{eq}}(\mathbf{k})$$

$$⟨A⟩(t)=⟨A⟩(-t).$$

$$S_A(t)=S_A(-t)$$

$$⟨A⟩(t)=\sum _{ab}⟨{\psi }_a|A|{\psi }_b⟩{\Psi }_a{\Psi }_b\left[\cos ({\omega }_{ab}t)+i\sin ({\omega }_{ab}t)\right]$$

$$|\Psi (0)⟩=\sum _a{\Psi }_a|{\psi }_a⟩$$

$$|\Psi (t)⟩=\sum _p{\lambda }_p(t)|a_p(t)⟩\otimes |b_p(t)⟩$$

$$S_i^x$$

$$S_i^y$$

$$S_i^z$$

$${\mathcal{H}}^{(2)}=Js\sum _{i,j}(b_i^{†}{A}_{ij}{b}_j+b_i{A}_{ij}{b}_j^{†}+b_i{B}_{ij}{b}_j+b_i^{†}{B}_{ij}{b}_j^{†})$$

$$A_{ij}$$

$$B_{ij}$$

$${\mathcal{H}}^{(2)}=(1/2)\sum _k\left(\begin{array}{cc}{b}_{\mathbf{k}}^{†}& b_{-\mathbf{k}}\end{array}\right)\left(\begin{array}{cc}{A}_{\mathbf{k}}& B_{\mathbf{k}}\\ B_{\mathbf{k}}& A_{\mathbf{k}}\end{array}\right)\left(\begin{array}{c}{b}_{\mathbf{k}}\\ b_{-\mathbf{k}}^{†}\end{array}\right)$$

$${\mathcal{H}}_{\mathrm{eff}}=\sum _{\mathbf{k}}{\omega }_{\mathbf{k}}{\beta }_{\mathbf{k}}^{†}{\beta }_{\mathbf{k}}.$$

$$b_{\mathbf{k}}(t)$$

$$⟨b_{\mathbf{k}}^{†}{b}_{\mathbf{k}}⟩(t)$$

$$⟨b_{\mathbf{k}}{b}_{-\mathbf{k}}⟩(t)$$

$$\frac{c(\alpha )}{2Js}=\underset{\delta k\to 0}{\lim }\frac{1}{\delta k}\sqrt{{\gamma }_0}\sqrt{{\gamma }_0-{\gamma }_{(\delta k,0,\mathrm{0..})}}$$

$$=\underset{\delta k\to 0}{\lim }\frac{1}{\delta k}\sqrt{{\gamma }_0}{\left(\frac{(\delta k{)}^2}{2}\sum _{\mathbf{r}}\frac{x^2}{r^{\alpha }}+...\right)}^{1/2}\sim \sqrt{\frac{1}{\alpha -d-2}}$$

$$r(t)=\frac{1}{2sN}\sum _{\mathbf{k}}⟨b_{\mathbf{k}}^{†}{b}_{\mathbf{k}}⟩(t)\ll 1.$$

$$r(\infty )$$

$$\approx$$

$$\frac{1}{N}\sum _{\mathbf{k}}(u_{\mathbf{k}}^2+v_{\mathbf{k}}^2)⟨{\beta }_{\mathbf{k}}^{†}{\beta }_{\mathbf{k}}⟩\approx \frac{k_{B}T}{N}\sum _{\mathbf{k}}\frac{A_{\mathbf{k}}}{{\omega }_{\mathbf{k}}^2}$$

$$S_k^{zz}$$

$$S_k^{yy}(t)=\frac{s}{2}\left[1+\frac{{\gamma }_{\mathbf{k}}}{{\gamma }_0-{\gamma }_{\mathbf{k}}}{\sin }^2({\omega }_{\mathbf{k}}t)\right]$$

$$\mathcal{A}={\int }_0^{\beta }d\tau \int d^{d}x\left(|{\partial }_{\tau }\psi {|}^2+|\nabla \psi {|}^2+r|\psi {|}^2+u|\psi {|}^4\right).$$

$$\int d\tau \sum _{ij}{\psi }_i^{\ast }(\tau ){\tilde{J}}_{ij}^{-1}{\psi }_j(\tau )$$

$${\tilde{J}}_{ij}^{-1}=\frac{1}{N}\sum _{\mathbf{k}}\frac{e^{i\mathbf{k}\cdot ({\mathbf{r}}_i-{\mathbf{r}}_j)}}{K+{\gamma }_{\mathbf{k}}}$$

Full text

Introduction. The relaxation of the pure state of a generic quantum many-body system following a quantum quench is a complex phenomenon: there each subsystem sees its reduced state reach thermal equilibrium (or a non-equilibrium steady-state) via the drastic rearrangement of correlations with its exterior, acting as a bath lindenetal2009 ; dalessionetal2016 . A particularly dramatic case – yet highly relevant experimentally – concerns the choice of a factorized initial state, in which each subsystem admits a pure-state (zero-entropy) description. The buildup of entropy and fluctuations towards equilibration requires therefore the growth of entanglement and correlations between subsystems, namely the flow of information from each subsystem to its exterior. The post-quench dynamics of entanglement and correlations strongly depends a priori on the peculiar system under investigation. Yet, in generic systems with short-range (SR) interactions and translational invariance, a seemingly universal picture emerges for sufficiently small quenches, as shown both theoretically (see Ref. CalabreseC2016 and refs therein) as well as experimentally cheneauetal2012 ; Langen2013 . In such a picture (see Fig. 1(a) for a sketch) correlations and entanglement establish because of the ballistic propagation of quasiparticle excitations (QP) triggered by the quench. As a consequence 1) correlations reorganize over distances which grow linearly in time (within a so-called “causal cone” or “light cone”) consistent with Lieb-Robinson bounds liebR1972 on the propagation of signals in quantum systems – and, as a consequence, subsystem fluctuations (which are integrals of correlation functions) grow linearly in time as well. The velocity of correlation fronts is twice the maximum group velocity of the QP – this is consistent with correlations among two points being established by counter-propagating QP originating from the midpoint; 2) entanglement entropy (EE) grows linearly in time as a consequence of ballistic QP traversing a subsystem calabreseC2005 .

Very interestingly, the above universal picture of post-quench dynamics breaks down in the presence of long-range (LR) interactions which decay as a power-law $1/r^{\alpha}$ of the distance $r$ – stirring a considerable recent activity both in theory haukeT2013 ; schachenmayeretal2013 ; eisertetal2013 ; hazzardetal2014 ; foss-feigetal2015 ; schachenmayeretal2015 ; cevolanietal2015 ; buyskikhetal2016 ; maghrebietal2016 and experiments richermeetal2014 ; Jurcevic2014 . In particular, in contrast with SR interactions, one generally observes that 1) correlation fronts can move in a super-ballistic fashion (namely faster than linearly in time) for sufficiently small values of $\alpha$; 2) the growth of EE can be sub-linear in time (and as slow as logarithmic for small $\alpha$). This very rich phenomenology lacks universal traits, and it escapes so far a unifying picture; as an example, rigorous generalizations of the Lieb-Robinson bounds to LR interactions Hastings2006 ; Gongetal2014 ; foss-feigetal2015 reveal to be loose compared to the actual observations on specific systems. At the same time the importance of closed, nonequilibrium LR-interacting systems is augmented by their experimental realization in atomic physics via trapped ions richermeetal2014 ; Jurcevic2014 , magnetic atoms dePazetal2013 , Rydberg atoms Browaeysetal2016 and dipolar molecules Moses2017 .

Inspired by the potential of the above cited experimental setups to realize quantum spin models with LR interactions, we provide in this work a complete picture of the quench dynamics in paradigmatic XX spin models with $1/r^{\alpha}$ interactions in arbitrary ($d$) dimensions. Indeed we show that a QP description – via linear spin-wave (LSW) theory, corroborated in $d=1$ by density-matrix renormalization group (DMRG) calculations – allows one to achieve a unifying picture for the dynamics of correlations, fluctuations and entanglement. The contradiction between the fast dynamics of correlations and fluctuations and the slow growth of entanglement is resolved by a scenario of multi-speed prethermalization 222We use the term prethermalization here as our QP treatment of the system can only describe dephasing between independent modes, while true thermalization is achieved over longer time scales due to the non-linear coupling terms among QP modes. (see sketch in Fig. 1(b)), in which LR interactions lead to the divergence of the maximum group velocity, so that the velocities of QPs moving with different wavevectors span an infinite range. Different Fourier components of local fluctuations possess relaxation times which scale as the inverse of the corresponding group velocities, allowing one to reconstruct (both theoretically as well as experimentally) the multi-speed nature of the prethermalization process. The growth of the EE, on the other hand, is sensitive to the dynamics of fluctuations at all wavelengths, and as a consequence it is much slower than the evolution of the fastest fluctuations.

Model and methods.— We center our analysis on the quantum XX model with $1/r^{\alpha}$ interactions (hereafter denoted as $\alpha$XX model), describing quantum $s$-spins localized at the nodes of a $d-$dimensional lattice, interacting with a power-law decaying exchange interaction :

[TABLE]

Here $S_{i}^{\beta}$ ($\beta=x,y,z$) are quantum spin operators defined on the sites of a $d$-dimensional hypercubic lattice. The results presented here are restricted to the case $s=1/2$, although they are immediately generalizable to arbitrary $s$ values. Special emphasis is put on the case $\alpha=3$, relevant to dipolar systems dePazetal2013 ; hazzardetal2014 ; Browaeysetal2016 . We focus our study on the unitary quench dynamics $|\Psi(t)\rangle=\exp(-i{\cal H}_{\alpha\rm XX}t)|\Psi(0)\rangle$ generated by the $\alpha$XX Hamiltonian starting from the factorized state $|\Psi(0)\rangle=|\uparrow_{x}\rangle^{\otimes N}$ of spins polarized along the $x$ axis, which is both of fundamental interest, as well as of practical relevance for experiments 333In the DMRG calculations we consider the above state projected onto the $S_{\rm tot}^{z}=0$ subspace, in order to decrease the size of the Hilbert space visited by the evolved state.. When considering the case of ferromagnetic interactions $J>0$, one immediately recognizes that $|\Psi(0)\rangle$ is the mean-field ground state of the Hamiltonian – which is expected to be an increasingly good approximation to the exact ground state of the $\alpha$XX model for large $d$ and small $\alpha$ paper1-XX_LR . As a consequence, the energy injected by the quench leads the system to relax towards a low-temperature regime which exhibits long-range ferromagnetic order in the $xy$ plane – this is true in $d=3$ for any $\alpha$, while it requires $\alpha<2d$ in $d=1,2$ to stabilize long-range order at low but finite temperatures SuppMat 444In fact in $d=2$ even for $4<\alpha<\infty$ the system exhibits quasi-long-range order at low temperature, so that our treatment remains justified.. Under these circumstances, a description of the Hamiltonian in terms of harmonic fluctuations around the mean-field ground state – namely LSW theory – is very well justified, and its quantitative validity can be monitored systematically during the time evolution. Moreover, as discussed in the Supplementary Material (SM - SuppMat ) the time-reversal symmetry of the $\alpha$XX Hamiltonian and of the initial state leads to the remarkable property that the forward- and backward-evolved states $|\Psi(t)\rangle$ and $|\Psi(-t)\rangle$ give the same expectation values for time-reversal symmetric observables as well as for the EE: as a consequence all of our results apply as well to the evolution driven by the antiferromagnetic Hamiltonian ($J<0$), for which the initial state we choose is a superposition of highly excited eigenstates.

The central output of LSW theory paper1-XX_LR ; SuppMat is the reduction of the Hamiltonian Eq. (1) to that of free bosonic QP, ${\cal H}\approx\sum_{\bm{k}}\omega_{\bm{k}}\beta_{\bm{k}}^{\dagger}\beta_{\bm{k}}$, with the dispersion relation

[TABLE]

where $\gamma_{k}=\sum_{\bm{r}\neq 0}e^{i{\bm{k}\cdot\bm{r}}}|{\bm{r}}|^{-\alpha}$ is the Fourier transform of the interaction. The long-wavelength behavior of this dispersion relation is strongly affected by the value of the $\alpha$ exponent: indeed for small $k$ one finds paper1-XX_LR that $\omega_{\bm{k}}\approx c(\alpha)k^{z}$, where the exponent $z$ has value $z=1$ for $\alpha\geq d+2$ – continuously connecting with the well-known case of nearest-neighbor interactions Mattis-book – while $z=(\alpha-d)/2<1$ for $d\leq\alpha\leq d+2$, and $z=0$ for $\alpha\leq d$. When $z=1$ the QP excitations have a well-defined, $\alpha$-dependent maximum group velocity $v_{\rm max}=c$, which also sets a characteristic time scale $t^{*}(l)=l/(2c)$ for relaxation of observables in a subsystem of linear size $l$, in accordance with the universal SR picture of post-quench dynamics. On the other hand, the regime with $z<1$ is characterized by a divergent group velocity as $k\to 0$, $v_{\rm max}\sim k^{z-1}$, suggesting the breakdown of the SR picture. This result is on par with the field-theoretical (i.e. long-wavelength) description of the Hamiltonian Eq. 1, corresponding to a 2-component $(d+1)$-dimensional critical field theory with dynamical critical exponent $z$ SuppMat . In the following we shall particularly focus on the intermediate regime $d\leq\alpha\leq d+2$, exemplified by the experimentally relevant case of dipolar interactions ($\alpha=3$) in $d=2$. 555The regime $\alpha\leq d$ is continuously connected to the case of infinite-range interactions, describing the dynamics of the collective spin – we shall postpone the study of this regime to future work.

From straight to curved light cones. A fundamental trait of the post-quench dynamics is represented by the characteristic buildup of spin-spin correlations starting from the uncorrelated, factorized state – in the following we shall concentrate on the correlation function for the $z$ spin component $C^{zz}(\bm{r},t)=\langle\delta S_{i}^{z}\delta S_{i+\bm{r}}^{z}\rangle(t)$ where $\delta O=O-\langle O\rangle$; (see SuppMat for the details of its LSW calculation). In the case of $\alpha>d+2$ (so that $z=1$) our LSW results in $d=2$, as well as the DMRG ones in $d=1$, show a ballistic propagation of the correlation front – in the LSW calculation the speed of propagation is consistent with 2$c(\alpha)$. When instead the dynamical exponent $z$ becomes different from 1, the long-wavelength dynamics of the system loses the Lorentz invariance characteristic of the $z=1$ case, and time and space enter very differently in the low-energy description. This insight suggests that the distance over which correlations establish at time $t$ advances as $r\sim t^{1/z}$, namely super-ballistically, as a result of the divergence of the group velocity of the fastest QP. This is precisely what is observed in both the $d=2$ LSW data [see Fig. 2(a,c)] and the $d=1$ DMRG data [Fig. 2(b,d)]: the space-time correlation fronts acquire a finite curvature which is quantitatively explained by the non-linear relationship between space and time established by the LR interactions. The centrality of the $z$ exponent in determining the shape of correlation fronts has been very recently established for the linear response to local quenches in the field theoretical description of the critical $1d$ Ising model with LR interactions, and of $1d$ free fermions with LR hopping and pairing maghrebietal2016 . Our results extend the role of the $z$ exponent to the case of global quenches beyond linear response, and within the microscopic description of bosonic critical and LR-interacting models.

Multi-speed evolution of finite-$k$ fluctuations. The appearance of a curved space-time correlation front, while revealing the existence of a divergent group velocity for QP, remains silent about the fact that this group velocity spans the whole interval from $\infty$ to 0 when moving from the center to the edge of the Brillouin zone. But a deeper analysis of the spatial structure of correlations inside the curved light cone can capture the extreme dispersion of group velocities of QPs. The guiding idea is that the group velocity governs the evolution of wave packets – namely of superpositions of Fourier modes with frequencies tightly distributed around a reference value. For the quench under examination, LSW predicts that the evolution of the spin structure factor $S^{zz}(\bm{k})=N^{-1}\sum_{ij\in N}e^{i{\bm{k}}\cdot(\bm{r}_{i}-\bm{r}_{j})}\langle\delta S_{i}^{z}\delta S_{j}^{z}\rangle_{t}$ contains a constant part plus a term oscillating at frequency $2\omega_{k}$

[TABLE]

where $g^{(1)}_{\bm{k}}=s\gamma_{\bm{k}}/4\gamma_{0}$ and $g^{(0)}_{\bm{k}}=s/2-g^{(1)}_{\bm{k}}$. In order to superpose these Fourier modes into a “wave packet” one can consider the evolution of the structure factor on a subsystem $A$ of linear size $l_{A}$ (obtained by replacing $N^{-1}\sum_{ij\in N}\to l_{A}^{-d}\sum_{ij\in A}$ in the above expression), which takes the form of a convolution sum in Fourier space

[TABLE]

Here $f_{A}$ is the Fourier transform of the characteristic function of subsystem $A$ SuppMat , which, for $l_{A}\gg 1$, has support roughly corresponding to the tight interval $[-\pi/l_{A},\pi/l_{A}]^{\otimes d}$. As a result the subsystem structure factor takes the approximate form

[TABLE]

where $t_{\bm{k},a}^{*}=l_{A}/[2v_{G,a}(\bm{k})]$ and $v_{G,a}(\bm{k})=\partial_{k_{a}}\omega_{\bm{k}}$ is the $a$-th component of the group velocity, and ${\rm sinc}(x)=\sin(x)/x$. Hence the subsystem structure factor exhibits the characteristic time dependence of a wave packet, with oscillations at frequency $2\omega_{\bm{k}}$ contained within an envelope that decays over a characteristic time scale $t_{\bm{k}}^{*}=\min_{a=1\dots d}t_{\bm{k},a}^{*}$. Otherwise stated, $S_{A}^{zz}(\bm{k},t)$ is a sum of close-by frequencies which dephase over a characteristic time scale $t_{\bm{k}}^{*}$, directly related to the group velocity. As shown in the SM SuppMat , Eq. (5) provides an excellent approximation to the LSW data for $l_{A}\sim{\cal O}(10)$. This implies that the subsystem structure factor exhibits a peculiar form of multispeed pre-thermalization, in which each Fourier component relaxes over a strongly ${\bm{k}}$-dependent time scale, which in turn is linear in the subsystem size for finite $\bm{k}$ – see Fig. 3(a). On the other hand for $k=0$ the group velocity diverges, which for a subsystem of size $l_{A}$ implies that $v_{G}(0)\sim c~{}l_{A}^{1-z}$. As a consequence the $k=0$ component of the structure factor, related to the fluctuations of the uniform subsystem magnetization as $\langle\delta^{2}S_{A}^{z}\rangle/l_{A}^{d}$, relaxes over a characteristic time scale $t_{0}^{*}\sim l_{A}^{z}/(2c)$, scaling sub-linearly with subsystem size. This is consistent with the super-ballistic propagation of correlations – of which $\langle\delta^{2}S_{A}^{z}\rangle$ simply represents the integral over a subsystem – and it is clearly exemplified by the data of Fig. 3(b) for a 2$d$ system with dipolar interactions. Furthermore, following this line of thinking one can introduce the Fourier transform of the spin on a subsystem $A$, $S^{z}_{\bm{k},A}=l_{A}^{-d/2}\sum_{i\in A}e^{i{\bm{k}}\cdot{\bm{r}}_{i}}S_{i}^{z}$, and the ${\bm{k}}$-dependent correlation function $C^{zz}_{\bm{k}}(A_{i},A_{j})=\langle\delta S^{z}_{\bm{k},A_{i}}\delta S^{z}_{-\bm{k},A_{j}}\rangle$, where $A_{i}$ and $A_{j}$ are two subsystems at distance ${\bm{r}}_{ij}$ – so that $S_{A}^{zz}(\bm{k})=\sum_{A_{i},A_{j}\in A}C^{zz}_{\bm{k}}(A_{i},A_{j})$, where subsystems $A_{i}$’s tile regularly (namely without overlaps) the subsystem $A$. As shown in the SM SuppMat , the $\bm{k}$-dependent correlation function $C^{zz}_{\bm{k}}$ after the quench exhibits a correlation front which advances at a speed $2v_{G,a}(\bm{k})$ in the $a$-th direction (or a ${\bm{k}}$-dependent light cone, by virtue of the strong group velocity dispersion), and consistently its integral on a region of linear size $l_{A}$ relaxes on the above-cited characteristic time scales $t_{\bm{k}}^{*}$.

Entanglement dynamics. Finally we turn to the dynamics of the entanglement entropy (EE), defined as $E_{A}=-{\rm Tr}\rho_{A}\log\rho_{A}$, where $\rho_{A}$ is the reduced density matrix of subsystem $A$. The EE is sensitive to the relaxation of all forms of fluctuations, and therefore, due to the extreme dispersion in relaxation timescales for the Fourier components of the fluctuations, it can be expected to exhibit a much slower relaxation than long-wavelength fluctuations. Fig. 3(c) shows indeed that the EE obeys a scaling form $E_{A}(t)=l_{A}^{d}f_{E}(t/l_{A})$, exhibiting a characteristic time scale which grows linearly with $l_{A}$; its detailed time dependence shows a fast, sublinear rise at short times, followed by significant slowdown resulting in a very broad shoulder (the broader the smaller $\alpha$). The EE dynamics can be quantitatively understood by conjecturing that it stems from the relaxation of independent QP modes, each contributing an entropy $s_{\rm eq}(\bm{k})=(1+n_{\bm{k}})\log(1+n_{\bm{k}})-n_{\bm{k}}\log n_{\bm{k}}$ in the stationary state (generalized Gibbs ensemble, GGE VidmarR2016 ), where $n_{\bm{k}}=\langle\beta^{\dagger}_{\bm{k}}\beta_{\bm{k}}\rangle$ is the conserved number of QP in each mode fixed by the initial state. One can then postulate the following form

[TABLE]

where the sum is restricted to modes ${\bm{k}}$ whose wavelength is commensurate with subsystem $A$ (a similar form can be proved to be exact for quenches in one-dimensional free fermions fagottiC2008 ). A simple choice of the $f$ function, $f(x)=\min(x,1)$, produces a very convincing agreement with the LSW data, as shown in Fig. 3(c). This result finally reconciles the slow EE dynamics with the ultrafast dynamics of correlations and long-wavelength fluctuations in LR interacting systems, by exposing the dependence of the EE on the slow evolution of short-wavelength fluctuations.

Conclusions and outlook. Based on linear spin-wave and density-matrix renormalization group calculations, we have unveiled a scenario of multi-speed prethermalization for lattice spin models with long-range (LR) interactions. Our results reconcile the seemingly contradictory aspects of very fast evolution of correlations and long-wavelength fluctuations on a subsystem on the one hand, and of the very slow buildup of entanglement entropies on the other hand. This multi-speed scenario, obtained explicitly for the LR interacting XX model of magnetism, can be easily translated to other models of LR interactions studied in the recent past (especially the transverse field Ising model haukeT2013 ; schachenmayeretal2013 ; eisertetal2013 ; hazzardetal2014 ; foss-feigetal2015 ; schachenmayeretal2015 ; cevolanietal2015 ; maghrebietal2016 ), allowing one to reconcile the discrepancy between the correlation dynamics and entanglement dynamics already observed in those systems. In order to reconstruct the variety of characteristic time scales governing the dynamics we offer a detailed diagnostic tool (namely the analysis of the subsystem structure factor) – this is of immediate interest to current experiments on quantum gas microscopes, based on trapped ions richermeetal2014 ; Jurcevic2014 or neutral atoms cheneauetal2012 ; Labuhnetal2016 , which have direct access to the diagnostics in question. As an example, the multi-speed dynamics is realized by dipolar interactions in $d=2$, which are implemented by magnetic atoms dePazetal2013 and resonantly interacting Rydberg atoms Browaeysetal2016 . Hence our results pave the way for the reconstruction of most salient features of LR interaction dynamics to be searched in new-generation microscope experiments, as well as in systems beyond atomic physics (such as nano-engineered magnets Khajetooriansetal2015 , ensembles of nitrogen-vacancy centers Kucskoetal2016 and nuclear spins Alvarezetal2015 ).

Acknowledgments We thank A. Browaeys, G. Carleo, and M. Kastner for fruitful discussions and M. Fagotti for correspondence. This work is supported by ANR (“ArtiQ” project).

Appendix A Theorem on time-reversal symmetry of unitary evolution

A.1 General theorem

In this section we prove the following general

Theorem. Be ${\cal H}$ an Hamiltonian admitting a representation as a real-valued matrix on a given basis $|\phi\rangle$ of Hilbert space, and $|\Psi(0)\rangle$ an initial state which has all real coefficients when decomposed on the same basis. Then, for all observables $A$ which are represented as real-valued matrices on the basis, the evolution of the expectation value $\langle A\rangle(t)=\langle\Psi(t)|A|\Psi(t)\rangle$ is invariant under time reversal, namely

[TABLE]

Moreover, if the basis $|\phi\rangle$ is factorized among subsystems $A$ and $B$, $|\phi\rangle=|\phi^{(A)}\rangle\otimes|\phi^{(B)}\rangle$, then

[TABLE]

where $S_{A}$ is the entanglement entropy of subsystem $A$.

Proof. The proof of the above result is readily achieved by inserting resolutions of the identities on the basis $|\psi_{a}\rangle$ of eigenstates of ${\cal H}$ with eigenenergies $E_{a}$:

[TABLE]

where we have introduced the frequencies $\omega_{ab}=(E_{a}-E_{b})/\hbar$, and the decomposition of the initial state on the Hamiltonian eigenbasis:

[TABLE]

with $\Psi_{a}=\langle\psi_{a}|\Psi(0)\rangle$. The matrix elements and scalar products entering in Eq. (9) are real by hypothesis (because they can always be reduced to products of real coefficients by inserting resolutions of the identities on the $|\phi\rangle$ basis). Hence the product $\langle\psi_{a}|A|\psi_{b}\rangle\langle\Psi_{a}\Psi_{b}\rangle$ is symmetric under the exchange of the $a$ and $b$ indices, and therefore the $\sin(\omega_{ab}t)$ term vanishes, leaving out only the term $\cos(\omega_{ab}t)$ term, which is even in $t$.

Given that $|\Psi(t)\rangle=\sum_{a}\Psi_{a}(t)|\psi_{a}\rangle$ with $\Psi_{a}(t)=\exp(-iE_{a}t/\hbar)\Psi_{a}$ and $\Psi_{a}\in\mathbb{R}$, one has $\Psi_{a}(-t)=\Psi^{*}_{a}(t)$ – namely time reversal amounts to simple conjugation of the coefficients of the decomposition on the Hamiltonian eigenbasis, or, equivalently, to conjugation of the coefficients of the decomposition on the basis $|\phi\rangle$. Consider then the Schmidt decomposition of the evolved state

[TABLE]

where $\lambda_{p}\in\mathbb{R}$. If the basis $|\phi\rangle$ admits the factorized form $|\phi_{A}\rangle\otimes|\phi_{B}\rangle$, then the basis of the Schmidt decomposition $|a_{p}(t)\rangle$ (resp. $|b_{p}(t)\rangle$) can be decomposed on the basis $|\phi_{A}\rangle$ (resp. $|\phi_{B}\rangle$) as $|a_{p}(t)\rangle=\sum_{l}\varphi^{(A)}_{p,l}(t)|\phi_{l}^{(A)}\rangle$ (resp. $|b_{p}(t)\rangle=\sum_{l}\varphi^{(B)}_{p,l}(t)|\phi_{l}^{(B)}\rangle$). The fact that $|\Psi(-t)\rangle$ is related to $|\Psi(t)\rangle$ by conjugation of the coefficients of the decomposition on the $|\phi\rangle$ basis implies that $\varphi^{(A)}_{p,l}(-t)=(\varphi^{(A)}_{p,l})^{*}(t)$ and $\varphi^{(B)}_{p,l}(-t)=(\varphi^{(B)}_{p,l})^{*}(t)$, but that $\lambda_{p}(-t)=\lambda_{p}(t)$. Given that the entanglement entropy only depends on the Schmidt coefficients $\lambda_{p}$ as $S_{A}=-\sum_{p}\lambda_{p}\log\lambda_{p}$, we conclude that it is invariant under time reversal.

A.2 Discussion

The hypotheses of the above theorem are satisfied by the initial state $|\Psi(0)\rangle$ we choose in this work, as well as by the $\alpha$XX Hamiltonian ${\cal H}_{\alpha{\rm XX}}$, due to its time-reversal symmetry - the basis $|\phi\rangle$ is just provided by the computational basis of eigenstates of the $S_{i}^{z}$ operators. All the observables of interest to this work - namely correlation functions - also satisfy the hypothesis of real valuedness on the same basis. Finally the computational basis is fully factorized, namely it is factorized for any spatial $A$-$B$ bipartition.

The above result amounts to say that the Hamiltonians ${\cal H}_{\alpha{\rm XX}}$ and $-{\cal H}_{\alpha{\rm XX}}$ generate the same dynamics for the observables of interest to this work starting from the chosen state $|\Psi(0)\rangle$. This is a remarkable fact, given that the two Hamiltonians have quite different low-energy properties: ${\cal H}_{\alpha{\rm XX}}$ is ferromagnetic and not frustrated, while $-{\cal H}_{\alpha{\rm XX}}$ is antiferromagnetic and strongly frustrated due to the long-range nature of the couplings.

Appendix B Details on the spin-wave approximation

In this section we provide additional details on the semi-classical approximation used to study the dynamics of the $\alpha$-XX model of Eq. (1) of the main text.

B.1 Assumptions behind spin-wave theory

The quench dynamics we studied starts from the mean-field ground state of the $\alpha$-XX model, $|\Psi(0)\rangle=|\uparrow_{x}\rangle^{\otimes N}$. Under the mean-field approximation $S_{i}^{\beta}S_{j}^{\beta}\approx S_{i}^{\beta}\langle S_{j}^{\beta}\rangle+\langle S_{i}^{\beta}\rangle S_{j}^{\beta}-\langle S_{i}^{\beta}\rangle\langle S_{j}^{\beta}\rangle$, the chosen initial state does not evolve, because $\langle S_{j}^{y}\rangle=0$, so that the mean-field Hamiltonian is proportional to the $S_{i}^{x}$ operator. Our semi-classical treatment of the ensuing dynamics (based on linear spin-wave - LSW - theory) goes beyond the mean-field approximation, but it relies on the hypothesis that the full many-body state remains close to the initial state throughout the evolution – in other words, the quantum fluctuations building up on top of the mean-field state are weak at all times. Since the initial state breaks the rotational symmetry of the system - and the evolved state inherits the symmetry breaking - this approach tacitly implies that symmetry breaking is observed at sufficiently low energy in the system, corresponding to the excitation energy of the mean-field ground state with respect to the actual ground state of the system. We shall explicitly test this assumption below.

Finally we observe that any finite spin system, although prepared in an initial state that breaks the symmetry of the Hamiltonian, is expected to undergo a depolarization dynamics which restores this symmetry, and this aspect is obviously not captured by our calculations. Nonetheless the explicit object of our study is not the evolution of the average collective spin - which is absent within our approach - but rather the evolution of the fluctuations around the instantaneous average value, which is therefore tacitly assumed to be essentially independent of the average value; or alternatively it is assumed to occur on much faster timescales than the dynamics of the collective spin. This latter assumption is well justified in the limit $N\to\infty$, as the collective spin relaxes over time scales that grow linearly with system size – reflecting a characteristic energy scale of order ${\cal O}(1/N)$ for the collective-spin excitations, related to the so-called Anderson’s tower of states [37]. On the other hand, the frequencies of the LSW excitations do not scale with system size, so that the dynamics of fluctuations around the average value can completely decouple from that of the average in the $N\to\infty$ limit.

B.2 Details of the spin-wave approach

To describe weak quantum fluctuations around the mean-field state one can introduce the linearized Holstein-Primakoff (HP) transformation, mapping spins onto bosons, and describing small deviations around the mean-field state:

[TABLE]

The linearization of the HP transformation to treat the dynamics of the system requires therefore that $\langle b_{i}^{\dagger}b_{i}\rangle(t)\ll 2s$ at all times. To study the dynamics we rewrite the Hamiltonian in terms of the HP bosons dropping all terms beyond quadratic in the $b^{(\dagger)}_{i}$’s operators, and obtaining (up to a constant energy that we neglect):

[TABLE]

where

[TABLE]

Assuming translational invariance we introduce the Fourier transformed Bose operators $b_{\bm{k}}=N^{-1/2}\sum_{j}e^{-i{\bm{k}}\cdot{\bm{r}}_{j}}b_{j}$, where $N$ is the number of spins, together with the Fourier transform of the interaction $\gamma_{\bm{k}}=\sum_{j}e^{i{\bm{k}}\cdot{\bm{r}}_{j}}r_{0j}^{-\alpha}$. With these notations, the Hamiltonian in momentum space takes the form:

[TABLE]

with $A_{\bm{k}}=2sJ(\gamma_{0}-\gamma_{\bm{k}}/2)$ and $B_{\bm{k}}=sJ\gamma_{\bm{k}}$. The above Hamiltonian is then diagonalized by a Bogoliubov transformation [38] : $b_{\bm{k}}=u_{\bm{k}}\beta_{\bm{k}}-v_{\bm{k}}\beta_{-\bm{k}}^{\dagger}$, with Bogoliubov coefficients $u_{\bm{k}}={\rm sign}(b_{\bm{k}})[(A_{\bm{k}}/\omega_{\bm{k}}+1)/2]^{1/2}$ and $v_{\bm{k}}=\sqrt{u_{\bm{k}}^{2}-1}$, where we have introduced $\omega_{\bm{k}}=\sqrt{A_{\bm{k}}^{2}-B_{\bm{k}}^{2}}=2sJ\gamma_{0}\sqrt{1-\frac{\gamma_{\bm{k}}}{\gamma_{0}}}$. Since $\omega_{k=0}=0$, the symmetry-restoring $k=0$ mode cannot be treated on the same footing as the others. We have simply discarded the $k=0$ sector in all our calculations, a procedure which does not affect the results in the thermodynamic limit (an alternative procedure would correspond to gap the $k=0$ mode with a term which vanishes in the thermodynamic limit, providing similar results). Finally, we obtain an effective Hamiltonian (dropping a constant ground-state energy):

[TABLE]

To obtain the physical properties at a given time $t$, we work in Heisenberg representation:

[TABLE]

where on the last line we used $u_{\bm{k}}^{2}=(1/2)(A_{\bm{k}}/\omega_{\bm{k}}+1)$, $v_{\bm{k}}^{2}=(1/2)(A_{\bm{k}}/\omega_{\bm{k}}-1)$ and $2u_{\bm{k}}v_{\bm{k}}=B_{\bm{k}}/\omega_{\bm{k}}$. Since the initial state is the vacuum of the $b$’s operators, it trivially satisfies Wick’s theorem [38] for the $b^{(\dagger)}(0)$ operators. As $b^{(\dagger)}(t)$’s operators are linear combinations of the $b^{(\dagger)}(0)$’s, Wick’s theorem is then satisfied at any time. As a consequence, in a translationally invariant system all the information required to reconstruct the physical properties is contained in the one-body correlators:

[TABLE]

B.3 Maximum group velocity for $d+2\lesssim\alpha$

Here we calculate explicitly the divergence of the group velocity around $k=0$ when $\alpha$ approaches $d+2$ from above.

[TABLE]

where $kx={\bm{k}}\cdot{\bm{r}}$, and the lattice sum has been approximated with an integral. The prefactor of the $1/\sqrt{\alpha-d-2}$ divergence depends on the dimension $d$. The above calculation shows the expected divergence for $\alpha\to(d+2)^{+}$.

B.4 Validity condition of spin-wave theory

The LSW approach to the dynamics is justified under the condition

[TABLE]

Since $\langle b_{\bm{k}}^{\dagger}b_{\bm{k}}\rangle(t)$ does not depend on $s$, one immediately sees that the spin-wave approximation becomes exact in the limit of large $s$.

The factor $r(t)$ is generally found to grow with $t$ (see Fig. 4(a)), so that the worst-case scenario to test the consistency of LSW theory can be found in the limit $t\to\infty$. In this limit, as well as in the $N\to\infty$ one, the $r$ factor takes the form:

[TABLE]

where we have replaced $\sin^{2}(\omega_{\bm{k}}t)$ by $1/2$, since oscillations at nearby frequencies are dephased at long times and only the average survives. The momentum integration runs over the Brillouin zone (BZ).

Fig. 4(b) shows $r(\infty)$ for finite-size lattices as a function of $\alpha$ in dimensions $d=1,2,3$, and in the extreme quantum case $s=1/2$. To put these results in perspective one can notice that:

for $\alpha<d$, $B_{\bm{k}}/\omega_{\bm{k}}=(\gamma_{\bm{k}}/2\gamma_{0})/\sqrt{1-\gamma_{\bm{k}}/\gamma_{0}}$ goes to zero when $N\to\infty$ at any nonzero $k$, since $\gamma_{0}\to\infty$ while $\gamma_{\bm{k}}$ is finite. On the other hand, $B_{\bm{k}}/\omega_{\bm{k}}$ has a finite limit for $k=0$. As a consequence $r(\infty)$ vanishes for $\alpha<d$, corresponding to a mean-field regime, identical (in the thermodynamic limit) to the $\alpha=0$ limit, where the dynamics occurs only in the collective spin $k=0$ sector; 2. 2.

for $\alpha>d$, if the exact value of the integral in Eq. (23) must be evaluated numerically, one can easily identify the condition under which it is convergent in the infrared, validating the LSW approach. Since $B_{k=0}$ is finite, the condition requires $1/\omega_{\bm{k}}^{2}$ to grow at small $k$ more slowly than $k^{-d}$. If $\omega_{\bm{k}}\sim k^{z}$ with $z$ the dynamical exponent, the convergence requires $z<d/2$. As discussed in the main text, $z=1$ for $\alpha\geq d+2$ and $z=(\alpha-d)/2$ for $d\leq\alpha\leq d+2$ [26]. The convergence is hence always guaranteed in $d=3$, in $d=2$ for $\alpha<4$ (with a logarithmic divergence for $\alpha\geq 4$), and in $d=1$ for $\alpha<2$. Due to the weak logarithmic divergence of $r(\infty)$, LSW remains quantitative on sufficiently small sizes for $\alpha=2$ in $d=1$, and for $\alpha>4$ in $d=2$.

We notice that the above conditions of validity of LSW theory in the pre-thermalized state are more restrictive than the condition that the theory has to satisfy at zero temperature - which reads $z<d$, always satisfied in $d=2$ and 3, and satisfied for $\alpha<3$ in $d=1$. On the other hand, the condition $z<d/2$ is the same as the one guaranteeing the validity of LSW theory at finite temperature $T$ – which requires the thermal correction to the order parameter

[TABLE]

to be finite. Even though the pre-thermalized state is not a thermal one, it shares with a thermal state the property of being at finite energy density.

B.5 Correlation functions and entanglement in LSW theory

Here we provide more details on the LSW expression for the various correlation functions studied in this work. Making use of Eqs. (19) and (12), the structure factor for the $z$ component of the magnetization takes the form:

[TABLE]

which corresponds the expression given in Eq. 3 of the main text. A similar calculation for the $yy$ correlations gives:

[TABLE]

The calculation of the entanglement entropy within LSW theory is well established, and based on the evaluation of the correlation matrix – we refer the reader to Ref. [39] for the details.

Appendix C Quantum field-theory description of the long-range XX model

The XX model is a fundamental model of strongly interacting bosons with translational as well as particle-hole symmetry. As it displays a gapless spectrum and algebraic correlations, its long-wavelength properties are expected to be correctly described via a quantum field theory in continuum space. In the case of short-range interactions the quantum-field theory corresponds to the quantum O(2) model for a complex scalar field $\psi({\bm{r}},\tau)$, whose action reads

[TABLE]

The above action can be obtained explicitly [40, 41] by mapping the quantum spins onto bosons via an HP transformation (with quantization axis along $z$) giving rise to a Bose-Hubbard-like model; using the coherent-state path-integral representation for the partition function of the bosonic model; and applying a Hubbard-Stratonovich transformation via the introduction of the auxiliary field $\psi$. The case of long-range interactions decaying as $1/r^{\alpha}$ in the spin model corresponds to long-range hopping in the bosonic model. The Hubbard-Stratonovich transformation leads to an effective action for the auxiliary field with a quadratic part given by [41]

[TABLE]

where

[TABLE]

where $K$ is a positive constant such that $K+\gamma_{\bm{k}}>0$, namely $K>\max_{\bm{k}}|\gamma_{\bm{k}}|$, and $\gamma_{\bm{k}}$ is the Fourier transform of the $1/r^{\alpha}$ interaction. A Taylor expansion of the numerator gives then

[TABLE]

In the last term the first part (containing a single spatial integral) reproduces the short-range part of the coupling for the auxiliary field, while the double integral captures the long-range part. Therefore, upon proper rescaling of the coefficients, the inclusion of long-range interactions amounts to adding the $v$ term to the effective action of Eq. (27).

Moving to $k$-space and Matsubara frequencies, the action of the long-range interacting system takes the form [42]

[TABLE]

where ${\cal A}_{\rm int}$ contains the terms beyond quadratic, and $\sigma=\alpha-d$. Far from its quantum critical point, the action can be treated within the Gaussian approximation [43] (which represents the field-theory analog to the LSW approach we relied upon); in this case a RG treatment by simple dimensional analysis shows that the long-range interaction term $v^{\prime}q^{\sigma}$ is relevant (and the $q^{2}$ term irrelevant) when $\sigma<2$, namely for $\alpha<d+2$. Under this circumstance the quantum field theory loses therefore Lorentz invariance, and $z=\sigma/2<1$ acquires the role of new dynamical critical exponent.

As it is not our purpose to investigate the quantum critical point of the long-range quantum O(2) model, we can stop our analysis at the level of the Gaussian approximation. In the case of the long-range quantum O(1) model – describing the quantum phase transition of the long-range Ising model in a transverse field – a perturbative RG analysis can be found in Ref. [16].

Appendix D Pre-thermalization dynamics in a subsystem

The LSW approach, neglecting interactions among the spin waves, is unable to describe a true thermalization process. This is explicitly manifested in the expressions in Eq. (19), (25) and (26): the different Fourier components of the structure factor show persistent oscillations at a frequency $2\omega_{\bm{k}}$ without any form of relaxation. Nonetheless, if one looks at the structure factor in a subsystem, the dephasing of oscillations at nearby frequencies leads to a *pre-*thermalization process, driving the subsystem towards a generalized Gibbs ensemble [31], with a different (generalized) “temperature” for each $k$ component. It is the purpose of this section to provide additional details on this pre-thermalization dynamics.

D.1 Relaxation of the fluctuations in a subsystem

Considering a subsystem $A$, it is convenient to introduce a weight function $w_{A}(i)$, weighting sites $i$ belonging to $A$. If $A$ contains $N_{A}$ sites, the most obvious choice is $w_{A}(i)=1/\sqrt{N_{A}}$ if $i\in A$ and [math] otherwise, but the following expressions are more general. The structure factor for the fluctuations of the spin component $S^{\beta}$ in the subsystem $A$ then reads :

[TABLE]

where we have introduced the form factor

[TABLE]

Typically, if $A$ has linear size $l_{A}$, $f_{A}(\bm{q})$ selects wavectors in the range $k\pm\Delta k$ with $\Delta k\approx\pi/l_{A}$.

Let us first consider the case $d=1$, postponing the case $d>1$ to the next subsection. If $w_{A}(i)$ selects sites inside a segment of length $l_{A}$, one has

[TABLE]

and $f_{A}(k=0)=l_{A}$. The relaxation of the fluctuations at wavevector $k$ in a subsystem then results from the dephasing of nearby frequencies around the “carrier” $\omega_{k}$. One expects such a dephasing mechanism to be effective after a timescale given by the inverse bandwidth of the frequencies involved in the subsystem dynamics, namely

[TABLE]

with $v_{G}(k)=d\omega_{k}/dk$ the group velocity. One thus recovers the typical timescale for the damping of subsystem fluctuations at wavevector $k$ as discussed in the main text.

More quantitatively, in the limit $N\to\infty$, one has

[TABLE]

For a sufficiently large $A$ subsystem, $f_{A}(q)$ is peaked around $q=0$, with $\int_{\rm BZ}(dq/2\pi)f_{A}(q)=1$. It is therefore legitimate to approximate $f_{A}(q)$ by a simple step function $f_{A}(q)=l_{A}$ for $-\pi/l_{A}<q<\pi/l_{A}$ and $f_{A}(q)=0$ otherwise, as illustrated on Fig. 5.

Within this approximation, one obtains

[TABLE]

From Eqs. (25) and (26), one sees that the structure factor $S^{\beta\beta}(k)$ has the form

[TABLE]

E.g. for $\beta=z$ we obtain:

[TABLE]

As the strongest dependence of $S^{\beta\beta}(k+q,t)$ on $q$ comes from the cosine term, in the integral over $q$ we can approximate all other $q-$dependent terms with their value at $q=0$. This reduces the calculation of Eq. (37) to the following integral:

[TABLE]

This leads to Eq. (5) of the main text (for the case $d=1$). We observe that the damping term ${\rm sinc}(\pi t/t_{k}^{*})$ comes from the dephasing of nearby frequencies, leading to an effective pre-thermalization within the subsystem.

As illustrated in Fig. 6, Eq. (5) of the main text reproduces perfectly the relaxation dynamics up to $t=t^{*}_{k}$. At $t>t^{*}_{k}$, the approximate expression of Eq. (5) predicts small revivals of the oscillations which are not observed with the complete expression of Eq. (32).

D.2 Multi-speed relaxation for $d>1$

The above considerations generalize readily to the case $d>1$. Considering a (hyper-)tetragonal subsystem of size $l_{A,1}\times l_{A,2}\times\dots$, one obtains Eq. (5) of the main text, which we repeat here for completeness:

[TABLE]

with $t_{\bm{k},a}^{*}=l_{A,a}/(2\partial_{k_{a}}\omega_{\bm{k}})$. The sinc term leads therefore to damping of the oscillations on a characteristic time scale $t^{*}_{\bm{k}}=\min_{a}t_{\bm{k},a}^{*}$. An important qualitative difference with respect to the $1d$ case is that, even if the dispersion relation is linear at low frequency, $\omega_{\bm{k}}=ck$ (as it is found for $\alpha>d+2$), the group velocity has several components $v_{G,a}({\bm{k}})=\partial_{k_{a}}\omega_{\bm{k}}=ck_{a}/k$, and the minimal $t_{\bm{k},a}$ depends on the orientation of the ${\bm{k}}$ vector with respect to the coordinate axes of the subsystem. In the case of an hypercube $l_{A}^{d}$, one has

[TABLE]

which, depending on the orientation of ${\bm{k}}$, ranges between $l/2c$ and $(l/2c)\sqrt{d}$.

D.3 Scaling forms for the evolution of the collective spin on a subsystem

We have seen in the main text that the evolution of the subsystem structure factor at ${\bm{k}}=0$ for the $z$ spin component is compatible with the following Ansatz

[TABLE]

for the evolution of the variance of the collective spin component $S_{A}^{z}$ – namely, a volume-law scaling at all times (including the initial state).

On the other hand, the variance of $S_{A}^{y}$ is found to relax towards a super-extensive value $\propto l_{A}^{d+2z}$ – while it obeys a volume-law scaling in the initial state at $t=0$. At thermal equilibrium LSW predicts the same super-extensive scaling, which is therefore expected to hold for temperatures compatible with spontaneous symmetry breaking (namely, the hypothesis underlying LSW theory). The following scaling Ansatz

[TABLE]

regroups together this observation as well as the divergence of the group velocity, and it is found to be verified by the LSW data, as illustrated on Fig. 7 for $\alpha=3$ in $d=2$. Incidentally, this scaling form demands that $2z$ be smaller than $d$, as it is impossible that the collective spin fluctuations scale faster than the square of the volume, which the same criterium of stability of LSW approximation, required by the convergency of the $r_{\infty}$ in Eq. (22).

D.4 Mode-decomposition Ansatz and short-time behavior for the entanglement entropy

In this section we provide more details on the mode-decomposition Ansatz for the entanglement entropy, Eq. (6) of the main text. Its strict validity - even in the context of the LSW approximation - can only be assumed in the limit in which subsystem $A$ becomes very large, so that the subsystem excitations can be identified with those excitations of the bulk system whose wavelength is commensurate with the linear dimensions of the subsystem itself; and the populations of the subsystem modes can be deduced from those of the bulk excitations, as fixed by the initial conditions of the quench.

In order to minimize the finite-size effects of subsystem $A$, its ideal geometry is that of a slice of a $L_{1}\times L_{2}^{d-1}$ hyper-torus of size $l_{A}\times L_{2}^{d-1}$ ($l_{A}<L_{1}$). In this case finite-size limitations to the applicabiity of Eq. (6) only appear in the ”short” direction (of length $l_{A}$) which is perpendicular to the $A$-$B$ cut - this particular geometry is also very helpful in the numerical calculation of the entanglement entropy, as the correlation matrix necessary to reconstruct it [39] is block-diagonal in the wavevector longitudinal to the cut. In this geometry, for $d>1$ the time scale $t^{*}_{\bm{k}}$ associated to the buildup of the entropy in mode ${\bm{k}}$ is to be modified to $t^{*}_{\bm{k},\perp}=l_{A}/2|v_{\perp}|$, where $v_{\perp}=\partial_{k_{\perp}}\omega_{\bm{k}}$ is the component of the group velocity perpendicular to the $A$-$B$ cut; indeed only the propagation of the modes perpendicular to the cut is responsible for the entangling dynamics of the $A$ and $B$ subsystems. Figs. 8 and 9 confirm the validity of the decomposition Ansatz, as well as the choice of the scaling function $f$, for $d=1$ and $d=2$ respectively.

The linear dependence of the function $f(t/t^{*}_{\bm{k}})$ on time may naïvely suggest that the overall time dependence of the entanglement entropy at short times is always linear. This is generally true if the time scale $t^{*}_{\bm{k},\perp}$ admits a non-zero lower bound $t^{*}_{\rm min}>0$: for $t<t^{*}_{\rm min}$ all modes contribute a linearly growing entropy up to time $t^{*}_{\rm min}$. Such a condition is verified when $z=1$, namely for $\alpha>d+2$: there $t^{*}_{\rm min}=l_{A}/(2c)$. In this case, Eq. (6) of the main text shows indeed that the entropy grows linearly in time, with a prefactor proportional to the boundary of $A$, $S_{A}(t)\sim l_{A}^{d-1}t$. On the other hand when $\alpha<d+2$ $c(\alpha)\to\infty$ so that the lower bound on the time scales vanishes, $t^{*}_{\rm min}\to 0$. This means that for arbitrarily short times (in the thermodynamic limit) one can always find some modes ${\bm{k}}$ that contribute non-linearly to the entanglement entropy, as their contribution has already saturated (for $t>t_{\bm{k},\perp}$). Therefore in this regime one can expect a sublinear growth of entanglement.

The above expectations appear to be verified in Fig. 10. There we observe that the short-time evolution of the entanglement entropy in a $d=2$ system appears generally compatible with a power-law behavior $S_{A}\sim t^{x}$. If $\alpha>4$, a linear growth of entanglement entropy is clearly observed at short times; in particular, when rescaling time to $\tau=l_{A}/(2Js\gamma_{0})$ – where $2Js\gamma_{0}$ sets the bandwidth of the dispersion relation – the short time behavior becomes nearly independent of $\alpha$. On the other hand for $\alpha<4$ our data clearly indicates that $x<1$.

Appendix E Buildup of correlations among distant extended subsystems

In this section, we provide additional details concerning the space and time dependence of the ${\bm{k}}$-dependent correlation function $C^{zz}_{\bm{k}}(A_{i},A_{j})=\langle\delta S^{z}_{\bm{k},A_{i}}\delta S^{z}_{-\bm{k},A_{j}}\rangle$. The latter was introduced in the main text in order to probe the correlations among the ${\bm{k}}$ Fourier components of the magnetization profile on two equivalent subsystems $A_{i}$ and $A_{j}$ of linear size $l_{A}$, as a function of their spatial separation $r_{ij}$ – see Fig. 11 for a sketch.

Along the same line as the calculations discussed in Appendix D.1, one finds that

[TABLE]

where $f_{A}({\bm{q}})$ is the form factor of the two equivalent subsystems. Replacing the latter with a window function in the limit $l_{A}\gg 1$ (see Fig. 5), the momentum integration delivers

[TABLE]

The product of the ${\rm sinc}$ functions implies that $C^{zz}_{\bm{k}}(\bm{r},t)$ is peaked around $\bm{r}=\pm 2\bm{v}_{G}(\bm{k})t$, namely it exhibits correlation fronts moving ballistically in space at $\pm 2\bm{v}_{G}(\bm{k})$. This result agrees with the physical intuition that correlations at wavevector $\bm{k}$ between two subsystems appear when they are reached by a pair of quasiparticles of wavevector $\bm{k}$, emitted halfway between them, and counter-propagating ballistically at speeds $\pm\bm{v}_{G}(\bm{k})$.

This result shows that the dynamics of the ${\bm{k}}$-dependent correlations among two extended subsystems exhibits a multiple light-cone structure, related to the dispersion of group velocities. Each correlator $C^{zz}_{\bm{k}}(\bm{r},t)$ possesses a different light-cone dynamics, and the “causal cone” of correlations spreads only along the vector direction of the group velocity. All these aspects are captured by Fig. 12, showing moreover that, for $z<1$, the ${\bm{k}}$-dependent light cone becomes “curved” for $k=0$, revealing the nonlinear dispersion relation $\omega\sim k^{z}$; and moreover the correlation fronts acquire an approximate spherically symmetric structure in the long-time limit.

Appendix F Time dependent DMRG simulations

In this section we provide more details about the time-dependent simulations based on the density-matrix renormalization group (DMRG). The time-dependent DMRG [44, 45, 46, 47] is an efficient algorithm to study the static and dynamical properties of quantum systems in one dimension. Our calculations have been performed for a chain of $L=30$ sites, taking as initial state for the evolution the ground state of the fully connected Hamiltonian ($\alpha=0$):

[TABLE]

which has quantum numbers $S_{\rm tot}=N/2$ and $S^{z}_{\rm tot}=0$ for the collective spin operator. This state is a projection onto the $S^{z}_{\rm tot}=0$ sector of Hilbert space of the perfectly polarized state $|\uparrow_{x}\rangle^{\otimes N}$ used in the LSW calculations. This choice facilitates significantly the DMRG calculations, allowing us to restrict the time evolution to the $S^{z}=0$ sector. We then studied the unitary dynamics generated by the Hamiltonian ${\cal H}_{\alpha{\rm XX}}$ for several values of $\alpha$ in the range $[1.5,2.5]$. The time-evolution is based on a Runge-Kutta 4th order scheme, with a time step of $\delta t/J=0.05$. We used a variable number of states (up to 800) in order to keep the truncation error smaller than $10^{-5}$ at each time step.