Lieb Robinson bounds and out of time order correlators in a long range spin chain

Lieb Robinson bounds quantify the maximal speed of information spreading in nonrelativistic quantum systems. We discuss the relation of Lieb Robinson bounds to out of time order correlators, which correspond to different norms of commutators $C(r,t) = [A_i(t),B_{i+r}]$ of local operators. Using an exact Krylov space time evolution technique, we calculate these two different norms of such commutators for the spin 1/2 Heisenberg chain with interactions decaying as a power law $1/r^\alpha$ with distance $r$. Our numerical analysis shows that both norms (operator norm and normalized Frobenius norm) exhibit the same asymptotic behavior, namely a linear growth in time at short times and a power law decay in space at long distance, leading asymptotically to power law light cones for $\alpha<1$ and to linear light cones for $\alpha>1$. The asymptotic form of the tails of $C(r,t)\propto t/r^\alpha$ is described by short time perturbation theory which is valid at short times and long distances.


I. INTRODUCTION
One of the most general concepts to study dynamical properties of quantum many-body systems is the dynamics of quantum information, generalizing the spreading of all possible types of correlations in the system. Of particular interest is the dynamical spreading of local operators 1-3 , which contains information about all correlation functions composed of these operators. While in relativistic systems, the spreading of information is limited by the speed of light, there is no such strict limit in nonrelativistic quantum mechanics. However, it was shown by Lieb and Robinson 4 that quantum systems with short range interactions exhibit a similar, non-universal speed limit, implying a causal structure. This emergent "light cone" has important consequences for the behavior of many-body interacting systems such as the area law of entanglement 5,6 , the decay of correlations 7 and stability of topological order 8 , but also for the timescales of thermalization [9][10][11] . Recently, a lot of progress has been made in establishing similar speed limits for the spreading of information in systems with long range interactions and it is clear that also in such systems, information can not spread infinitely fast [12][13][14][15][16][17][18][19][20][21][22][23] , however it is not always clear whether so far established analytical bounds on information spreading are tight for experimentally relevant lattice models. Power law decaying interactions are present in several quantum simulator platforms such as trapped ions 24,25 , Rydberg atoms 26 , ultra cold atoms 27 and superconducting qubits 28 and it is therefore important to obtain tight bounds on information scrambling and thermalization timescales.
A useful measure to quantify the spreading of an initially local operatorÂ i (t) is the commutator with another local operatorB i+r C(r, t) = ||[Â i (t),B i+r ]|| (1) whereB i+r serves as a probe and the operatorsÂ i and * dluitz@pks.mpg.de; lcolmena@pks.mpg.deB i+r act only on sites i and i + r of the system respectively.Â i (t) = exp (iĤt)Â i exp (−iĤt) is the operatorÂ i under time evolution in the Heisenberg picture and ||.|| denotes any matrix norm. Vanishing C(r, t) indicates that no information has traveled from site i to i + r at time t. It should be noted here that this generally depends very little on the choice of the operator in chaotic quantum many-body systems and only fine tuned situations exist, where a dependence on the operator at short times can be observed 12 .
In systems with short range interactions, C(r, t) is bounded within a "light-cone" region t > r/v where v is a velocity that depends on the microscopic model. This bound does not represent a strict cutoff, since exponential tails exist outside the light cone 4,29,30 .
Similarly, analytical bounds have been derived in systems with long range interactions decaying as a power law 1/r α with distance r. Hastings and Koma suggested a logarithmic bound t ∼ log(r) 13 for any α.
In the case of strongly long ranged systems with small α < D (D is the spatial dimension), the logarithmic bound is dominant 16,31 . Polynomial light cones have been proposed 14,17,32 of the form t ∼ r (α−2D)/(α−D+1) in the regime α > 2D which consistently recovers the linear light cone in the short range limit α → ∞. This bound was tightened to t ∼ r (α−2D)/(α−D) 18 . In general α large but finite is consistently found to exhibit asymptotically short-ranged behavior 12,19,21,33 . It was argued by Gong et al. that a linear light cone structure persists for α > D 15 , which is also supported by numerical simulations 12 . A stochastic model of operator spreading in long range interacting systems points to linear light cones for α ≥ D + 1 2 34 . For general quantum state transfer protocols, only a weaker bound for a linear light cone for α > 2D + 1 is valid 35,36 .
The analysis of analytical bounds 4,[13][14][15][16]18,31,32,35,37 is mostly concerned with the operator norm ||Ĉ(r, t)|| 2 (the largest singular value of the commutator matrixĈ(r, t) in Eq. (1)) because it encodes the "worst case" scenario, namely the fastest spreading modes in the system. On the other hand, numerical simulations of C(r, t) usually employ the square of the normalized Frobenius norm arXiv:2005.10257v1 [cond-mat.str-el] 20 May 2020 Ĉ (r, t) 2 F 11,12,33 , which is the average over the square of its singular values and is directly related to the out of time order correlator (OTOC) [38][39][40][41] as shown in Eq. (7). In Ref. 36 a bound on the Frobenius norm was established and found to be different from the bound on the operator norm. The Frobenius norm, associated with typical states, was shown to exhibit linear light cones for α > 5/2, while for the operator norm a weaker bound of α > 3 was found.
In a previous numerical study of long range interacting spin chains 12 , the asymptotic shape of the light cone of the OTOC (normalized Frobenius norm) was considered. In the present work, we are interested instead in the behavior of the operator norm of the commutator to compare the spreading of the fastest mode to that of typical modes (Frobenius norm). Interestingly, our analysis suggests, that both the average and the largest singular value ofĈ(r, t) have the same asymptotic behavior: We find a linear growth in time of Ĉ (r, t) 2 at short times, and a power law decay with distance at long distances with the exponent α, which can be understood from perturbation theory.

II. MODEL AND METHOD
We study the isotropic one-dimensional Heisenberg XXX model with long range interactions: where S γ i = σ γ i /2 are spin 1/2 operators acting on site i, with γ = x, y, z (σ γ i are the corresponding Pauli matrices). The interaction exponent α controls the range of the interactions and we use J = 1 throughout this paper. We do not use a rescaling of the coupling constant with system size L to make the energy extensive for small α, since this essentially only rescales our units of time. The model (2) conserves the total magnetization S z tot = iŜ z i and we focus on the largest magnetization sector S z tot = 0 for even L and S z tot = 1 2 for odd L, to maximize the accessible system sizes. The limit α → 0 corresponds to all-to-all interactions and α → ∞ is the nearest neighbor limit, which are both integrable points of the model. There is also a special integrable point at α = 2, the so called Haldane-Shastry model 42,43 .
For concreteness, we consider the dynamical spreading of the localŜ z i (t) operator, probed by the commutatorŝ We note that due to the SU(2) symmetry of the model, all S x i , S y i , S z i operators spread in the same way.
A. Matrix norms of the commutator In order to quantify the growth of the commutator norm C(r, t), we use two different matrix norms. The (normalized) Frobenius norm is defined as where s i are the singular values of the commutatorĈ(r, t) (and consequently s 2 i the eigenvalues ofĈ(r, t) †Ĉ (r, t)). The operator norm Ĉ (r, t) 2 , or 2-norm is defined by the largest singular value Therefore, the normalized Frobenius norm is always smaller than (or equal to) the operator norm Where for simplicity we have denoted C 2 (r, t) = Ĉ (r, t) 2 and C F (r, t) = Ĉ (r, t) F .

B. Out of time order correlator and relation to Frobenius norm of the commutator
Expanding the definition of the normalized Frobenius norm (4) for the commutatorĈ(r, t) = [Ŝ z i (t),Ŝ z i+r ] yields The correlation function 1 N Tr Ŝ z i (t)Ŝ z i+rŜ z i (t)Ŝ z i+r is known as the out of time order correlator (OTOC) and can be viewed as an infinite temperature four point function, where the partition function is given by the dimension of the Hilbert space Z = N .
In order to study the long distance behavior of this quantity, it is crucial to access large enough system sizes to ensure the convergence of our results in the thermodynamic limit and we therefore use dynamical typicality 44 for computing the trace which appears in the Frobenius norm Tr(C(r, t) † C(r, t)). This method consists of replacing the trace operation by the expectation value ψ|Ĉ(r, t) †Ĉ (r, t)|ψ where ψ is a random vector in the Hilbert space drawn from the Haar measure 45 , and averaging over random vectors |ψ . Eq. (7) is then boiled down to Tr Ĉ (r, t) †Ĉ (r, t) = ψ |ψ , where |ψ = C(r, t)|ψ , up to an error exponentially small in the system size L, requiring a very small number of random vectors (typically 1 . . . 100) for large enough systems. The operation C(r, t)|ψ is performed as a sequence of matrix-vector multiplications and several time propagations e −iHt |ψ of (intermediate) wave functions |ψ . These propagations can be performed efficiently using massively parallel sparse matrix Krylov space techniques. Technical details of this method for the calculation of the OTOC are discussed in Refs. 11, 12, and 46.

C. Operator norm of the commutator
In the present paper, our main focus is on the operator norm (2-norm) of the commutator which corresponds to the largest eigenvalue (equivalent to the largest singular value of C(r, t)) of the Hermitian form of the commutator iC( We use a matrix free implementation of the matrix vector product |ψ ← iC(r, t)|ψ , such that we never have to deal with dense matrices and use the Lanczos algorithm to obtain the largest eigenvalue of iC(r, t).
This means that we calculate Here, we have used the replacements |ψ 2 =Ŝ z i+r |ψ , . The matrix-free matrix-vector product involves again forward |ψ(t) = e −iĤt |ψ and backward |ψ(−t) = e iĤt |ψ real time evolution of the wavefunction, very similarly to the case of the OTOC 11,12 , for which we employ a Krylov space technique for the matrix exponential [47][48][49] . Matrix vector multiplications ofŜ z i operators with wave functions are trivial, since these operators are diagonal in the computational S z basis, and the entire algorithm thus requires only storage of a few vectors. This method gives access to the largest eigenvalue of the commutator with controlled accuracy up to system size L = 22 (N = 705432 in the zero magnetization sector). We note that for larger α and short distances r, the convergence of the Lanczos algorithm is particularly challenging due to small gaps in the spectrum. Lastly, for treating small systems L < 18, the calculations were performed using full exact diagonalization.
Throughout this paper, we fix the position of the spreading operator to i = 3 (the left most is indexed i = 0) in such a way that distances r = 0, 1, . . . , L − 4 are accessible (using open boundaries) and the reflection of the left information front does not interfere with propagation of the right one (which is the one we study in detail).

III. RESULTS
In the following, we analyze in detail the space-time profile of the operator norm of the commutator Ĉ (r, t) 2 of the long range XXX chain (2) and compare it to the case of the normalized Frobenius norm (OTOC), for which a very detailed analysis can be found in Ref. 12. We provide additional complementary data for the XYZ chain in Appendix B.

A. Causal space time region
We start our analysis by providing a qualitative comparison of the two norms of the commutatorsĈ(r, t) = [Ŝ z 3 (t),Ŝ z 3+r ] for different range of the interaction α and all distances r as a function of time. The synopsis of these results is shown in Fig. 1, where the top row shows the operator norm and the bottom row the normalized Frobenius norm (OTOC), while columns correspond to different ranges of the interaction α = 0.3, 0.7, 1.2, 2.4. Both norms are shown on the same color scale. Full lines show contour lines of the space time profile for various thresholds θ, extracted from the solution of the equation C(r, t) = θ to obtain the light "cone" t θ (r). It is clear already from a visual inspection of the two norms that the essential behavior is identical. Both norms reveal a clear causal space time region outside of which the commutator is very small, which means that almost no quantum information is communicated at short times and long distances for all α.
The contours are calculated for the same set of threshold values θ (indicated as vertical lines in the colorbar for clarity), clearly showing that the operator norm reaches a fixed threshold earlier than the Frobenius norm due to the property given by Eq. (6). Since no signal can travel faster than governed by the operator norm, it strictly limits the amount of quantum communication outside the causal region. The comparison between the operator and the Frobenius norm shows that typical modes in the system [singular values ofĈ(r, t)] travel significantly slower than the fastest mode (maximal singular value), an effect which is particularly pronounced at small α as can be seen from a comparison of the contour lines between the two norms.
The overall shape of the contour lines appears to be identical (with different prefactors). For large α, both norms are consistent with asymptotically linear light cones.
For intermediate α = 1.2, and large thresholds (black contour lines), we observe a "bump" in the case of the operator norm, which is likely nonuniversal and stems from the reflection at the left edge of the system, therefore we focus on smaller thresholds in these cases, where reflection does not (yet) interfere due to the observed causality.

B. Early time growth
In Fig. 2 we analyze the growth of the operator norm C 2 (r, t) = Ĉ (r, t) 2 as a function of time t for fixed distances r. The results are shown for a system of size L = 22 and reveal very clearly that the operator norm grows linearly in time (linear growth shown by dashed black lines for comparison). We have checked that this short time behavior is converged in system size. This is identical to the behavior of the normalized Frobenius norm (OTOC) 12 and is expected in the short time perturbative regime when t 1. We show how this linear growth arises from short time perturbation theory in Sec. III C, where it becomes also clear that nearest neighbor interactions lead to a different (power law) short time behavior.

C. Perturbation theory in the short time limit
At short times, we can use the Baker-Campbell-Hausdorff (BCH) formula The commutator Eq. (3) can then be written as For systems with long range interactions, the commutator to linear order [[Ĥ,Ŝ z i ],Ŝ z i+r ] is nonzero for any distance r, and is therefore the leading order at short times, leading to a dominant term linear in t. For the Therefore, the operator norm C 2 (r, t) to leading order in t reads For any finite α the operator norm Eq. (8) grows linearly in time and scales as r −α at long distance and short times. We note that this perturbative behavior is true for any choice of the norm. In Fig. 3 i ) which has support only on sites i − 1, i, i + 1 and therefore [[Ĥ,Ŝ z i ] 1 ,Ŝ z i+r ] vanishes as long as r > 1. Higher order terms in the BCH formula for r > 1 become only nonzero if a string of nontrivial Pauli matrices of length r is generated between sites i and i + r, and thus the leading order in the BCH formula reads for nearest neighbor interactions whereÔ(1) is given by the operator norm of a sum of Pauli strings of length r. In other words, at short times and outside the light cone the operator norm grows as a power law in time with an exponent given by the distance between the two operators. This is a quite general result, valid for any pair of local operators that are separated by a distance r larger than the support of the most extended term in the Hamiltonian. In Fig. 6 the exact time evolution of Eq. (8) for the XXX short range model is compared to the perturbation theory result Eq. (15). The power law growth C 2 (r, t) ∝ t r is in excellent agreement with the exact calculation at short times. For large values of α, there is a significant speedup of the growth of the commutator norm and a clear departure from the linear growth of C 2 (r, t) at intermediate times (cf. Fig. 2 lower panels). On the other hand, when α is large enough the commutator C 2 (r, t) is expected to exhibit a similar behavior as a short range interacting system, which corresponds to the limit α → ∞ in Eq. (2). The Hamiltonian of the long range model contains the nearest neighbor part plus longer distance couplings, decaying as 1/r α , which are strongly suppressed for α 1. Therefore, a dominant effect of the nearest neighbor part is expected for large α 12,22 . In Fig. 3 we compare C 2 (r, t) for the long range model (full colored lines) to the nearest neighbor model (red dashed) at fixed distances r. At short times, the long range model shows the perturbative growth r −α t/2 for all α, and speeds up at intermediate times. The initial growth is significantly faster than in the case of nearest neighbor interactions. For nearest neighbor interactions, the operator norm C 2 (r, t) grows much faster due to the large power law ∝ t r . Therefore, at later times and for α > 1, the nearest neighbor part catches up and dominates the overall growth of the commutator and leads to an asymptotic linear light cone.
Focusing only on large α, the operator norm of the commutator exhibits two kinds of grows: linear at short times (see Fig. 2) and short-ranged like at intermediate times (see Fig. 3). The short-range time evolution of C 2 (r, t) is well characterized by (t/r) r 50 , while the longrange part is described by Eq (14) t/r α in the limit t 1. These two results can be combined into a single expression: At short times the linear term on the right hand side is always dominant, at intermediate time the second term become dominant and the dynamics is short-range like.
At long distances there is a clear transition from a linear light cone for α > 1 to a power law light cone at α < 1, which can be understood with the following reasoning. The asymptotic form of C 2 (r, t) in Eq. (16) grows monotonically and the two terms compete. The light cone is given by the set of times t θ (r) as a function of distance r, for which C 2 (r, t) reaches a threshold value θ, i.e. C 2 (r, t θ (r)) = θ. It is clear that t θ (r) ≤ t c (r) = 2θr α , since this is the time the first (linear in t) term needs to reach the threshold. If the second (power law in t) term reaches the threshold first, we get a linear light cone, otherwise we get a power law light cone. We can estimate the power law term at long distances and t ≤ t c by t r /r r ≤ t r c /r r = (2θr α ) r /r r . Therefore, at t c , this term diverges for α > 1 and r → ∞, and overwhelms the linear term, leading to a linear light cone t θ (r) ∝ r. For α < 1, and r → ∞, this term is irrelevant and we are left with a power law light cone t θ (r) ∝ r α . Analogously, the linear term is bounded by the time when the short range front reaches the threshold, i.e. t/2r α ≤ rθ 1/r /2r α . The bounding term vanishes in the limit r → ∞ and α > 1 and diverges otherwise, which agrees with the bound to the short range term. This behavior is consistent with the numerical observation in Fig. 1.  4) computed from the discrete logarithmic derivative β = ∂ log r log C2(r, t) as a function of distance r for different system sizes L at fixed time t = 0.67. The asymptotic exponent β converges towards the interaction exponent α (dashed horizontal lines) for all values of α at long distances.
In Fig. 4 we analyze the behavior of C 2 (r, t) at long distances r outside the "light cone". It falls off as a power law at long distances, with an exponent that asymptotically approaches the interaction exponent α (see quantitative analysis using a discrete logarithmic derivative in Fig. 5). The same behavior was found previously 12 for the normalized Frobenius norm. This power law decay r −α is in perfect agreement with the prediction from perturbation theory in the short time limit given by Eq. (14), and seems to be valid asymptotically outside the causal region.
This analysis confirms the validity of short time perturbation theory (which is valid for any matrix norm) and shows that the asymptotic shape of the tails (outside the causal region) of long range interacting spin chains is given by C 2 (r, t) ∝ t/r α .

IV. CONCLUSIONS
The operator norm C 2 (r, t) of the commutator Eq. 1 has been examined in the XXX chain with long range interactions falling off as a power law r −α with distance r and interaction exponent α. In order to reach large enough systems to check the convergence of our results with system size, we introduce a Krylov space method for the direct calculation of the operator norm of the commutator C(r, t). We find a linear growth in t at early time t and a long distance decay outside the causal region given by r −α . Both the normalized Frobenius norm (directly related to OTOCs) and operator norm which correspond to the average and fastest information spreading modes in the system have the same asymptotic behavior (with different prefactors) of t/r α at long distance and short time, which is strikingly different from systems with nearest neighbor interactions, which instead exhibit a leading growth as a power law in time ∝ t r at short times and long distances.
For α > 1, the information front is dominated by the contribution from the nearest neighbor part of the Hamiltonian, which overtakes the initial linear growth of C 2 (r, t), inducing a linear light cone. These results are confirmed using a slightly different XYZ model in appendix B.
We conclude that the findings from the study of out of time order correlators in Ref. 12 provide information about Lieb Robinson bounds and agree with the behavior of the fastest spreading information mode in the system.

V. ACKNOWLEGMENTS
We acknowledge financial support from the Deutsche Forschungsgemeinschaft (DFG) through SFB 1143 (project-id 247310070). In the main text, section III C, the Baker-Campbell-Hausdorff (BCH) formula was employed for treating the time evolution of Eq. (1) at short times. From this analysis, the leading order of the operator norm C 2 (r, t) ≈ t r /r! was obtained. In Fig. 6 the exact time evolution C 2 (r, t) is shown along with the leading order in the BCH formula (dashed lines), with excellent agreement at short times. In order to test the universality of the results presented in the main text, we have also performed similar calculations for the long range XYZ Heisenberg model: with parameters J x = 0.9, J y = 1.2, J z = 0.7. The XXX model study in the main text is recovered by setting J x = J y = J z = 1. The XYZ model does not have U (1) symmetry, therefore the Hamiltonian is not block diagonal and we must deal with the full Hilbert space dimension N = 2 L . Analogously to the main text, we study the operator norm of the commutator The difference compared to Eq. (8) lies in the static operator that is nowŜ x i+r . The time evolution and operator norm computation are carried out using exact diagonalization.
We analyze the behavior of Eq. (B2) as function of both r and t. In Fig. 7 the causal regions of C 2 (r, t) are shown. For all values of α the overall shape of the causal region is the same as for the XXX long range model (see Fig. 1). Small α exhibit fast spreading with power law  causal regions and large α approach the linear light cone limit α → ∞. The long distance decay is also similar to what has been found in the main text (see Fig. 8 and 4) outside the light "cone" there is a power law decay of C 2 (r, t). Time evolution at fixed distance is displayed in Fig. 9 and is compatible with linear growth at short times t < 1 which was found also in the XXX version (see Fig. 2). In conclusion, the main features of both operator norm discussed in the main text are the same when considering a different commutator [Ŝ z i (t),Ŝ x i+r ] and a different long range model, namely the long range XYZ model. As expected only the interaction exponent α is crucial for characterizing the dynamics of the commutator Eq. (1) Applying perturbation theory Eq. (12) up to first order we get the following expression for the commutator when r > 0: yielding the following asymptotic form for the operator norm: In Fig. 9 this asymptotic form is compared with the exact time evolution yielding very good agreement. In conclusion, the BCH formula also predicts the short time behavior for the XYZ model.