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

TL;DR

This paper investigates universal scaling laws governing how correlations spread in quantum systems with both short- and long-range interactions, revealing a double-structured, velocity-dependent behavior after a quench.

Contribution

It introduces a unified framework for correlation spreading in quantum systems, deriving universal exponents and velocities applicable to various interaction ranges and excitation spectra.

Findings

Correlation edge spreads with velocity twice the maximum group velocity.

For long-range interactions, correlation edge spreads slower than ballistic.

Maxima of correlations propagate faster or ballistically depending on system gap.

Abstract

We study the spreading of information in a wide class of quantum systems, with variable-range interactions. We show that, after a quench, it generally features a double structure, whose scaling laws are related to a set of universal microscopic exponents that we determine. When the system supports excitations with a finite maximum velocity, the spreading shows a twofold ballistic behavior. While the correlation edge spreads with a velocity equal to twice the maximum group velocity, the dominant correlation maxima propagate with a different velocity that we derive. When the maximum group velocity diverges, as realizable with long-range interactions, the correlation edge features a slower-than-ballistic motion. The motion of the maxima is, instead, either faster-than-ballistic, for gapless systems, or ballistic, for gapped systems. The phenomenology that we unveil here provides a unified framework, which encompasses existing experimental observations with ultracold atoms and ions. It also paves the way to simple extensions of those experiments to observe the structures we describe in their full generality.

Figures

Equations

$$\hat{H}=\sum _{\mathbf{R}}h(\mathbf{R}){\hat{K}}_1(\mathbf{R})+\sum _{\mathbf{R},{\mathbf{R}}^{\prime }}J(\mathbf{R},{\mathbf{R}}^{\prime }){\hat{K}}_2(\mathbf{R},{\mathbf{R}}^{\prime }),$$

$$G\left(\mathbf{R},t\right)=g(\mathbf{R})-{\int }_{\mathcal{B}}\frac{d\mathbf{k}}{{\left(2\pi \right)}^D}\mathcal{F}\left(\mathbf{k}\right)\frac{{\text{e}}^{i\left(\mathbf{k}\cdot \mathbf{R}+2E_{\mathbf{k}}t\right)}+{\text{e}}^{i\left(\mathbf{k}\cdot \mathbf{R}-2E_{\mathbf{k}}t\right)}}{2},$$

$$2V_{\text{g}}(k_{\text{sp}})=\pm R/t,$$

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

$$G_{\text{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],$$

$$t^{\star }\propto R^{{\beta }_{\text{CE}}},{\beta }_{\text{CE}}=\chi /\gamma .$$

$$t_{\text{m}}\propto R^{{\beta }_{\text{m}}},{\beta }_{\text{m}}=z.$$

$$\hat{H}=-J\sum _{⟨\mathbf{R},{\mathbf{R}}^{\prime }⟩}\left({\hat{a}}_{\mathbf{R}}^{†}{\hat{a}}_{{\mathbf{R}}^{\prime }}+\text{H.c.}\right)+\frac{U}{2}\sum _{\mathbf{R}}{\hat{n}}_{\mathbf{R}}\left({\hat{n}}_{\mathbf{R}}-1\right),$$

$$E_k\simeq 2\sqrt{2J{\sin }^2\left(k/2\right)\left[2J{\sin }^2\left(k/2\right)+nU\right]}.$$

$$\mathcal{F}(k)=\frac{2n^{2}J{U}_{\text{f}}(U_{\text{f}}-U_{\text{i}})}{E_{k,\text{i}}{E}_{k,\text{f}}^2}{\sin }^2(k/2),$$

$$2E_k\simeq U-2(2n+1)J\cos (k).$$

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

$$\hat{H}=\sum _{\mathbf{R}\ne {\mathbf{R}}^{\prime }}\frac{J/2}{|\mathbf{R}-{\mathbf{R}}^{\prime }{|}^{\alpha }}\left[-\left({\hat{S}}_{\mathbf{R}}^x{\hat{S}}_{{\mathbf{R}}^{\prime }}^x+{\hat{S}}_{\mathbf{R}}^y{\hat{S}}_{{\mathbf{R}}^{\prime }}^y\right)+ϵ{\hat{S}}_{\mathbf{R}}^z{\hat{S}}_{{\mathbf{R}}^{\prime }}^z\right].$$

$$S_{\mathbf{R}}^x$$

$$S_{\mathbf{R}}^y$$

$$S_{\mathbf{R}}^z$$

$$E_k=\frac{J{P}_{\alpha }(0)}{2}\sqrt{1-\frac{P_{\alpha }(k)}{P_{\alpha }(0)}},$$

$$P_{\alpha }(k)\approx P_{\alpha }(0)+P_{\alpha }^{\prime }{k}^{\alpha -D},$$

$$\mathcal{F}\left(k\right)=\frac{{ϵ}_{\text{i}}}{8}\frac{P_{\alpha }\left(k\right)}{P_{\alpha }\left(0\right)}\sqrt{\frac{P_{\alpha }(0)-P_{\alpha }\left(k\right)}{P_{\alpha }(0)+{ϵ}_{\text{i}}{P}_{\alpha }\left(k\right)}}.$$

$$\hat{H}=\sum _{\mathbf{R}\ne {\mathbf{R}}^{\prime }}\frac{2J}{|\mathbf{R}-{\mathbf{R}}^{\prime }{|}^{\alpha }}{\hat{S}}_{\mathbf{R}}^x{\hat{S}}_{{\mathbf{R}}^{\prime }}^x-2h\sum _{\mathbf{R}}{\hat{S}}_{\mathbf{R}}^z.$$

$$S_{\mathbf{R}}^x$$

$$S_{\mathbf{R}}^y$$

$$S_{\mathbf{R}}^z$$

$$E_k=2\sqrt{h\left[h+J{P}_{\alpha }(k)\right]}.$$

$$E_k=\Delta +c|k{|}^z,$$

$$\mathcal{F}(k)=\frac{h\left(J_{\text{i}}-J_{\text{f}}\right)P_{\alpha }\left(k\right)}{8\left[h+J_{\text{f}}{P}_{\alpha }\left(k\right)\right]\sqrt{h\left[h+J_{\text{i}}{P}_{\alpha }(k)\right]}}.$$

Full text

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 lieb1972 . This bound implies ballistic propagation of equal time-correlation functions bravyi2006 that has been observed experimentally cheneau2012 ; fukuhara2013 and characterized numerically manmana2009 ; barmettler2012 ; carleo2014 ; dellanna2016 ; kormos2017 . Generalized LR bounds have been derived for long-range systems where the interactions decay algebraically, $1/R^{\alpha}$, with the distance $R$, see Ref. hastings2006, ; foss-feig2015, . The related experiments and numerical investigations have, however, lead to conflicting pictures hauke2013 ; eisert2013 ; jurcevic2014 ; richerme2014 ; cevolani2015 ; cevolani2016 ; buyskikh2016 . For instance, experiments on ion chains richerme2014 and numerical simulations within truncated Wigner approximation schachenmayer2015b 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 jurcevic2014 . 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. schachenmayer2013, ; hauke2013, ; eisert2013, ; cevolani2015, ; cevolani2016, ; buyskikh2016, .

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. kastner2015, ) 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 bendkowsky2009 ; weimer2010 ; schausz2012 ; browaeys2016 , nonlinear optical media firstenberg2013 , polar molecules micheli2006 ; yan2013 ; moses2017 , magnetic atoms griesmaier2005 ; beaufils2008 ; lu2011 ; baier2016 ; lahaye2009 , superconductors houck2012 , ion chains deng2005 ; islam2011 ; lanyon2011 ; schneider2012b ; NaturePhysicsInsight2012blatt , and solid-state defects childress2006 ; balasubramanian2009 ; dolde2013 .

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

[TABLE]

where $\mathbf{R}$ and $\mathbf{R}^{\prime}$ 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}\left(\mathbf{R},t\right)-G_{0}\left(\mathbf{R},0\right)$ with $G_{0}\left(\mathbf{R},t\right)\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 cheneau2012 ; fukuhara2013 ; richerme2014 ; jurcevic2014 . 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

[TABLE]

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. barmettler2012, ; natu2013, ; hauke2013, ; cevolani2015, ; cevolani2016, ; buyskikh2016, ; frerot2018, 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 haegeman2012 ; nakatani2014 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=\textrm{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,

[TABLE]

where $V_{\textrm{g}}=\nabla_{k}E_{k}$ is the group velocity. Since the latter is upper bounded by some value $V_{\textrm{g}}^{\star}$, Eq. (3) has a solution only for $R/t<2V_{\textrm{g}}^{\star}$. The correlation function then reads SeveralMax

[TABLE]

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

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^{\star}$, which move at $V_{\textrm{g}}^{\star}$, contribute to the correlations. There the maxima (m), defined by the equation $k^{\star}R-2E_{k^{\star}}t=\textrm{const}$, propagate at the velocity $V_{\textrm{m}}=2V_{\varphi}^{\star}\equiv 2E_{k^{\star}}/k^{\star}$, i.e. twice the phase velocity at the maximum of the group velocity, $k^{\star}$. 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^{\star}<\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_{\textrm{\tiny CE}}\simeq 2V_{\textrm{g}}^{\star}$ (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^{\star}$, $V_{\textrm{m}}\simeq 2V_{\varphi}^{\star}$ (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_{\textrm{\tiny CE}}$) and local maxima ($V_{\textrm{m}}$) on the one hand, to twice the group ($V_{\textrm{g}}^{\star}$) and phase ($V_{\varphi}^{\star}$) velocities at $k^{\star}$ on the other hand.

The distinction between the edge propagating at $V_{\textrm{g}}^{\star}$ and the maxima emanating from it and propagating at $V_{\varphi}^{\star}$ permits to understand previously unexplained observations. The propagation velocity extracted from t-VMC calculations in Ref. carleo2014, quantitatively agrees with the value $2V_{\varphi}^{\star}$ calculated here. Our analysis shows that it should thus be assimilated to the propagation of the local maxima, i.e.  the inner structure of the causal region NoDoubleInPrev . In contrast, the CE is determined by the raise of the envelop of these maxima, and moves at the velocity $2V_{\textrm{g}}^{\star}$. The same analysis applies to the results of t-DMRG calculations for a quench in the superfluid regime of the Fermi-Hubbard model manmana2009 .

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 cheneau2012 . 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}^{\star}>V_{\textrm{g}}^{\star}$, 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}^{\star}$) faster than the correlation cone (still at $V_{\textrm{g}}^{\star}$), and vanish when reaching it, see Figs. 1(b) and (d).

IV Long-range couplings

We now turn to long-range systems with power-law couplings, $J_{\mathbf{R},\mathbf{R}^{\prime}}\sim J/|\mathbf{R}-\mathbf{R}^{\prime}|^{\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_{\textrm{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. frerot2017, . 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. hauke2013, ; cevolani2015, ; cevolani2016, . 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 hastings2006 ; foss-feig2015 are all super-ballistic, we find a sub-ballistic CE, $t\sim R^{\beta_{\textrm{\tiny CE}}}$ with ${\beta_{\textrm{\tiny 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_{\textrm{\tiny m}}}$ with ${\beta_{\textrm{\tiny 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_{\textrm{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_{\textrm{sp}}=\left(2|c|zt/R\right)^{1/\left(1-z\right)}$. 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\geq 0$, into Eqs. (3) and (4), we find

[TABLE]

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{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

[TABLE]

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 ShortRangeLimit . Since $\chi=\gamma+D/2$ the CE is always sub-ballistic, ${\beta_{\textrm{\tiny CE}}}>1$. For the LRXY model and spin-spin correlations, we have $\nu=z=(\alpha-D)/2$, which yields ${\beta_{\textrm{\tiny CE}}}=1+\frac{D}{2\alpha}\left(2+D-\alpha\right)$. 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_{\textrm{\tiny CE}}^{\textrm{fit}}}}$, see details in App. D. For $\alpha=2.3$ and $D=1$ [Fig. 2(a)], we find ${\beta_{\textrm{\tiny CE}}^{\textrm{fit}}}\simeq 1.083\pm 0.013$, in good agreement with the theoretical value ${\beta_{\textrm{\tiny CE}}}\simeq 1.15$.

For the LRTI model, we have $z=\alpha-D$ and $\nu=0$. It yields the exponent ${\beta_{\textrm{\tiny 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_{\textrm{\tiny CE}}^{\textrm{fit}}}=1.28\pm 0.02$, in excellent agreement with the theoretical value ${\beta_{\textrm{\tiny CE}}}=1.3$. Note that the general formula for ${\beta_{\textrm{\tiny CE}}}$ matches the exact result of Ref. cevolani2015, 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. cevolani2016, 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

[TABLE]

The correlation maxima are thus always super-ballistic, ${\beta_{\textrm{\tiny m}}}<1$. For the LRXY model and $\alpha=2.3$ [Fig. 2(a)], we find the theoretical value ${\beta_{\textrm{\tiny 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_{\textrm{m}}=aR^{{\beta_{\textrm{\tiny m}}^{\textrm{fit}}}}+b$. For the parameters of Fig. 2(a), it yields ${\beta_{\textrm{\tiny m}}^{\textrm{fit}}}\simeq 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. richerme2014, and in rough agreement with the analysis of numerical calculations performed within the truncated Wigner approximation for 1D and 2D LRXY models schachenmayer2015b . The same result as Eq. (7) was found in Ref. frerot2018, , 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_{\textrm{\tiny 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_{\textrm{\tiny m}}^{\textrm{fit}}}\simeq 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 jurcevic2014 .

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 enveilling the physical information encoded in correlation spreading and how this can be extracted in the next generation of experiments (see also kormos2017 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 cheneau2012 , 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.

Acknowledgements.

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,

[TABLE]

is constructed using the particle operators $\hat{K}_{1}(\mathbf{R})\equiv\hat{n}_{\mathbf{R}}\left(\hat{n}_{\mathbf{R}}-1\right)$ and $\hat{K}_{2}(\mathbf{R},\mathbf{R}^{\prime})\equiv-\hat{a}_{\mathbf{R}}^{\dagger}\hat{a}_{\mathbf{R}^{\prime}}-\hat{a}_{\mathbf{R}^{\prime}}^{\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}^{\prime})=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. georges2010, .

A.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. (8) up to quadratic order. The resulting quadratic form is then diagonalized using standard Bogoliubov transformation (see for instance Refs. natu2013, ; cevolani2015, ). It yields the gapless dispersion relation

[TABLE]

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

After the quench, the connected one-body correlation function, $G(R,t)=\langle\hat{a}^{\dagger}_{R}(t)\hat{a}_{0}(t)\rangle-\langle\hat{a}^{\dagger}_{R}(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

[TABLE]

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

A.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. altman2002, ; huber2007, ; barmettler2012, . For $n\in\mathbb{N}^{*}$ and $U\gg Jn$, it yields low-energy excitations made of doublon-holon pairs of energy

[TABLE]

The maximum of the group velocity is at the center of the band, $k^{\star}=\pi/2$, where the group and phase velocities are $V_{\textrm{g}}^{\star}=(2n+1)J$ and $V_{\varphi}^{\star}=U/\pi$, respectively. Hence, for $U>\pi(2n+1)J$, the phase velocity exceeds the group velocity, $V_{\varphi}^{\star}>V_{\textrm{g}}^{\star}$. Note that this regime is well inside the Mott regime where Eq. (11) 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

[TABLE]

for the quench from $U_{\textrm{i}}=\infty$ to $U_{\textrm{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}^{\prime})$ 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}^{\prime})\equiv\hat{S}_{\mathbf{R}}^{x}\cdot\hat{S}_{\mathbf{R}^{\prime}}^{x}+\hat{S}_{\mathbf{R}}^{y}\cdot\hat{S}_{\mathbf{R}^{\prime}}^{y}$, $J(\mathbf{R},\mathbf{R}^{\prime})=-J/2|\mathbf{R}-\mathbf{R}^{\prime}|^{\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

[TABLE]

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 holstein1940 ; auerbach1994 ,

[TABLE]

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. (13) yields a quadratic Bose Hamiltonian, which can be diagonalized using canonical Bogoliubov transformations, see for instance Ref. frerot2017, . For $\epsilon=0$ (LRXY model), it yields the dispersion relation for $D=1$

[TABLE]

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

[TABLE]

where $P_{\alpha}(0)$ and $P_{\alpha}^{\prime}$ 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_{\textrm{g}}$ diverges in the infrared limit $k\rightarrow 0$.

The connected spin-spin correlation function along the $z$ direction for a quench from $\epsilon_{\textrm{i}}\neq 0$ to $\epsilon_{\textrm{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

[TABLE]

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

The linearization of the Holstein-Primakoff transformation holds for $|1/2-S^{x}_{\mathbf{R}}|\ll 1$, see Ref. auerbach1994, . For the calculations corresponding to Fig. 2 (a), we find $\max\{|1/2-S^{x}_{\mathbf{R}}|\}\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. frerot2017, 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}^{\prime})\equiv\hat{S}_{\mathbf{R}}^{x}\cdot\hat{S}_{\mathbf{R}^{\prime}}^{x}$ with a uniform magnetic field $h(\mathbf{R})=-2h$ and the algebraically decaying exchange amplitude $J(\mathbf{R},\mathbf{R}^{\prime})=2J/|\mathbf{R}-\mathbf{R}^{\prime}|^{\alpha}$, which yields

[TABLE]

The LRTI has two phases koffel2012 . In the $z$-polarized phase, the disprsion relation can be found using again the Holstein-Primakoff transformation

[TABLE]

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

[TABLE]

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

[TABLE]

where the gap $\Delta=2\sqrt{h\left[h+JP_{\alpha}(0)\right]}$ is finite, $c=\sqrt{\frac{h}{h+JP_{\alpha}(0)}}JP_{\alpha}^{\prime}$, and $z=\alpha-D$, see Ref. cevolani2016, . 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

[TABLE]

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

For the calculations corresponding to the Ising model in Fig. 2 (b), we find $\max\{|1/2-S^{z}_{\mathbf{R}}|\}=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.

D.1 LRXY model

For the LRXY model, we consider the time evolution of the connected spin-spin correlation function $G\left(R,t\right)$ 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^{\star}(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_{\textrm{\tiny CE}}}\simeq 1.15$ shown as a solid green line on the figure. We have also fitted a power law function, $t^{\star}\propto R^{{\beta_{\textrm{\tiny CE}}^{\textrm{fit}}}}$, to the blue points for $20<R<175$ (not shown on the figure). It yields ${\beta_{\textrm{\tiny CE}}^{\textrm{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 zhang2017 . Fitting the same algebraic function to the correlation edge in the range $10<R<30$ and $t<10/J$ yields ${\beta_{\textrm{\tiny CE}}^{\textrm{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_{\mathrm{m}}=aR^{{\beta_{\textrm{\tiny m}}^{\textrm{fit}}}}+b$ (dashed blue line). The fit yields ${\beta_{\textrm{\tiny m}}^{\textrm{fit}}}=0.634\pm 0.014$ in excellent agreement with the theoretical exponent ${\beta_{\textrm{\tiny m}}}=0.65$.

D.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_{\textrm{\tiny CE}}}}$ with ${\beta_{\textrm{\tiny CE}}}=1.3$ (solid green line with a fitted multiplicative factor). Moreover, fitting a power law function to the blue points for $4<R<30$ yields ${\beta_{\textrm{\tiny CE}}^{\textrm{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 sufficently 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 fro 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_{\textrm{m}}=aR^{{\beta_{\textrm{\tiny m}}^{\textrm{fit}}}}+b$ to the date (dashed blue line on the figure). The fit gives ${\beta_{\textrm{\tiny m}}^{\textrm{fit}}}=1.0045\pm 0.0003$ in excellent agreement with the theoretical exponent ${\beta_{\textrm{\tiny m}}}=1$.

We have also performed similar fits on the other local maxima, which confirms the asymptotic ballistric 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_{\textrm{\tiny m}}^{\textrm{fit}}}=1.0039\pm 0.0002$, which is also in excellent agreement with the theoretical prediction.