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

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

Another dramatic effect of long-range interactions is the breakdown of the notion of causality in non-relativistic quantum models.For a wide class of lattice models with shortrange interactions, the Lieb-Robinson bound implies that correlations decay exponentially beyond a limit set by the system's maximum group velocity [37][38][39].This effective lightcone effect has been demonstrated for various models in both experiments [40][41][42] and numerics [43][44][45][46][47][48][49].In contrast, for sufficiently long-range interactions, causality breaks down and information can propagate arbitrarily fast [50,51].This is consistent with the absence of known Lieb-Robinson bounds [38,52,53] and the vanishing of the characteristic dynamical time scale in the thermodynamic limit [54].Conversely, algebraic interactions that decay fast enough are effectively short range, and one recovers ballistic propagation of information [55][56][57].The intermediate regime is, however, strongly debated, and whether a form of causality emerges remains an open question.Known bounds in principle allow for super-ballistic propagation [38,52,53] but they are challenged by numerical simulations, which point towards a significantly slower propagation [9,50,[58][59][60][61][62].The latter, however, re-ported different propagation scaling laws.Microscopic meanfield theory suggests that these apparent contradictions may be attributed to the coexistence of several signals governed by different dynamical scaling laws [63].This prediction, however, relies on a generic but non-universal form of the correlation functions and ignores beyond-mean-field effects.Experiments performed with trapped ions have reported bounded propagation [64,65] but they are limited to very small systems, which prevents extraction of the scaling laws and closure of the debate.
The aim of this work is to characterize the spreading of quantum correlations in the intermediate regime of a longrange spin system numerically.Specifically, we determine the scaling laws for the propagation of a variety of observables in the long-range transverse Ising (LRTI) chain using matrix product state simulations.We find that a weak form of causality emerges, characterized by generic algebraic scaling forms  ∼   , where the specific value of the exponent  depends on the observable.The spin-correlation and local-spin causal edges are both sub-ballistic, with the same dynamical exponent ( CE =  SE > 1).In the vicinity of the edge, however, the local maxima propagate differently, i.e. ballistically ( m = 1) for spin-spin correlations and super-ballistically ( m < 1) for local spins.In contrast, the Rényi entanglement entropies always propagate ballistically, irrespective of the range of interactions, in both the local and quasi-local regime.The analytic quasi-particle picture, based on linear spin wave theory (LSWT), accurately reproduces the numerics and provides a clear interpretation of the numerical results.The different algebraic space-time patterns of correlation functions provide an unambiguous fingerprint of the dynamical regimes of the model, suggesting the emergence of a dynamical phase-diagram.These correlation patterns can be directly measured in state-of-the-art experiments.Model and approach.-Thedynamics of the LRTI chain is governed by the Hamiltonian where Ŝ   (  = , , ) are the spin-1/2 operators on lattice site  ∈ [0,  − 1],  is the system size,  > 0 is the coupling energy, and ℎ is the transverse field.It can be realized on various quantum simulation platforms, including cold Rydberg gases [24,66,67] and artificial ion crystals, where the exponent  can be controlled via light-mediated interactions [16,18,19,64,65].At equilibrium, the phase diagram of the LRTI chain comprises two gapped phases separated by a second order quantum phase transition, see Fig. 1 and Ref. [7].For low fields ℎ and rather short-range couplings (high values of ), the nearest-neighbor anti-ferromagnetic couplings dominate and the system forms a staggered Néel-ordered phase along the  direction.For a large field ℎ and long-range couplings (low values of ), the spin-field interaction is favoured and a -polarized phase is formed.Out of equilibrium, the LSWT predicts three dynamical regimes (shaded by different colors in Fig. 1) [50,58]: For  ≥ 2 (local regime), the spin wave excitations are regular with bounded energies   and group velocities  g () =     ( is the momentum).This regime is reminiscent of the short-range case where correlations spread at finite speed, giving rise to a linear causality cone [37,38].For  < 1 (instantaneous regime),   features an algebraic, infrared divergence.There is no characteristic time scale and correlations spread arbitrarily fast.Finally, for 1 ≤  < 2 (quasi-local regime),   is bounded but  g () diverges as  g () ∼ 1/ 2− .Whether some form of causality emerges in the quasi-local regime remains debated [50,54,58,60,64,65,68], and, in the following, we mainly focus on this case.
Except where otherwise indicated, the results discussed below are obtained using time-dependent variational principle (TDVP) simulations within a matrix-product state (MPS) framework [7,50,69].Convergence of the calculations with the bond dimension has been systematically checked.Our results are summarized on the dynamical phase diagram of Fig. 1.The main features of the observable that we consider are all described by an algebraic space-time dependence of their edges and maxima, shown in the insets.
Spin correlations and global quench.-Wefirst consider the spreading of the connected spin correlations,   (, ) = Ŝ  0 () , along the directions  = , . Figure 2(a) shows a typical TDVP result for   (, ), for a global quench in the quasi-local regime,  = 1.7, from (ℎ/) i = 50 to (ℎ/) f = 1, both in the -polarized phase.The initial state of the system is the ground state of the Hamiltonian with (ℎ/) i .Similar results are found when scanning the values of (ℎ/) i , (ℎ/) f , and 1 <  < 2. The correlation pattern shows a series of maxima which propagate algebraically [notice the log-log scale in Fig. 2(a)].In the asymptotic long-time and long-distance limits, they are well fitted by the scaling law  ∼   m (dashed blue lines).The correlation edge (CE), which sets the causality horizon, cannot, however, be inferred from the behavior of the maxima [47].To find the CE, we track the points of the  − plane where the correlations reach a fraction  of their maximum.Scanning  ∈ [0.01, 0.1], we find that the CE is well fitted by the algebraic scaling law  ∼   CE (solid line), with  CE nearly independent of  [70].Figure 2(b) shows the results of the fits (empty symbols) versus the exponent .
The numerical results for   (, ) have an extra checkerboard structure inside the causal region due to the antiferromagnetic coupling in Eq. (1).Once this structure is taken into account, we can identify a CE and local maxima, and extract the corresponding exponents.The results are similar to those found for   (, ), see filled symbols in Fig. 2(b) [70].
Comparing the fitted dynamical exponents  m and  CE to the predictions of the LSWT [63], we find an excellent agreement, as shown in Fig. 2(b).While the maxima spread ballistically,  m 1, the CE is sub-ballistic with  CE  LSWT = 3 −  > 1.This is characteristic of gapped long-range models, such as the LRTI model in the z-polarized phase [63].It confirms the emergence of a weak, slower-than-ballistic form of causality for 1 <  < 2.
In the local regime,  ≥ 2, we find that both the maxima and the CE spread ballistically,  m  CE 1.The spreading velocities are, however, different from each other, which is characteristic of a non-phononic excitation spectrum.The CE velocity is  CE 2 g ( * ), i.e. twice the maximum group velocity of the spin waves, while that of the maxima is  m 2  ( * ), i.e. twice the phase velocity   () =   / at the quasi-momentum  * where the group velocity is maximum, see Fig. 2(c).
We have also performed calculations for large quenches.While the LSWT is well justified only for weak quenches, we have found that the dynamical exponents it predicts are extremely robust, even for quenches across the critical line.We found dynamical scaling laws in very good agreement with those reported on Figs.2(b) and (c), although the signal is blurred as compared to weak quenches owing to the proliferation of quasi-particles when the quench amplitude increases [70].
Local magnetization and local quench.-Wenow consider the dynamics of another quantity, namely the local magnetization Ŝ  () , and perform a local quench.We initialize the system in the ground state of the -polarized phase and flip the central spin.Figure 3(a) shows a typical TDVP result for the quantity 1/2 − Ŝ  () versus the time  and the distance  from the flipped spin, in the quasi-local regime.The result displays, as in the case of global quenches, a twofold algebraic structure.Fitting the maxima (dashed blue line) and the spin edge (SE, solid green line) as previously, we extract the dynamical exponents  m and  SE plotted in Fig. 3(b) (blue disks and green diamonds, respectively).The SE follows the same sub-ballistic scaling law as the CE in the previous case,  SE 3 −  (solid green line).In contrast, the maxima of the local spin spreads faster than those of the spin correlations, and here we find  m < 1, corresponding to a super-ballistic propagation.
To understand this behavior, let us first recall that the standard LSWT generically predicts a super-ballistic propagation of the maxima in gapless systems [63]  |↑ ↑ ... ↑ , and thus lives in the first excited manifold.The ground state and the gap are thus irrelevant, and we may expect super-ballistic propagation of the maxima consistent with the TDVP results.This argument applies to any observable for such a local quench.More specifically, the LSWT applied to Ŝ  () for the local quench yields  m =  − 1, in very good agreement with the TDVP results, see Fig. 3(b) [70].
In the local regime ( > 2) we find that Ŝ R () propagates ballistically.Both this property and the SE velocity extracted from the TDVP calculations are in good agreement with the LSWT analysis, see Fig. 3(c).Note that in this regime, we do not observe maxima propagating at a different velocity.This is also consistent with the LSWT analysis.In the local regime, the quantity Ŝ  () is the sum of several contributions, each with a twofold structure but that are in phase quadrature and cancel each other.A similar effect has been found for density correlations deep in the Mott insulator phase of the Bose-Hubbard model [47,70].
Entanglement entropy.-Wefinally study the spreading of quantum information after the same local quench.It may be measured via the Rényi entropies of the reduced density matrix of a block, Figure 4(a) shows a typical TDVP result for the von Neumann entropy,  → 1 [71], and the same local quench as in Fig. 3(a).Similar results are found for the other Rényi entropies [70].We find that S  (, ) is a monotonic function of both position and time, and no local maxima such as those found previously for the correlation function and the magnetization are observed.This is consistent with a causal spreading of information and the expectation that the entanglement entropy decreases with the size of the smaller partition ().The entanglement edge (EE) is clearly visible in Fig. 4(a) and a fit to the algebraic scaling law  ∼    EE allows us to extract the exponent.Varying the threshold in a wide range,  ∈ [0.2, 0.8], we find   EE ≈ 1 within error bars, irrespectively of the interaction range, see cyan points in Fig. 4(b).The error bars correspond to the variation of   EE with , which is due to finite-size effects (see below).Our results imply that the propagation of entanglement is close to ballistic in both the local and quasilocal regimes.This contrasts with the behavior of the oneand two-point observables considered thus far, which exhibit sub-ballistic causal edges.It is a direct consequence of the fact that the bipartite entanglement entropy is a highly non-local quantity, which takes into account all entangled pairs on either side of the bipartition  [39,[72][73][74].
More precisely, we have analytically computed the entanglement entropy within LSWT [70].The results are in good agreement with those extracted via TDVP for the same system sizes, see Fig. 4(b).The ballistic spreading of entanglement is further confirmed by the LSWT results obtained in much larger systems, with significantly smaller error bars, see Fig. 4(c).Bipartitioning the quenched initial state within LSWT yields an entanglement Schmidt rank bounded by 2. This is consistent with our TDVP results at all times and the saturation of the entanglement entropies to S  (,  → ∞) log(2) 0.69, see Fig. 4(c).The two eigenvalues of the reduced density matrix ρ ,  1 (, ) and  2 (, ) = 1− 1 (, ), can then be computed analytically.In the asymptotic limit and for not too small values of /, we find  2 (, ) ∝   [75,76].Hence, the -order Rényi entropy is a function of the ratio /, This confirms the ballistic propagation of the entanglement entropy (  EE = 1) consistent with the results of Fig. 4 for  = 1 and other Rényi orders  [70].
Conclusion.-Our results show the emergence of a weak form of causality in the intermediate regime of the long-range Ising model, characterized by algebraic propagation laws with exponents that depend on the observables and the range of interactions.While local spins and spin correlations both have a sub-ballistic propagation edge,  ∝   , with  > 1, the causal region is characterized by local maxima propagating super-ballistically and ballistically, respectively.The distinction between the causal edge and the local maxima, which can show drastically different dynamical behaviors, is thus pivotal in the characterization of causality in long-range quantum systems.In contrast, the propagation of entanglement is ballistic in both the local and quasi-local regimes, and the causal region is featureless.
These results call for future experimental and theoretical work.On the one hand, our predictions are directly relevant to quantum simulators using for instance trapped ions, where the interaction range can be controlled.While analysis of spin and correlation spreading in first experiments have been limited by finite-size effects [64,65], systems with more than 50 ions, comparable to the system size used in our simulations, are now accessible [77] and Rényi entanglement entropies can now be measured in trapped ion platforms [78,79] and digital quantum computers [80].On the other hand, it would be interesting to further test the robustness of the observed algebraic scaling laws by quantitatively investigating the dependence (if any) of the exponents on the strength of the quenches, the phases of the model and their values in different models, such as the long-range XY [9,63], Heisenberg [81,82], and Hubbard [44,58] models, as well as in dimensions higher than one.This extended analysis could allow the identification of dynamical universality classes, i.e. models which share the same algebraic laws for correlations out-of-equilibrium.

Linear spin wave theory
In the -polarized phase, the LRTI model can be diagonalized using a Holstein-Primakoff transformation [84], where â () is the annihilation operator of a boson at site  and time .Inserting these formulas into the LRTI Hamiltonian [Eq.( 1) of the main paper], we find a quadratic Bose Hamiltonian, see for instance Refs.[50,54] for details.The latter is readily diagonalized using a standard Bogoliubov transformation [85].The elementary are magnons (spin-wave excitations), characterized by the gapped spectrum where   () =     /  denotes the Fourier transform of the long-range term in the Hamiltonian.For a quench on the exchange amplitude  from  i to  f , the   (, ) correlation function can be cast in the form [63]   (, ) () where the index "i" refers to the pre-quench (initial) Hamiltonian and the index "f" to the post-quench (final) Hamiltonian.The function F () gives the weight of each quasi-particle according to their mode .It depends on the observable and the quench.
To characterize the asymptotic behavior, ,  → +∞ along a constant line /, one can then rely on the stationary-phase approximation, applied to Eq. (S4).It yields where  is a constant phase irrelevant to our study,  sp is the stationary-phase quasi-momentum, given by the solution of the equation 2 g ( sp ) = 2    sp = / > 0 with  g the group velocity associated to the spin-wave excitations.
Quasi-local regime -The quasi-local regime corresponds to the case where the quasi-particle energy   is finite over the whole first Brillouin zone,  ∈ [−, +], but the group velocity  g () presents a divergence.The space-time behavior of   in the vicinity of the CE (correlation edge) is dominated by the quasi-particles propagating with the highest group velocities.Due to the infrared divergence, only the limit behavior in  → 0 is relevant.There the long-range term reads as   () ≈   (0) +   | | −1 , with   (0) > 0,   < 0, and the excitation energy is   Δ − | |  , where Δ = 2 √︁ ℎ[ℎ +   (0)] > 0 is the gap,  =  − 1 ≥ 0 is the dynamical exponent, and  = √︁ ℎ/[ℎ +   (0)] |  | > 0 is a prefactor.Then, computing  sp and injecting it into Eq.(S5), we find with where  ≥ 0 is the scaling exponent of the amplitude function in the infrared limit, F () ∼ | |  .It follows from Eq. (S4) and the approximation of   in the infrared limit that  = 0 for the   spin correlation function.On the one hand, the CE is found by imposing the condition that the prefactor is constant.The latter leads to the algebraic form  ∝   CE with  CE = / = 3 − .Since 1 ≤  < 2 in the quasi-local, the CE is always sub-ballistic, i.e.  CE > 1.On the other hand, the spreading law of the local extrema is determined by the equation The maxima are thus ballistic, i.e.  ∼   m with  m = 1. Figure S4 shows the same calculations, now performed in the local regime,  = 3.Similar conclusions as for the quasi-local regime hold.The   (, ) correlation function is rapidly blurred.In contrast, the   (, ) is more robust and allows us to fit the dynamical exponents  m and  CE , yielding values nearly independent of (ℎ/) i .In fact, for the   (, ) correlation function, we find no significant impact of increasing the strength of the quench.where  (, ) = ∞ =0 ( + ) − is the Hurwitz zeta function [75,76] and  = −2 −1 = −3 −2 .Using the series representation of the Hurwitz zeta function in the limit of large  and positive , this may be written as where   are the Bernoulli numbers, and Γ() = ´∞ 0 e −  −1 d is the gamma function [76].In the quasi-local regime,  > 2 and so in the limit of large  we may neglect the second and third term which contain contributions of the order equal to  − and higher, and only consider the leading term of order  1− .It yields which scales in  and  with the same power.Note that the stationary phase approximation is not valid when taking the limit  → ∞ alone, and the above equation does not hold in this limit.Consequently, any analytic function of  2 , particularly Rényi entropies, will scale in  and  with the same power, so that S  (, ) This yields a dynamical exponent of the EE of unity,     = 1 ∀ , in close agreement with our numerical results using both TDVP and LSWT.The additional finite-size corrections neglected in moving from Eq. (S9) to Eq. (S10) are likely responsible for the deviation from ballistic behavior seen in our simulations on small system sizes.
Local regime -The calculation in the local regime  > 2 proceeds similarly.However, both the quasi-particle energy   and group velocity  g () are finite over the whole Brillouin zone such that there exists a quasi-momentum  ★ where the group velocity is maximum, max    () =   ( ★ ).The stationary phase condition in the local regime thus has only a solution for / ≤   ( ★ ).The eigenvalues of the reduced density matrix are then dominated by the momenta  ∼  sp =  ★ , which is independent of  and , and thus scale ballistically.Finally, any analytic function of the eigenvalues will also scale ballistically.

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

Figure 2 .
Figure 2. Spreading of spin correlations in a global quench, with system size  = 48.(a) TDVP results for   and a quench from (ℎ/) i = 50 to (ℎ/) f = 1 in the quasi-local regime at  = 1.7.The solid green and dashed blue lines are fits to the CE and the extrema, respectively.(b) Dynamical exponents  CE (green diamonds) and  m (blue disks), fitted form results as in panel (a) for   (empty symbol) and   (full symbol), and comparison to the LSWT predictions (solid green and dashed blue lines).(c) Spreading velocities  CE (green diamonds) and  m (blue disks), in the local regime and comparison to the LSWT predictions (solid green and dashed blue lines).

Figure 3 .
Figure 3. Spreading of the local magnetization for a local quench, with system size  = 48.(a) TDVP results for the quantity 1/2 − Ŝ  () versus the time and the distance from the flipped spin for ℎ/ = 50 and  = 1.8.(b) Dynamical exponents  SE (green diamonds) and  m (blue disks), fitted form results as in panel (a) and comparison to the LSWT predictions (solid green and dashed blue lines).(c) Spreading velocity  SE (green diamonds) and comparison to the LSWT prediction in the local regime (solid green line).

Figure 4 .
Figure 4. Spreading of the von Neumann entanglement entropy S =1 (, ) for the same quench as in Fig. 3(a).(a) TDVP results for ℎ/ = 50,  = 1.8, and a system size  = 96.The solid green line marks a power law fit to the EE with  = 0.5, yielding     = 0.899 ± 0.005.(b) Dynamical exponents for the EE obtain via TDVP (cyan downwards triangles) and LSWT (magenta upwards triangles), with error bars corresponding to the variations with respect to the threshold .(c) Dynamical exponents obtained via LSWT for a larger system,  = 512.

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

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