Interplay of soundcone and supersonic propagation in lattice models with power law interactions

Traditionally, the study of lattice models has focused on Hamiltonians where interactions and/or hopping is restricted to a few neighboring sites. Only recently there has been a surge of interest in long-range interacting systems where interaction strengths or hopping amplitudes decay like a power law ${r}^{-\alpha }$ at large distances $r$. This interest was triggered on the experimental side by progress in the control of ultra-cold atoms, molecules, and ions, which led to the realization of a variety of long-range systems. Examples include magnetic atoms [1], polar molecules [2], trapped ions [3–6], Rydberg atoms [7], and others. On the theoretical side, intriguing physical effects and properties have been predicted for long-range interacting quantum systems, including nonequivalent statistical ensembles and negative response functions [8, 9], equilibration time scales that diverge with system size [10–12], prethermalization [13, 14], and others.

In this article we study the propagation in time and space of various physical quantities, and this is another topic where long-range interactions lead to peculiar behavior. A number of papers devoted to this topic have appeared in the past two years, reporting results on the spreading of correlations, information, or entanglement [15–21]. In short-range systems, all these quantities are known to propagate approximately within a soundcone, reminiscent of the lightcone in relativistic theories, with only exponentially small effects outside the cone. This behavior is termed quasilocality and was rigorously proved by Lieb and Robinson for a class of short-range interacting lattice models [22]. In the presence of long-range interactions this picture is altered significantly: the concept of a group velocity breaks down, and the spreading of correlations, information, or entanglement may speed up dramatically. This, in turn, has a bearing on all kinds of dynamical properties, and one might hope to harness long-range interactions for fast information transmission, improved quantum state transfer, or other applications.

Much of our understanding of propagation in long-range systems comes from analytical or numerical studies of model systems, where for example correlations or entanglement between lattice sites i and j are calculated as functions of time t and spatial separation ${\text{d}}(i,j)$. Typical examples of such results, similar to some of those in [15–21], are shown in figure 1 for a number of different models, physical quantities, and exponents α. For larger α (figure 1 right), the behavior is reminiscent of the short-range case, with only small effects outside a cone-shaped region. For small α (figure 1 left), correlations propagate faster than any finite group velocity would permit, and are mostly confined to a region with power law-shaped boundaries. For intermediate α (figure 1 center), a crossover from cone-like to faster-than-cone behavior is observed. While these three regimes seem to be typical and occur in many of the models studied, notable exceptions (some of which will be discussed further below) do occur and lead to a more complicated overall picture.

Zoom In Zoom Out Reset image size

Besides model calculations, Lieb–Robinson-type bounds have contributed significantly to our understanding of propagation in long-range interacting models. The first result of this kind

$$\begin{eqnarray}&&\parallel \left[{O}_{A}(t),{O}_{B}(0)\right]\parallel \leqslant C\parallel {O}_{A}\parallel \parallel {O}_{B}\parallel {\displaystyle \frac{| A| | B| ({e}^{v| t| }-1)}{{[{\text{d}}(A,B)+1]}^{\alpha }}},\end{eqnarray} \tag{ 1 }$$

valid for exponents α larger than the lattice dimension D, was reported by Hastings and Koma [23]. Here, $A,B\subset \Lambda$ are non-overlapping regions of the lattice Λ, and ${O}_{A}(0)$ and ${O}_{B}(0)$ are observables supported only on the subspaces of the Hilbert space corresponding to A and B, respectively. $\parallel \cdot \parallel$ denotes the operator norm, and ${\text{d}}(A,B)$ is the graph-theoretic distance between A and B ⁴ . The relevance of the bound (1) lies in the fact that a number of physically interesting quantities, like equal-time correlation functions, can be related to the operator norm of the commutator on the left-hand side of (1), so that similar bounds hold also for these physical quantities [24, 25]. For any α, a contour plot of the bound equation (1) looks qualitatively like the plot in figure 1 (left), although with logarithmic contour lines instead of power laws. This implies that, while correct as a bound for all $\alpha \gt D$, the shape of the propagation front (figure 1 center and right) is not correctly reproduced by (1) for intermediate or large values of α. Another bound put forward in [26] improves the situation for the case of large α, but turns out to be weaker than (1) for smaller values ⁵ . Summarizing the situation, the existing Lieb–Robinson-type bounds struggle to reproduce the transition from cone-like to faster-than-cone propagation for intermediate α as in figure 1 (center) ⁶ . For small α, no bounds have been published so far.

In this article we prove general bounds, complemented by model calculations, that help to establish a comprehensive and consistent picture of the various kinds of propagation behavior that occur in long-range interacting lattice models. We extend Lieb–Robinson-type bounds to strong long-range interactions where $\alpha \lt D$. This is complemented by model calculations showing that, even in the regime $\alpha \lt D$ of strong long-range interactions, cone-like propagation may be a dominant feature. We also prove that faster-than-cone propagation can occur for all $\alpha \lt 2$ in any spatial dimension, and this answers a question put forward in [28]. For intermediate exponents α, we advocate the use of a Lieb–Robinson-type bound in the form of a matrix exponential, which is tight enough to capture the transition from a cone-like to a faster-than-cone propagation as in figure 1 (center), and is also computationally efficient.

For deriving analytical results in the regime $\alpha \lt D$, an understanding of the time scales of the dynamics turns out to be crucial. The presence of strong long-range interactions is known in many cases to cause a scaling of the relevant time scales with system size [10, 11, 14, 29–31]. For long-range quantum lattice models the fastest time scale $\mathcal{T}\propto {N}^{-q}$ was found to shrink like a power law with increasing system size N, where q is a positive exponent [30, 31]. This observation makes clear why previous attempts to derive a Lieb–Robinson-type bound for $\alpha \lt D$ failed: in the large-N limit the dynamics becomes increasingly faster, and hence propagation is not bounded by any finite quantity. Considering evolution in rescaled time $\tau ={{tN}}^{q}$ can resolve this problem and allows us to obtain a finite bound in the thermodynamic limit.

On an arbitrary D-dimensional lattice Λ with N sites we consider the Hilbert space

$$\begin{eqnarray}&&\mathcal{H}=\underset{i=1}{\overset{N}{{\displaystyle \otimes }}}{\mathcal{H}}_{i}\end{eqnarray} \tag{ 2 }$$

with finite-dimensional local Hilbert spaces ${\mathcal{H}}_{i}$. On $\mathcal{H}$ a generic Hamiltonian

$$\begin{eqnarray}&&H={\displaystyle \sum }_{X\subset \Lambda }{h}_{X}\end{eqnarray} \tag{ 3 }$$

with n-body interactions is defined, with local Hamiltonian terms h_X compactly supported on the finite subsets $X\subset \Lambda$. The Hamiltonian is required to satisfy the following two conditions.

• (i)  
  Boundedness

  $$\begin{eqnarray}&&{\displaystyle \sum }_{X\ni i,j}\parallel {h}_{X}\parallel \leqslant {\displaystyle \frac{\lambda }{{[1+{\text{d}}(i,j)]}^{\alpha }}}\end{eqnarray} \tag{ 4 }$$

  with a finite constant $\lambda \gt 0$. This condition, also used in [23], is a generalization of the definition of power law-decaying interactions, and it reduces to the usual definition in the case of pair interactions, i.e., when X consists only of the two elements i and j.
• (ii)  
  Reproducibility

  $$\begin{eqnarray}&&{\mathcal{N}}_{\Lambda }{\displaystyle \sum }_{k\in \Lambda }{\displaystyle \frac{1}{{[1+{\text{d}}(i,k)]}^{\alpha }{[1+{\text{d}}(k,j)]}^{\alpha }}}\leqslant {\displaystyle \frac{p}{{[1+{\text{d}}(i,j)]}^{\alpha }}}\end{eqnarray} \tag{ 5 }$$

  for finite $p\gt 0$, with

  $$\begin{eqnarray}&&{\mathcal{N}}_{\Lambda }=1/\underset{i\in \Lambda }{\mathrm{sup}}{\displaystyle \sum }_{j\in \Lambda \setminus \{i\}}{\displaystyle \frac{1}{{[1+{\text{d}}(i,j)]}^{\alpha }}}.\end{eqnarray} \tag{ 6 }$$

The lattice-dependent factor ${\mathcal{N}}_{\Lambda }$ is the same that is frequently used to make a long-range Hamiltonian extensive [10, 32], but we use it here for a different purpose. Asymptotically for large regular lattices, one finds [10]

$$\begin{eqnarray}{\mathcal{N}}_{\Lambda }\sim \left\{\begin{array}{ccc}{c}_{1}{N}^{\alpha /D-1} & \mathrm{for} & 0\leqslant \alpha \lt D,\\ {c}_{2}/\mathrm{ln}N & \mathrm{for} & \alpha =D,\\ {c}_{3} & \mathrm{for} & \alpha \gt D,\end{array}\right.\end{eqnarray} \tag{ 7 }$$

with α-dependent positive constants c₁, c₂, and c₃. Equation (5) is a modified version of one of the requirements for the proof in [23], but due to the modification by the factor ${\mathcal{N}}_{\Lambda }$ the condition is satisfied for a larger class of models, including regular D-dimensional lattices with power law-decaying interactions with arbitrary positive exponents α [33]. For the above described setting we derive in appendix A.2 the Lieb–Robinson-type bound

$$\begin{eqnarray}&&\parallel [{O}_{A}(\tau {\mathcal{N}}_{\Lambda }),{O}_{B}(0)]\parallel \leqslant C\parallel {O}_{A}\parallel \parallel {O}_{B}\parallel {\displaystyle \frac{| A| | B| ({e}^{v| \tau | }-1)}{p{[{\text{d}}(A,B)+1]}^{\alpha }}}\end{eqnarray} \tag{ 8 }$$

in rescaled time

$$\begin{eqnarray}&&\tau =t{/\mathcal{N}}_{\Lambda }.\end{eqnarray} \tag{ 9 }$$

This bound reproduces qualitative features of supersonic propagation (as in figure 1 left), and also accounts for the system-size dependence of the time scale of propagation for exponents $\alpha \lt D$. While the bound ensures well-defined dynamics in rescaled time τ in the thermodynamic limit, it describes a speed-up in physical time t of the propagation with increasing lattice size, as illustrated in figure 2.

Zoom In Zoom Out Reset image size

For long-range models with intermediate exponents, in the range $3\lt \alpha \lt 6$ or even a bit larger, one observes an interplay of cone-like and supersonic propagation (figure 1 center). This is the most relevant regime for experimental realizations of long-range interactions by means of cold atoms or molecules, but a theoretical description of the shape of the propagation front turns out to be challenging. Existing bounds [26] are discussed in appendix A.3. Here we report bounds that capture the features of the propagation front as observed in long-range models with intermediate exponents, showing a clear and sharp crossover from cone-like to supersonic propagation.

As in section 2, our setting is a D-dimensional lattice Λ consisting of N sites and a Hilbert space (2) with finite-dimensional local Hilbert spaces. We consider a generic Hamiltonian with pair interactions

$$\begin{eqnarray}&&H={\displaystyle \frac{1}{2}}{\displaystyle \sum }_{k\ne l}^{N}{h}_{{kl}},\end{eqnarray} \tag{ 10 }$$

where the pair interactions h_kl are bounded operators supported on lattice sites k and l only. As observables O_A and O_B we consider bounded operators that are supported on single sites $A=\{i\}$ and $B=\{j\}$. In this setting, we prove in appendix A.1 a bound in the form of an N × N matrix exponential,

$$\begin{eqnarray}&&\parallel [{O}_{i}(t),{O}_{j}(0)]\parallel \leqslant 2\parallel {O}_{i}\parallel \parallel {O}_{j}\parallel \left(\mathrm{exp}{[2\kappa J| t| ]}_{i,j}-{\delta }_{i,j}\right),\end{eqnarray} \tag{ 11 }$$

where J is the interaction matrix with elements

$$\begin{eqnarray}&&{J}_{k,l}=\parallel {h}_{{kl}}\parallel \end{eqnarray} \tag{ 12 }$$

and $\kappa ={\mathrm{sup}}_{n\in \Lambda }{{\displaystyle \sum }}_{k}{J}_{n,k}$. In one-dimensional homogeneous lattice models the interaction matrix J is of Toeplitz type and thus (11) can be evaluated in $\mathcal{O}({N}^{2})$ time using the Levinson algorithm [34]. For translationally invariant one-dimensional systems, J is a circulant matrix, which permits an analytical solution of (11) by means of Fourier transformation.

The bound (11) is tighter than the bounds in [23, 25, 26], and the crossover from cone-like to supersonic propagation is nicely captured (see figure 3). Due to its form as a matrix exponential, the bound is less explicit than others, and asymptotic properties are not easily read off. But since the calculation of a matrix exponential scales polynomially in the matrix dimension N (like $\mathcal{O}({N}^{3})$ or even faster [35]) the bound can easily be evaluated for large lattices up to $\mathcal{O}({10}^{4})$ on a desktop computer. This is orders of magnitude larger than the sizes that can be treated by exact diagonalization, and covers the system sizes that can be reached for example with state-of-the-art ion trap based quantum simulators of spin systems [36]. Different from other bounds of Lieb–Robinson-type, our matrix exponential bound is computed for the exact type of interaction matrix realized in a specific experimental setup. This improves the sharpness of the bound, and can make it a useful tool for investigating all kinds of propagation phenomena in lattice models of intermediate system size.

Zoom In Zoom Out Reset image size

The bounds discussed in sections 2 and 3 are valid for arbitrary initial states, and therefore it may well happen that propagation for a given model and some, or even most, initial states is significantly slower than what the bound suggests. Indeed, linear (cone-like) propagation was observed in model calculations even for moderately large exponents like α = 3 [15–19]. But, as we show in the following, such cone-like propagation can, for suitably chosen initial states, even persist into the strongly long-range regime $0\lt \alpha \lt D$. In this and the next section we analyze free fermions on a one-dimensional lattice with long-range hopping, which is arguably the simplest model to illustrate cone-like propagation in long-range models and to explain the observation on the basis of dispersion relations and density of states. While strictly speaking such a long-range hopping model does not meet the conditions under which Lieb–Robinson bounds have been proved, it proves helpful for understanding the conditions under which cone-like propagation may or may not be observed in other long-range interacting models.

4.1. Long-range hopping model

Consider a free fermionic hopping model in one-dimension with periodic boundary conditions

$$\begin{eqnarray}&&H=-{\displaystyle \frac{1}{2}}{\displaystyle \sum }_{j=1}^{N}{\displaystyle \sum }_{l=1}^{N-1}{{\text{d}}}_{l}^{-\alpha }\left({c}_{j}^{\dagger }{c}_{j+l}+{c}_{j+l}^{\dagger }{c}_{j}\right),\end{eqnarray} \tag{ 13 }$$

where ${c}_{j}^{\dagger }$, c_j are fermionic creation and annihilation operators at site j. We choose long-range hopping rates $\propto \;{{\text{d}}}_{l}^{-\alpha }$, where

$$\begin{eqnarray}{{\text{d}}}_{l}=\{\begin{array}{ccc}l & \mathrm{if} & l\leqslant \ N/2,\\ N-l & \mathrm{if} & l\gt N/2,\end{array}\end{eqnarray} \tag{ 14 }$$

is the shortest distance between two sites on a chain with periodic boundary conditions. A Fourier transformation brings the Hamiltonian into diagonal form

$$\begin{eqnarray}&&H={\displaystyle \sum }_{k}\epsilon (k){a}_{k}^{\dagger }{a}_{k}\end{eqnarray} \tag{ 15 }$$

with

$$\begin{eqnarray}&&{c}_{j}={\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \sum }_{k}{\rm{e}}^{\rm{i}{kj}}{a}_{k}.\end{eqnarray} \tag{ 16 }$$

and dispersion relation

$$\begin{eqnarray}&&\epsilon (k)=-{\displaystyle \sum }_{l=1}^{N-1}{\displaystyle \frac{\mathrm{cos}\left({kl}\right)}{{d}_{l}^{\alpha }}},\end{eqnarray} \tag{ 17 }$$

where $k=2\pi m/N$ with $m=1,...,N$.

4.2. Propagation from staggered initial state

We choose a staggered initial state $| 1010...\rangle$ in position space, i.e., initially every second site is occupied. For simplicity of notation we assume the number N of lattice sites to be even. A straightforward calculation, similar to that in [37], yields

$$\begin{eqnarray}&&\langle {n}_{j}(t)\rangle ={\displaystyle \frac{1}{2}}-{\displaystyle \frac{{(-1)}^{j}}{2N}}{\displaystyle \sum }_{n=1}^{N}\mathrm{cos}[t\Delta (k)]\end{eqnarray} \tag{ 18 }$$

for the time-dependence of the occupation number at lattice site j, where

$$\begin{eqnarray}&&\Delta (k):= \epsilon (k+\pi )-\epsilon (k)=2{\displaystyle \sum }_{l=1}^{N/2}{\displaystyle \frac{\mathrm{cos}[k(2l-1)]}{{d}_{2l-1}^{\alpha }}}\end{eqnarray} \tag{ 19 }$$

and $k=2\pi m/N$ with $m=1,...,N$. In figure 4 (left) the time evolution of the occupation number $\langle {n}_{j}(t)\rangle$ is plotted for different values of α, showing that the time it takes to relax to the equilibrium value of $1/2$ increases dramatically for small α (note the logarithmic timescale). This may seem counterintuitive, as a longer interaction range may naively be expected to lead to faster propagation. The effect can be understood from figure 4 (right), showing the spectrum of the frequencies Δ in the cosine terms of equation (18). As α decreases, the majority of these frequencies lie within a small window around zero, implying very slow dephasing of the cosine terms.

Zoom In Zoom Out Reset image size

A more refined picture of the propagation behavior can be obtained by studying the spreading of correlations. Starting again from a staggered initial state, a straightforward calculation similar to that in [37, 38], and similar to the one leading to (18), yields

$$\begin{eqnarray}&&\langle {c}_{j+\delta }^{\dagger }(t){c}_{j}(t)\rangle ={\displaystyle \frac{1}{2}}{\delta }_{\delta ,0}-{\displaystyle \frac{{(-1)}^{j+\delta }}{2N}}{\displaystyle \sum }_{k}{\rm{e}}^{\rm{i}t[\epsilon (k+\pi )-\epsilon (k)]}{\rm{e}}^{-\rm{i}k\delta }.\end{eqnarray} \tag{ 20 }$$

Figure 5 (bottom) shows contour plots in the $(\delta ,t)$-plane of the absolute values of the correlations (20) for different values of α. For all α shown, a cone-like propagation front is clearly visible, even in the case of $\alpha =3/4\lt D$. Two properties of the cone can be observed to change upon variation of the exponent α: (i) the boundary of the cone is rather sharp for larger α (like α = 3), whereas correlations 'leak' into the exterior of the cone for smaller α (like α = 3/4). (ii) The velocity of propagation, corresponding to the inverse slope of the cone, decreases with decreasing α (see figure 6 (left)), confirming the counterintuitive observations of figure 4 (left). We will argue in section 4.4 that some of these features can be understood on the basis of the dispersion relation (17) and the density of states of the long-range hopping model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.3. Dispersion and group velocity

In the limit of large system size the dispersion relation takes the form

$$\begin{eqnarray}&&\epsilon (k)=-\left[{\mathrm{Li}}_{\alpha }\left({\rm{e}}^{\rm{i}k}\right)+{\mathrm{Li}}_{\alpha }\left({\rm{e}}^{-\rm{i}k}\right)\right],\end{eqnarray} \tag{ 21 }$$

where ${\mathrm{Li}}_{\alpha }$ is the polylogarithm [39], and this function is plotted in figure 7 (left). For α = 3 the dispersion is a smooth function of k, while it shows a cusp at k = 0 for α = 2, and a divergence at k = 0 for α = 1. Correspondingly, the derivative ${\epsilon }^{\prime }(k)$ as shown in figure 7 (right) is discontinuous at k = 0 for α = 2, and diverges at k = 0 for α = 1. More generally we can analyze $\epsilon \prime$ in the vicinity of k = 0 by considering the difference quotient between the zeroth and the first mode

$$\begin{eqnarray}&&\begin{array}{c}| {\displaystyle \frac{\epsilon (2\pi /N)-\epsilon (0)}{2\pi (1-0)/N}}| ={\displaystyle \frac{N}{2\pi }}{\displaystyle \sum }_{l=1}^{N-1}{\displaystyle \frac{| \mathrm{cos}\left(2\pi l/N\right)-1| }{{d}_{l}^{\alpha }}}\\ \geqslant {\displaystyle \frac{N}{2\pi }}{\displaystyle \sum }_{l=1}^{N-1}{\displaystyle \frac{{(2\pi l/N)}^{2}}{{d}_{l}^{\alpha }}}={\displaystyle \frac{4\pi }{N}}{\displaystyle \sum }_{l=1}^{N/2}{l}^{2-\alpha }.\end{array}\end{eqnarray} \tag{ 22 }$$

Zoom In Zoom Out Reset image size

In the large-N limit we approximate the sum by an integral

$$\begin{eqnarray}&&{\displaystyle \frac{4\pi }{N}}{{\displaystyle \int }}_{1}^{N/2}{l}^{2-\alpha }\rm{d}l={\displaystyle \frac{2\pi }{N(3-\alpha )}}\left[{(N/2)}^{3-\alpha }-1\right]\sim {N}^{2-\alpha }.\end{eqnarray} \tag{ 23 }$$

This implies that, for $\alpha \lt 2$, the derivative $\epsilon \prime$ diverges at k = 0 in the limit of infinite system size. Interpreting ${\epsilon }^{\prime }(0)$ as a group velocity, we infer that we have a finite group velocity only for $\alpha \gt 2$, whereas the concept of a group velocity breaks down for $\alpha \lt 2$ ⁷ . This finding can help us to understand figure 5: for $\alpha \gt 2$ a finite group velocity restricts the propagation to the interior of a cone, which makes this cone appear rather sharp. For $\alpha \lt 2$, although a cone is still visible, larger (and, in fact, arbitrarily large) propagation velocities may occur and are responsible for the 'leaking' of correlations outside the cone.

The threshold value α = 2 for supersonic propagation (i.e., propagation not bounded by any finite group velocity) is also found in a different context, by very different methods. In [15] it was proved that information can be transferred supersonically through a quantum channel with finite local Hilbert space dimension for any $\alpha \lt D+1$, while no such proof exists for $\alpha \gt D+1$ ⁸ , but this result requires the measuring of observables supported on semi-infinite sublattices, which is not the most physical scenario. In appendix B we use techniques similar to those in [15], but apply them to a model with translationally invariant interactions, to prove that supersonic transmission through a quantum channel can occur for any $\alpha \lt 2$, also for measurements performed on single lattice sites. This result can be seen as complementary to the experimental observations in [28], where supersonic propagation of correlations was observed for exponents up to $\alpha \approx 1.19$ in a one-dimensional lattice.

4.4. Density of states and typical propagation velocities

From figure 5 and the discussion in section 4.3 we have seen that, while supersonic propagation can occur for $\alpha \lt 2$, cone-like propagation is observed for these values of α at least for some initial states. In this section we argue that the qualitative features of the observed behavior can be understood on the basis of the density of states

$$\begin{eqnarray}&&\rho (v)={\displaystyle \frac{1}{2\pi }}{{\displaystyle \int }}_{0}^{2\pi }\delta \left(v-{\displaystyle \frac{\rm{d}\epsilon }{\rm{d}k}}\right)\rm{d}k\end{eqnarray} \tag{ 24 }$$

in the large system limit. Equation (24) can be rewritten as

$$\begin{eqnarray}&&\rho (v)={\displaystyle \frac{1}{2\pi }}{\displaystyle \sum }_{{k}_{0}}{{\displaystyle \int }}_{0}^{2\pi }\delta \left(k-{k}_{0}\right){| {\displaystyle \frac{{\rm{d}}^{2}}{\rm{d}{k}^{2}}}\epsilon ({k}_{0})| }^{-1}\rm{d}k,\end{eqnarray} \tag{ 25 }$$

where the sum is taken over all roots k₀ of the argument of the delta function. The polylogarithms that appear in the dispersion relation (21) can be analytically evaluated for certain integer values of α, yielding

$$\begin{eqnarray}\rho (v)={\displaystyle \frac{1}{\pi }}\left\{\begin{array}{cc}{\displaystyle \frac{1}{1+{v}^{2}}} & \mathrm{for}\ \alpha =1,\\ {\displaystyle \frac{1}{2}}\Theta (\pi -v)\Theta (\pi +v) & \mathrm{for}\ \alpha =2,\\ {\displaystyle \frac{1}{\sqrt{4-{v}^{2}}}} & \mathrm{for}\ \alpha \to \infty ,\end{array}\right.\end{eqnarray} \tag{ 26 }$$

where Θ is the Heaviside step function. For those three values of α, the density of states is plotted in figure 6 (right), but other cases can be evaluated numerically (not shown in the figure). Again, as for the group velocity in figure 7 and the classical information capacity in appendix B, we find a threshold value of α = 2, as explained in the following.

For $\alpha \lt 2$, the density of states ρ is nonzero for all v, implying that propagation is not bounded by any finite maximum velocity. The maximum of ρ, however, is at v = 0 for all $\alpha \lt 2$, and this gives an indication that slow propagation with a small velocity is favored, although larger velocities do occur (as in figure 5 (left and center)). The maximum at v = 0 becomes more sharply peaked when α approached zero, explaining the vanishing of the inverse slope of the cone in figure 5 in that limit, as shown in figure 6 (left).

For $\alpha \geqslant 2$, the density of states ρ is nonzero only on a finite interval $[-{v}_{\mathrm{max}},+{v}_{\mathrm{max}}]$, where v_max depends on α. For $\alpha \gt 2$ the density of states diverges, and therefore takes on its maximum, at $\pm {v}_{\mathrm{max}}$. This implies that the maximum velocity is favored, although smaller velocities also occur (as in figure 5 (right)).

In this paper we have studied, from several different perspectives, the nonequilibrium dynamics of lattice models with long-range interactions or long-range hopping, and in particular the propagation in space and time of correlations and other physical quantities. The focus of our work is on the competition between linear, cone-like propagation and faster-than-linear, supersonic propagation. We illustrate this competition in two regimes, both relevant for experimental realizations of long-range many-body systems in cold atoms, ions, or molecules:

• (i)  
  For small exponents $\alpha \lt 2$ we prove that supersonic propagation can occur. At the same time, in such systems cone-like spreading can be the dominant form of propagation, with supersonic effects appearing only as small corrections (as in figure 5 (center)).
• (ii)  
  For intermediate exponents (roughly between 3 and 8), propagation is observed to be linear initially, with supersonic effects setting in at larger times and distances (as in figures 1 (center) and 3 (right)).

To explain these observations, we provide model calculations as well as general bounds that provide a comprehensive and consistent picture of the various shapes of propagation fronts that can occur. Two of our results are Lieb–Robinson-type bounds, valid for large classes of models with long-range interactions. The first is a bound for models with exponents α smaller than the lattice dimension D, a regime for which hitherto no such bounds existed. Key to deriving the bound is the insight that for $\alpha \lt D$ the propagation speed in general scales asymptotically like a power law with the system size, and a meaningful bound therefore has to be derived in rescaled time τ as defined in (9). In physical time t, the bound then describes the increase of the propagation speed with increasing lattice size, as illustrated in figure 2. The second Lieb–Robinson-type bound we report is essentially a cheat, as we stop half way through the derivation of a 'conventional' Lieb–Robinson bound. Specializing this result to single-site observables and Hamiltonians with pair interactions only, we obtain an expression that can be evaluated numerically in an efficient way, easily reaching system sizes of $\mathcal{O}({10}^{4})$. This bound (11) is sharp enough to capture cone-like as well as supersonic behavior. In experimental studies of propagation in long-range interacting lattice models [28, 41], the currently feasible lattice sizes are small and measured data can be compared to results from exact diagonalization. However, experimental work on systems of larger size is in progress, and exact diagonalization will not be feasible in that case. We expect that the matrix exponential bound (11) can provide guidance and sanity checks when analyzing the results of such experiments.

In the second half of the paper we complemented the bounds with results of one of the simplest long-range quantum models, namely a fermionic long-range hopping model in one-dimension. We observed that cone-like propagation fronts can be a dominant feature also for small values of α, and we explain the opening angle of such a cone, as well as the interplay of cone-like and supersonic features, on the basis of the dispersion relation combined with the density of states. These results indicate that it will depend crucially on the k-modes occupied whether cone-like or supersonic propagation is dominant. We expect that such an improved understanding can provide guidance for optimizing experimental efforts to harness long-range interactions in a variety of quantum information and signaling tasks.

The authors acknowledge helpful discussions with J Eisert, F Essler, M Foss-Feig, A Gorshkov, S Kehrein, S Manmana, R Sweke and D Vodola. D S acknowledges financial support by the Studienstiftung des deutschen Volkes; M K by the National Research Foundation of South Africa through the Incentive Funding and the Competitive Programme for Rated Researchers.

A.1. Derivation of the matrix exponential bound

As in section 2, we consider a D-dimensional lattice Λ consisting of N sites, a Hilbert space (2) consisting of finite-dimensional local Hilbert spaces, and a generic Hamiltonian with pair interactions (10). Let O_A and O_B be two bounded linear operators compactly supported on subsets $A,B\subset \Lambda$ with $A\cap B=\varnothing$. Under these conditions, similar to the derivation of equation (2.12) of [25], one can derive the upper bound

$$\begin{eqnarray}&&\parallel \left[{O}_{A}(t),{O}_{B}(0)\right]\parallel \leqslant 2\parallel {O}_{A}\parallel \parallel {O}_{B}\parallel {\displaystyle \sum }_{n=1}^{\infty }{\displaystyle \frac{{(2| t| )}^{n}}{n!}}{a}_{n}.\end{eqnarray} \tag{ A.1 }$$

For pair interactions, and considering observables O_A and O_B supported on single lattice sites only (i.e., $A=\{i\}$ and $B=\{j\}$), the coefficients a_n are upper bounded by

$$\begin{eqnarray}&&{a}_{n}\leqslant {\kappa }^{n}{\displaystyle \sum }_{{k}_{1},...,{k}_{n-1}}{J}_{i,{k}_{1}}{J}_{{k}_{1},{k}_{2}}...{J}_{{k}_{n-1},j}={\kappa }^{n}{\left({J}^{n}\right)}_{i,j},\end{eqnarray} \tag{ A.2 }$$

where J is the interaction matrix with elements ${J}_{k,l}=\parallel {h}_{{kl}}\parallel$ and $\kappa ={\mathrm{sup}}_{q\in \Lambda }{{\displaystyle \sum }}_{k}{J}_{q,k}$. Then (A.1) can be written as

$$\begin{eqnarray}&&{\displaystyle \frac{\parallel \left[{O}_{A}(t),{O}_{B}(0)\right]\parallel }{2\parallel {O}_{A}\parallel \parallel {O}_{B}\parallel }}\leqslant {\left({\displaystyle \sum }_{n=1}^{\infty }{\displaystyle \frac{{(2\kappa | t| )}^{n}}{n!}}{J}^{n}\right)}_{i,j}=\mathrm{exp}{(2J\kappa | t| )}_{i,j}-{\delta }_{i,j},\end{eqnarray} \tag{ A.3 }$$

which proves (11).

For translationally invariant one-dimensional lattices, J is a circulant matrix and can be diagonalized by means of Fourier transformation. For the example of power law interactions ${h}_{{kl}}\propto {\text{d}}{(k,l)}^{-\alpha }$, the diagonal elements of the Fourier-transformed matrix J are given by

$$\begin{eqnarray}&&\epsilon (k)={\displaystyle \sum }_{n=1}^{\lfloor (N-1)/2\rfloor }{\displaystyle \frac{\mathrm{cos}({nk})}{{n}^{\alpha }}}\underset{\mathrm{if}\ N\ \mathrm{even}}{\underbrace{+{\displaystyle \frac{{\rm{e}}^{{ikN}/2}}{2{(N/2)}^{\alpha }}}}}\end{eqnarray} \tag{ A.4 }$$

with $k=2\pi m/N$, $0\leqslant m\lt N$. Using these eigenvalues, J can be exponentiated in the diagonal basis, followed by an inverse Fourier transformation to evaluate the Lieb–Robinson bound (11).

We envisage the bound (11) to be particularly useful for finite systems of intermediate size where the matrix exponential can be computed numerically with reasonable effort. However, since (11) is sharper than the bounds in [23, 25], a thermodynamic limit will exist (at least) under the same conditions required in those proofs, and in particular for D-dimensional regular lattices with power law interactions and exponents $\alpha \gt D$.

A.2. Lieb–Robinson bounds in rescaled time

As in section 3, we consider a D-dimensional lattice Λ consisting of N sites, a Hilbert space (2) with finite-dimensional local Hilbert spaces, and a generic Hamiltonian H with n-body interactions (3). We require that H satisfies conditions (4) and (5). For the proof of a Lieb–Robinson-type bound, we follow the general strategy of [25], augmented with the ${\mathcal{N}}_{\Lambda }$-rescaling taken from [33].

As a shorthand we introduce

$$\begin{eqnarray}&&\chi (t)=\left[{O}_{A}(t),{O}_{B}\right].\end{eqnarray} \tag{ A.5 }$$

for the commutator on the left-hand side of (8). Differentiating with respect to t yields

$$\begin{eqnarray}&&{\chi }^{\prime }(t)=\rm{i}\left[\chi (t),{I}_{A}(t)\right]+\rm{i}\left[{O}_{A}(t),\left[{I}_{A}(t),{O}_{B}\right]\right],\end{eqnarray} \tag{ A.6 }$$

where ${I}_{A}={{\displaystyle \sum }}_{Z:Z\cap A\ne \varnothing }{h}_{Z}$ is the set of local Hamiltonian terms that have non-zero overlap with A. Using the boundedness of O_A(t) we apply lemma A.1 of [25] to the norm-preserving first term of (A.6), yielding

$$\begin{eqnarray}&&\parallel \chi (t)\parallel -\parallel \chi (0)\parallel \leqslant 2\parallel {O}_{A}\parallel {\displaystyle \sum }_{Z:Z\cap A\ne \varnothing }{{\displaystyle \int }}_{0}^{| t| }\parallel \left[{h}_{Z}(s),{O}_{B}\right]\parallel \rm{d}s.\end{eqnarray} \tag{ A.7 }$$

Next we define

$$\begin{eqnarray}&&{C}_{{O}_{B}}(A,t):= \underset{{O}_{A}\in {\mathcal{O}}_{A}}{\mathrm{sup}}{\displaystyle \frac{\parallel \chi (t)\parallel }{\parallel {O}_{A}\parallel }},\end{eqnarray} \tag{ A.8 }$$

where ${\mathcal{O}}_{A}$ is the set of observables compactly supported on A. Making use of this definition, (A.7) can be rewritten as

$$\begin{eqnarray}&&{\displaystyle \frac{{C}_{{O}_{B}}(A,t)-{C}_{{O}_{B}}(0,t)}{2}}\leqslant {\displaystyle \sum }_{Z:Z\cap A\ne \varnothing }{{\displaystyle \int }}_{0}^{| t| }\parallel {h}_{Z}\parallel {C}_{{O}_{B}}(Z,s)\rm{d}s.\end{eqnarray} \tag{ A.9 }$$

Equation (A.9) can be applied recursively to show that

$$\begin{eqnarray}&&{C}_{{O}_{B}}(A,t)\leqslant 2\parallel {O}_{B}\parallel {\displaystyle \sum }_{n=1}^{\infty }{\displaystyle \frac{{(2| t| )}^{n}}{n!}}{a}_{n}\end{eqnarray} \tag{ A.10 }$$

with coefficients

$$\begin{eqnarray}&&{a}_{n}={\displaystyle \sum }_{{\displaystyle \genfrac{}{}{0em}{}{{Z}_{1}\subset \Lambda }{{Z}_{1}\cap A\ne \varnothing }}}{\displaystyle \sum }_{{\displaystyle \genfrac{}{}{0em}{}{{Z}_{2}\subset \Lambda }{{Z}_{2}\cap {Z}_{1}\ne \varnothing }}}...{\displaystyle \sum }_{{\displaystyle \genfrac{}{}{0em}{}{{Z}_{n}\subset \Lambda }{{Z}_{n}\cap {Z}_{n-1}\ne \varnothing }}}{\displaystyle \prod }_{l=1}^{n}\ \parallel {h}_{{Z}_{l}}\parallel {\delta }_{B}({Z}_{n}),\end{eqnarray} \tag{ A.11 }$$

where

$$\begin{eqnarray}{\delta }_{B}(Z)=\{\begin{array}{cc}0\quad & \mathrm{if}Z\cap B\ne \varnothing ,\\ 1\quad & \mathrm{otherwise}.\end{array}\end{eqnarray} \tag{ A.12 }$$

Under the conditions (4) and (5) these coefficients can be bounded by

$$\begin{eqnarray}&&{a}_{n}\leqslant {\displaystyle \frac{{p}^{n-1}{\lambda }^{n}}{{\mathcal{N}}_{\Lambda }^{n}{(1+{\text{d}}(A,B))}^{\alpha }}}.\end{eqnarray} \tag{ A.13 }$$

Inserting (A.13) into (A.10) and using the definition (9) of rescaled time τ, one obtains

$$\begin{eqnarray}&&{C}_{{O}_{B}}(A,t)\leqslant {\displaystyle \frac{2\parallel {O}_{B}\parallel | A| | B| }{p{(1+{\text{d}}(A,B))}^{\alpha }}}(\mathrm{exp}[2\;p\lambda | \tau | ]-1),\end{eqnarray} \tag{ A.14 }$$

and this implies the bound

$$\begin{eqnarray}&&\parallel \left[{O}_{A}(\tau {\mathcal{N}}_{\Lambda }),{O}_{B}\right]\parallel \leqslant {\displaystyle \frac{2\parallel {O}_{A}\parallel \parallel {O}_{B}\parallel | A| | B| }{p{(1+{\text{d}}(A,B))}^{\alpha }}}(\mathrm{exp}[2p\lambda | \tau | ]-1)\end{eqnarray} \tag{ A.15 }$$

in rescaled time τ, valid for power law interactions with exponents $\alpha \gt 0$.

A.3. Discussion of the bound in [12]

In [26] a Lieb–Robinson-type bound was derived whose functional form consists of a linear (cone-like) and a faster-than-linear (supersonic) contribution. This bound is a major improvement over that in [23] in the regime of large α, where the former becomes more and more similar to a nearest-neighbor bound, as it should. Here we scrutinize the applicability of the bound in [26] for describing the cone-like and supersonic features of long-range models with intermediate exponents α (roughly in the range $3\leqslant \alpha \leqslant 8$).

The bound in [26] is derived for Hamiltonians

$$\begin{eqnarray}&&H={\displaystyle \frac{1}{2}}{\displaystyle \sum }_{i\ne j}{h}_{{ij}}\end{eqnarray} \tag{ A.16 }$$

with two-body interactions h_ij satisfying

$$\begin{eqnarray}&&\parallel {h}_{{ij}}\parallel \leqslant {\displaystyle \frac{1}{{\text{d}}{(i,j)}^{\alpha }}}\end{eqnarray} \tag{ A.17 }$$

on D-dimensional regular cubic lattices. For exponents $\alpha \geqslant 1$ a bound of the form

$$\begin{eqnarray}&&{\displaystyle \frac{\parallel [A(t),B]\parallel }{2| A| | B| }}\leqslant {T}_{1}+{T}_{2}\end{eqnarray} \tag{ A.18 }$$

is obtained, where A and B are observables on lattice sites that are a distance δ apart, and

$$\begin{eqnarray}&&{T}_{1}={c}_{1}{\displaystyle \frac{{e}^{{v}_{1}t}-1}{{e}^{\mu \delta }}},\qquad {T}_{2}={c}_{2}{\displaystyle \frac{{e}^{{v}_{2}t}-1}{{[(1-\mu )\delta ]}^{\alpha }}},\end{eqnarray} \tag{ A.19 }$$

with ${c}_{1}={\lambda }^{-1}$, ${v}_{1}=2{\lambda }^{2}e$, ${c}_{2}={(\lambda {9}^{D})}^{-1}$, ${v}_{2}=2{\lambda }^{2}{9}^{D}$, $\lambda ={{\displaystyle \sum }}_{k}d{(i,k)}^{-\alpha }$, and $0\lt \mu \lt 1$. T₁ has the same functional form as the classic Lieb–Robinson bound for Hamiltonians with finite-range interactions [22], which is known to produce a linear, cone-shaped causal region. T₂ has the functional form of the bound originally derived by Hastings and Koma [23]. Both, T₁ and T₂ contain the free parameter μ, which determines, among other things, the slope of the linear soundcone. So the 'velocity' associated with the cone can be tuned to an arbitrary value, irrespectively of the physical behavior of the model studied.

Based on (A.18) and (A.19), the sharpest bound

$$\begin{eqnarray}&&B(\delta ,t)=\underset{\mu }{\mathrm{min}}\left({c}_{1}{\displaystyle \frac{{e}^{{v}_{1}t}-1}{{e}^{\mu \delta }}}+{c}_{2}{\displaystyle \frac{{e}^{{v}_{2}t}-1}{{[(1-\mu )\delta ]}^{\alpha }}}\right)\end{eqnarray} \tag{ A.20 }$$

is obtained by minimizing, separately for each value of δ and t, over the free parameter μ. From the contour plots of B in figure A1 it becomes clear that the 'linearity' of T₁ can be deceiving, as a linear, cone-like regime is not particularly prominent, not even for larger exponents like $\alpha =8$. Of course it is always possible to construct a linear-looking bound by weakening B, but this would be unrelated to the physical behavior of the class of models studied.

Zoom In Zoom Out Reset image size

Another, more sophisticated bound has recently been put forward in [27], but the form of the propagation front has not yet been analyzed and discussed (beyond the long-distance asymptotics).

In this appendix we prove that supersonic transmission through a quantum channel can occur for any $\alpha \lt 2$, also for measurements performed on single lattice sites. Like for the study of the group velocity of the long-range hopping model in section 4.3, we find a threshold value of $\alpha =2$ below which propagation becomes supersonic. The proof uses techniques from [15] and applies them to a slightly more involved model for which supersonic propagation is found to occur also for single-site measurements.

We consider a finite one-dimensional lattice $\Lambda =\{1,...,N\}$ consisting of N sites. To implement a quantum channel, we encode a signal on site 1, and measure the effect of that encoding after a time time at site N. On this lattice we define an Ising Hamiltonian with arbitrary couplings,

$$\begin{eqnarray}&&H={\displaystyle \sum }_{i\lt j}{J}_{{ij}}{\sigma }_{i}^{z}{\sigma }_{j}^{z}.\end{eqnarray} \tag{ B.1 }$$

Defining the sublattices $A=\{1\}$, $B=\{N\}$, and $S=\Lambda \setminus (A\cup B)$, the Hamiltonian can be rewritten as

$$\begin{eqnarray}&&H={H}_{{AS}}+{H}_{{AB}}+{H}_{{SB}}+{H}_{{SS}}\end{eqnarray} \tag{ B.2 }$$

with

$$\begin{eqnarray}&&{H}_{{XY}};= {\displaystyle \sum }_{i\in X}{\displaystyle \sum }_{j\in Y}{J}_{{ij}}{\sigma }_{i}^{z}{\sigma }_{j}^{z},\end{eqnarray} \tag{ B.3 }$$

where $X,Y\in \{A,S,B\}$. As an initial state we choose

$$\begin{eqnarray}&&\rho (0)={| \downarrow \rangle }_{1}{\;}_{1}\langle \downarrow | \underset{s\in S}{{\displaystyle \otimes }}{| \downarrow \rangle }_{s}{\;}_{s}\langle \downarrow | \otimes {| +\rangle }_{N}{\;}_{N}\langle +| \end{eqnarray} \tag{ B.4 }$$

with ${\sigma }_{i}^{z}{| \uparrow \rangle }_{i}={| \uparrow \rangle }_{i}$, ${\sigma }_{i}^{z}{| \downarrow \rangle }_{i}=-{| \downarrow \rangle }_{i}$, and ${| +\rangle }_{j}=\left({| \uparrow \rangle }_{j}+{| \downarrow \rangle }_{j}\right)/\sqrt{2}$. Initially all the spins are pointing down, except the one at $B=\{N\}$.

A binary quantum channel is implemented by starting the time evolution either with $\rho (0)$ (sending a '0'), or starting with ${U}_{A}\rho (0){U}_{A}^{\dagger }$ (sending a '1'), where U_A is a unitary supported on A only. The classical information capacity C_t can be bounded from below by the probability to detect, by measuring according to a positive operator valued measure ${\pi }_{B}$, a signal at B after a time t,

$$\begin{eqnarray}&&{C}_{t}\geqslant {p}_{t}=| \mathrm{Tr}\left\{{N}_{t}[\rho (0)]{\pi }_{B}\right\}-\mathrm{Tr}\left\{{T}_{t}[\rho (0)]{\pi }_{B}\right\}| ,\end{eqnarray} \tag{ B.5 }$$

with

$$\begin{eqnarray}&&{N}_{t}[\rho (0)];= {\mathrm{Tr}}_{\Lambda \setminus B}\left[{\rm{e}}^{-\rm{i}{Ht}}\rho (0){\rm{e}}^{\rm{i}{Ht}}\right],\end{eqnarray} \tag{ B.6 }$$

$$\begin{eqnarray}&&{T}_{t}[\rho (0)];= {\mathrm{Tr}}_{\Lambda \setminus B}\left[{\rm{e}}^{-\rm{i}{Ht}}{U}_{A}\rho (0){U}_{A}^{\dagger }{\rm{e}}^{\rm{i}{Ht}}\right].\end{eqnarray} \tag{ B.7 }$$

In the following we compute a lower bound on the right-hand side of (B.5), and study this bound as a function of the channel length, i.e., the distance between A and B.

We choose ${\pi }_{B}=| +{\rangle }_{N}{\;}_{N}\langle +|$ and ${U}_{A}=| \uparrow {\rangle }_{1}{\;}_{1}\langle \downarrow |$, where the latter is a spin flip operator on the first lattice site. For the time-evolved density operator in (B.6) we find

$$\begin{eqnarray}&&\rho (t)={\rm{e}}^{-\rm{i}{H}_{{AB}}t}{\rm{e}}^{-\rm{i}{H}_{{AS}}t}{\rm{e}}^{-\rm{i}{H}_{{SS}}t}{\rm{e}}^{-\rm{i}{H}_{{SB}}t}\end{eqnarray} \tag{ B.8 }$$

$$\begin{eqnarray*}&&\begin{array}{c}\times \left[\underset{s=1}{\overset{N-1}{{\displaystyle \otimes }}}{| \downarrow \rangle }_{s}{\;}_{s}\langle \downarrow | \;\otimes \;{| +\rangle }_{N}{\;}_{N}\langle +| \right]{\rm{e}}^{\rm{i}{H}_{{SB}}t}{\rm{e}}^{\rm{i}{H}_{{SS}}t}{\rm{e}}^{\rm{i}{H}_{{AS}}t}{\rm{e}}^{\rm{i}{H}_{{AB}}t}\\ =\underset{s=1}{\overset{N-1}{{\displaystyle \otimes }}}{| \downarrow \rangle }_{s}{\;}_{s}\langle \downarrow | \left\{\mathrm{exp}\left[\rm{i}t{\displaystyle \sum }_{r=1}^{N-1}{J}_{{rN}}{\sigma }_{N}^{z}\right]{| +\rangle }_{N}{\;}_{N}\langle +| \mathrm{exp}\left[-\rm{i}t{\displaystyle \sum }_{r=1}^{N-1}{J}_{{rN}}{\sigma }_{N}^{z}\right]\right\}.\end{array}\end{eqnarray*}$$

All the exponentials not supported on B add up to zero since the initial state prepared on $\Lambda \setminus B$ is an eigenstate of the Ising Hamiltonian. Taking the trace gives

$$\begin{eqnarray}&&\mathrm{Tr}\left\{{N}_{t}[\rho (0)]{\pi }_{B}\right\}={\displaystyle \frac{1}{2}}\left\{1+\mathrm{cos}\left[2t\left({\displaystyle \sum }_{r\in S}{J}_{{rN}}+{J}_{1N}\right)\right]\right\}.\end{eqnarray} \tag{ B.9 }$$

A similar calculation shows that

$$\begin{eqnarray}&&\mathrm{Tr}\left\{{T}_{t}[\rho (0)]{\pi }_{B}\right\}={\displaystyle \frac{1}{2}}\left\{1+\mathrm{cos}\left[2t\left({\displaystyle \sum }_{r\in S}{J}_{{rN}}-{J}_{1N}\right)\right]\right\}.\end{eqnarray} \tag{ B.10 }$$

The probability of detecting a signal in B at some time $t\gt 0$ is then given by

$$\begin{eqnarray}&&{p}_{t}={\displaystyle \frac{1}{2}}| \mathrm{cos}\left[2t\left({\displaystyle \sum }_{r\in S}{J}_{{rN}}+{J}_{1N}\right)\right]-\mathrm{cos}\left[2t\left({\displaystyle \sum }_{r\in S}{J}_{{rN}}-{J}_{1N}\right)\right]| \end{eqnarray} \tag{ B.11 }$$

$$\begin{eqnarray}&&=| \mathrm{sin}\left(2t{\displaystyle \sum }_{r\in S}{J}_{{rN}}\right)\mathrm{sin}\left(2{{tJ}}_{1N}\right)| .\end{eqnarray} \tag{ B.12 }$$

To derive a nontrivial (nonzero) lower bound on p_t, we target the regime before oscillatory behavior in (B.12) sets in. Using the inequality

$$\begin{eqnarray}&&\mathrm{sin}x\geqslant {\displaystyle \frac{2x}{\pi }}\qquad \mathrm{for}\ 0\leqslant x\leqslant \pi /2\end{eqnarray} \tag{ B.13 }$$

and assuming power law interactions ${J}_{{ij}}=| i-j{| }^{-\alpha }$, we obtain

$$\begin{eqnarray}&&{p}_{t}\geqslant {\displaystyle \frac{4t}{\pi }}{\displaystyle \frac{1}{{(N-1)}^{\alpha }}}{\displaystyle \frac{4t}{\pi }}{\displaystyle \sum }_{r=2}^{N-1}{\displaystyle \frac{1}{{(N-r)}^{\alpha }}},\end{eqnarray} \tag{ B.14 }$$

valid for times

$$\begin{eqnarray}&&t\leqslant {\displaystyle \frac{\pi }{4}}/{\displaystyle \sum }_{r=1}^{N-2}{\displaystyle \frac{1}{{r}^{\alpha }}}.\end{eqnarray} \tag{ B.15 }$$

Interpreting the sum in (B.14) as an upper Riemann sum, we have

$$\begin{eqnarray}&&{\displaystyle \sum }_{r=2}^{N-1}{\displaystyle \frac{1}{{(N-r)}^{\alpha }}}={\displaystyle \sum }_{r=1}^{N-2}{\displaystyle \frac{1}{{r}^{\alpha }}}\gt {{\displaystyle \int }}_{0}^{N-2}{\displaystyle \frac{\rm{d}r}{{(r+1)}^{\alpha }}}.\end{eqnarray} \tag{ B.16 }$$

Then we can bound p_t by

$$\begin{eqnarray}&&{p}_{t}\gt {\displaystyle \frac{16{t}^{2}}{{\pi }^{2}(\alpha -1)}}{\displaystyle \frac{1}{{(N-1)}^{\alpha }}}\left(1-{\displaystyle \frac{1}{{(N-1)}^{\alpha -1}}}\right)\;:=\;{\underline{p}}_{t}.\end{eqnarray} \tag{ B.17 }$$

For $\alpha \gt 1$ and large N the second term in the square bracket in (B.17) is much smaller than 1, and we obtain

$$\begin{eqnarray}&&{\underline{p}}_{t}\sim {\displaystyle \frac{16{t}^{2}}{{\pi }^{2}(\alpha -1)}}{\displaystyle \frac{1}{{(N-1)}^{\alpha }}}\end{eqnarray} \tag{ B.18 }$$

for the large-N asymptotic behavior of the bound ${\underline{p}}_{t}$. In our setting, $\delta =N-1$ is the distance between the regions A and B. To determine the shape of a contour line at which ${\underline{p}}_{t}$ is equal to some constant , we set

$$\begin{eqnarray}&&\epsilon ={\underline{p}}_{t}\propto {\displaystyle \frac{{t}^{2}}{{\delta }^{\alpha }}},\end{eqnarray} \tag{ B.19 }$$

and we can read off that

$$\begin{eqnarray}&&\delta \propto {t}^{2/\alpha }\end{eqnarray} \tag{ B.20 }$$

along any of those contour lines. Equation (B.20) describes faster-than-linear (supersonic) growth of δ for $\alpha \lt 2$. It is straightforward to extend the above calculation to more general initial conditions as well as to lattices of arbitrary dimension.