Spreading of Perturbations in LongRange Interacting Classical Lattice Models
Abstract
LiebRobinsontype bounds are reported for a large class of classical Hamiltonian lattice models. By a suitable rescaling of energy or time, such bounds can be constructed for interactions of arbitrarily long range. The bound quantifies the dependence of the system’s dynamics on a perturbation of the initial state. The effect of the perturbation is found to be effectively restricted to the interior of a causal region of logarithmic shape, with only small, algebraically decaying effects in the exterior. A refined bound, sharper than conventional LiebRobinson bounds, is required to correctly capture the shape of the causal region, as confirmed by numerical results for classical longrange chains. We discuss the relevance of our findings for the relaxation to equilibrium of longrange interacting lattice models.
In many nonrelativistic lattice systems, and despite the absence of Lorentz covariance, physical effects are mostly restricted to a causal region, often in the shape of an effective “light cone,” with only tiny effects leaking out to the exterior. The technical tool, known as LiebRobinson bounds Lieb and Robinson (1972); Nachtergaele and Sims (2010); *Review, to quantify this statement in a quantum mechanical context is an upper bound on the norm of the commutator , where and are operators supported on the subspaces of the Hilbert space corresponding to nonoverlapping regions and of the lattice. The importance of such a bound lies in the fact that a multitude of physically relevant results can be derived from it. Examples are bounds on the creation of equaltime correlations Nachtergaele et al. (2006), on the transmission of information Bravyi et al. (2006), and on the growth of entanglement Eisert and Osborne (2006), the exponential spatial decay of correlations in the ground state of a gapped system Hastings and Koma (2006), or a LiebSchultzMattis theorem in higher dimensions Hastings (2004). Experimental observations related to LiebRobinson bounds have also been reported Cheneau et al. (2012); *Langen_etal13.
The original proof by Lieb and Robinson Lieb and Robinson (1972) requires interactions of finite range. An extension to powerlawdecaying longrange interactions has been reported in Refs. Hastings and Koma (2006); Nachtergaele et al. (2006). In this case the effective causal region is no longer cone shaped, and the spatial propagation of physical effects is not limited by a finite group velocity Eisert et al. (2013). For “strong longrange interactions,” i.e., when the interaction potential decays proportionally to with an exponent smaller than the lattice dimension , the theorems in Hastings and Koma (2006); Nachtergaele et al. (2006) do not apply and no LiebRobinsontype results are known. This fact nicely fits into the larger picture that, for , the behavior of longrange interacting systems often differs substantially from that of shortrange interacting systems. Examples of such differences include nonequivalent equilibrium statistical ensembles and negative response functions LyndenBell and Wood (1968); *Thirring70; *ElHaTur00; *TouElTur04; *Kastner10; *KastnerJSTAT10, or the occurrence of quasistationary states whose lifetimes diverge with the system size Antoni and Ruffo (1995); *Kastner11; *Kastner12; Campa et al. (2009). The latter is a dynamical phenomenon, and it has been conjectured in Bachelard and Kastner (2013) that some of its properties are universal and in some way connected to LiebRobinson bounds.
In most cases the peculiarities of longrange interacting systems have been investigated in the framework of classical Hamiltonian systems Campa et al. (2009), but little is known about LiebRobinson bounds in classical mechanics. Exceptions are restricted to specific models with nearestneighbor interactions Marchioro et al. (1978); Buttà et al. (2007); *RazSims09. In the context of classical Hamiltonian mechanics, a LiebRobinson bound is an upper bound on the norm of the Poisson bracket , where and are phase space functions supported only on the subspaces of phase space corresponding to the nonoverlapping regions and of the lattice, respectively. The physical meaning of the norm of this Poisson bracket becomes evident from an expression put forward in Marchioro et al. (1978),
(1) 
where and are the (bounded) maxima of all the partial derivatives of and with respect to the phase space coordinates, and
(2) 
The partial derivatives on the righthand side of (2) quantify the effect that a variation of the initial momentum or position , at the lattice site have on the timeevolved momentum or position , at the lattice site . A classical LiebRobinson bound is therefore a measure for the spreading in time and space of an initial perturbation, with potential applications to a broad range of physical processes, including heat conduction, signal transmission, transfer of energy, or the approach to equilibrium.
In this Letter we study the spreading in time and space of initial perturbations in classical longrange interacting lattice models. We have proved a LiebRobinsontype result, providing an upper bound on in (2) [and hence on the Poisson bracket (1)] for a broad class of classical longrange interacting lattice models in arbitrary spatial dimension. Dipolar interactions in condensed matter systems are the prime example of such longrange interacting lattice systems Miloshevich et al. (2013), but many other examples exist Bachelard et al. (2010). To avoid the rather technical notation of the general result ^{1}^{1}1See Supplemental Material., we present the main result in this Letter for a specific class of systems, namely classical models in spatial dimensions with pair interactions that decay like a power law with the (1norm) distance between lattice sites and . The value of the exponent determines the range of the interaction, from meanfieldtype (distanceindependent) interactions at to nearestneighbor couplings in the limit .
Our study focuses on the influence of the interaction range on the spreading of perturbations, and we find pronounced quantitative and qualitative changes upon variation of . Different from systems with finiterange interactions, the effect of an initial perturbation is found to be effectively restricted to the interior of a causal region of logarithmic shape, with algebraically small effects in the exterior. Similar to the shortrange case, such a bound can be used to rigorously control finitesize effects in simulations of lattice models, exclude information transmission above a certain measurement resolution in the exterior of the effective causal region, and much more. Our analytical results are supplemented by numerical simulations of the time evolution of a longrange interacting chain. Besides confirming the validity of the bound, the numerical results reveal that the refined version of our bound, sharper than conventional LiebRobinsontype bounds, is required in order to correctly capture the shape of the propagation front.
model.—This model consists of classical spins (or rotors) attached to the sites of a dimensional hypercubic lattice . The phase space of a single rotor is , allowing to parametrize each rotor by an angular variable and by its angular momentum . On the phase space of the total system we define the Hamiltonian function
(3) 
For , the second sum on the righthand side of (3) is superextensive, i.e., asymptotically for large lattices it grows faster than linearly with the number of lattice sites. Our proof of a LiebRobinson bound requires the Hamiltonian to be extensive. We enforce extensivity also for by allowing the coupling constant to depend explicitly on the lattice,
(4) 
where is a real constant ^{2}^{2}2Introducing such a normalization is actually a harmless procedure and amounts to a rescaling of time. For a “physical” Hamiltonian without such a rescaling, the bounds we derived hold therefore in rescaled time . For finite systems, no rescaling is necessary..
Classical longrange LiebRobinson bound.—Upper bounds on the partial derivatives on the righthand side of (2) are given by
(5a)  
(5b)  
(5c)  
(5d) 
The positive coefficients are defined recursively by
(6) 
with and , and we use the constants
(7)  
(8) 
The proof of the bounds combines techniques for classical lattices with nearestneighbor interactions Marchioro et al. (1978) with those used for proving LiebRobinson bounds for longrange quantum systems Hastings and Koma (2006); Nachtergaele et al. (2006), with the additional refinement of allowing latticedependent coupling constants. In the thermodynamic limit of infinite lattice size, the quantity in (7) remains finite and hence the bounds remain meaningful Note1 ().
All the bounds on the righthand side of (5a)–(5d), and therefore also the norm of the Poisson bracket (1), grow exponentially for large times , and decay as a power law with the distance . For some , the effect of a perturbation is therefore smaller than outside a region in the plane specified, for large , by
(9) 
This effective causal region, in which the effect of an initial perturbation is nonnegligible, has a logarithmic shape, and differs in this respect from the linear (coneshaped) region derived in Marchioro et al. (1978) for shortrange interactions. As a consequence, no finite group velocity limits the spreading of perturbations in longrange interacting lattices, and supersonic propagation can occur. For a refined understanding of the spatiotemporal behavior, we go beyond the usual LiebRobinsontype estimates and consider the sharper bounds in (5a)–(5d), where the coefficients introduce an additional spatial dependence. The functional form of these bounds will be illustrated below, and we also show that only the sharper bounds correctly capture the shape of the propagation front.
The bounds (5a)–(5d) remain valid, with only minor modifications of the parameters and prefactors, under broad generalizations, including arbitrary graphs , multidimensional singleparticle phase spaces , manyparticle interactions, and very general forms of the interaction potential Note1 (). In the remainder of this Letter we subject the bounds (5a)–(5d) to a reality check, in the sense of testing their tightness and whether the form of the propagation front obtained numerically is faithfully reproduced.
Numerics.—We consider the model (3) on a ring, i.e., a onedimensional chain of sites with periodic boundary conditions. The partial derivatives in (5a)–(5d) are approximated by difference quotients,
(10) 
and similarly for the other derivatives. Here, is the timeevolved momentum obtained by starting from a certain initial condition, and is for a similar initial condition, but with the th initial position shifted by some small . The timeevolved momenta and are obtained by numerically integrating Hamilton’s equations using a sixthorder symplectic integrator McLachlan and Atela (1992). The numerical results for fluctuate strongly in time, obscuring the overall trend of the spreading. To reduce the fluctuating background, we compute the difference quotient (10) for 20 different (pairs of) initial conditions; details regarding the choice of initial conditions are given in the Supplemental Material Note1 (). Since our aim is to compare the numerical results to the upper bounds (5a)–(5d), we select, for any fixed time and lattice sites and , the largest of the 20 values. The resulting maximum is denoted by , and its time and distancedependence is shown in Fig. 1. The plots illustrate the supersonic propagation of perturbations in the presence of longrange interactions, as expected from the inequality (9).
For all distances and times , the numerical results are smaller than the bounds (5a)–(5d) and hence confirm their validity. What is more, the results nicely agree with the functional forms of the bounds, and only the prefactors are overestimated. This observation suggests fitting the function to the numerical data of (and similarly for the other derivatives), with and as fit parameters (see left and center plots of Fig. 2). Although the quality of this fit (having a residual sum of squares of ) is acceptable, it can be improved by about two orders of magnitude by using the fit function
(11) 
based on the sharper bound in (5c), with fit parameters and . This fit is of excellent quality, indicating that the distancedependence of the coefficients in (7) appreciably modifies the shape of the propagation front and correctly reproduces the actual spreading.
The sharper bounds in (5a)–(5d) inherit a systemsize dependence through the lattice dependence of and in the coefficients . As a result, (at fixed and fixed distance ) scales differently with for the cases , , and , respectively Note1 (). The switching from one regime to another at and nicely coincides with the different scaling regimes of equilibration times observed in Bachelard and Kastner (2013). Additionally to the dependence inherent to the bound, we find that, for , the optimal values for the fit parameters and in (11) show a strong dependence, wellcaptured by a power law for both parameters (Fig. 3 left and center). This scaling seems to originate from the dependence of the prefactor in (3) and (8), and is seen as an indication that the bounds could be further improved. For , in contrast, the dependence of and is negligible.
Comparison to the quantum mechanical bound.— LiebRobinson bounds were previously known for longrange interacting quantum systems Hastings and Koma (2006); Nachtergaele et al. (2006). The functional form of these bounds is , with a constant that depends on the observables considered. Asymptotically for long times and large distances , this functional form coincides with that of the weaker form of all four classical bounds (5a)–(5d). For small , however, the bounds differ, not only between the classical and the quantum case, but also between the four derivatives bounded in the classical case. The shorttime behavior is linear in for in (5c), quadratic in for and in (5a) and (5d), and cubic in for in (5b). The numerical results in Fig. 3 (right) confirm that the real shorttime dynamics of the model is correctly captured by these different functional forms of the bounds. The quantum mechanical bound, in contrast, increases linearly for short times , independently of the observables considered, although this may not reflect the actual behavior of expectation values in all cases.
Conclusions.— We have reported LiebRobinsontype inequalities, bounding the speed at which a perturbation can travel across the lattice, for a broad class of longrange interacting classical lattice models (including models on arbitrary graphs , with multidimensional singleparticle phase spaces , manyparticle interactions, and for rather general forms of the interaction potential). By a suitable rescaling, we extended the bounds to arbitrary nonnegative longrange exponents , deep into the regime of strong longrange interactions. The weaker bounds on the righthand side of (5a)–(5d) are direct analogs of the quantum mechanical version of LiebRobinson bounds for longrange interacting systems Hastings and Koma (2006); Nachtergaele et al. (2006). While our numerical results for chains confirm the validity of these bounds, they reveal that the shape of the propagation front is not correctly captured. Only the stronger versions of the bounds in (5a)–(5d), with an additional distancedependence introduced through the coefficients , reproduce the functional form of the propagation front. These findings are in contrast to the shortrange case, where already the weaker “conventional” form of the LiebRobinson bound yields the correct, coneshaped spatiotemporal behavior in agreement with the numerical results.
Since our results apply to arbitrary classical observables, potential applications cover a broad range of dynamical phenomena in longrange interacting classical lattice models, from heat conduction to information transmission, energy transfer, and the approach to equilibrium. In the latter context, different finitesize scaling properties of equilibration times had been observed in the regimes , , and , respectively Bachelard and Kastner (2013). These three regimes agree precisely with the different scaling regimes of the coefficients that enter and reflect in the bounds (5a)–(5d), providing a theoretical explanation of the numerical observations.
Sharper LiebRobinson bounds, similar in spirit to (5a)–(5d), can also be derived for quantum mechanical lattice models with longrange interactions and will be reported in a forthcoming paper.
R. B. acknowledges support by the Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), and the computational support of the Núcleo de Apoio a Óptica e Fotônica (NAPOFUSP). M. K. acknowledges financial support by the National Research Foundation of South Africa via the Incentive Funding and the Competitive Programme for Rated Researchers.
Supplemental Material
Appendix A A. Proof of Equations (5)–(7)
According to Hamilton’s equations, the equations of motion corresponding to the Hamiltonian (3) are
(A.1a)  
(A.1b) 
Integrating (A.1a) and (A.1b) over time and taking partial derivatives with respect to or , one obtains
(A.2a)  
(A.2b)  
(A.2c)  
(A.2d)  
These are the four derivatives occurring in (2), and we want to derive upper bounds on their absolute values. Here we show only the derivation of the bound on the absolute value of (A.2a). Bounds on the other derivatives can be obtained by the same strategy. For the sake of a compact notation, we introduce the definitions
(A.3a)  
(A.3b)  
(A.3c) 
We denote by the matrix with elements , and we define the vectors and . Inserting (A.2c) into (A.2a) and expressing the result in terms of the definitions (A.3a)–(A.3c) yields
(A.4) 
where denotes the th component of the vector . Integrating by parts we obtain
(A.5) 
fold iteration of this formula gives
(A.6) 
where is a vector of Kronecker deltas. To prove a bound on this series, and hence also its convergence, in the large limit, we proceed by constructing an upper bound on (the absolute value of the elements of) the matrix products in (A.6),
(A.7) 
with as defined in (6). For the proof of (A.7) we require that
(A.8) 
with as defined in (8). We postpone a detailed discussion of this condition to Sec. B.1, where we prove that (A.8) is satisfied for power law decaying longrange interactions with exponents . Provided (A.8) holds, we can prove (A.7) by mathematical induction in the number of matrix multiplications.
Induction basis: For we have
(A.9a)  
(A.9b) 
Induction hypothesis: Assume (A.7) holds for some .
Inductive step:
For and , the lefthand side of (A.7) can be bounded by
(A.10) 
For and , the matrix product is bounded by
(A.11)  
This completes the proof of (A.7) for all .
Making use of the timeindependent bound (A.7), the timedependence of the integrand in (A.6) becomes trivial. Hence the integration can be performed by elementary means, yielding
(A.12) 
which proves the stronger (middle) bound in (5c) for all with . (The case is not relevant here.) Here we have assumed convergence of the series in (A.6), and the following calculation of the weaker bound in (5c) will confirm that this assumption is indeed justified.
An upper bound on is obtained by taking the supremum over pairs of the recursion relation (6),
(A.13) 
which is solved by . Inserting the latter inequality into (A.10) yields
(A.14) 
Inserting this timeindependent bound into (A.6) and performing the integration we obtain
(A.15) 
for all with . The second term in the round brackets vanishes in the limit . This implies convergence of the series, and we obtain
(A.16) 
which proves the weaker (rightmost) bound in (5c).
Appendix B B. Bounds for general classical longrange systems
As a lattice, we consider a set of vertices and a set of edges connecting pairs of vertices. The graph distance (number of edges of the shortest path connecting the sites ) serves as a metric on the graph .
To each vertex we assign a dimensional manifold as the local configuration space. For any finite subset , the configuration space associated with is the product space , and the phase space associated with is the cotangent bundle . We define the Hamiltonian as a real function on phase space, and, for given initial conditions, it generates a flow on in the usual way via Hamilton’s equations. Within this setting, following Section 2 of Marchioro et al. (1978), one can show that Eqs. (1) and (2) of the main text hold for differentiable functions with bounded derivatives.
We consider Hamiltonian functions of a standard form,
(B.1) 
consisting of a kinetic term quadratic in the momenta , and a general interaction term depending on the coordinates via the differentiable vertex interactions , with the number of elements in the set . The normalization factor is chosen such that is extensive ^{3}^{3}3We allow both, and , to explicitly depend on . This can be necessary when, for example, 2body and 3body interactions are present on a dimensional regular lattice and both interactions decay like with the distance and some exponent smaller than . In this case the 2 and 3body terms require different normalization factors in order to guarantee extensivity of the Hamiltonian, and this difference can be accounted for by means of a dependent 3body interaction ..
For proving a LiebRobinsontype bound, we require the interactions to satisfy
(B.2a)  
(B.2b) 
For local configuration space manifolds of dimensions greater than 1, is a matrix. Here and for the remainder of this section, the vertical bars denote an elementwise matrix maximum norm, for a matrix with elements . Similar to Section 1.1 of Nachtergaele et al. (2006), the nonincreasing function has the following properties:

is uniformly summable over , i.e.,
(B.3) 
satisfies
(B.4)
Within this setting, we obtain the following bounds on the partial derivatives on the righthand side of (2),
(B.5a)  
(B.5b)  
(B.5c)  
(B.5d) 
with
(B.6) 
The coefficients are defined recursively by
(B.7) 
with and for ,
(B.8)  
(B.9) 
Proof.—The equations of motion corresponding to the Hamiltonian (B.1) are
(B.10a)  
(B.10b) 
Integrating over and take the partial derivative with respect to or , we obtain
(B.11a)  
(B.11b)  
(B.11c)  
(B.11d) 
where denotes the identity matrix and
(B.12) 
Inserting (B.11c) into (B.11a) and introducing the definition
(B.13) 
we obtain
(B.14) 
This equation is formally identical to (A.4), the only difference being that each entry of the matrix and of the vector is now a matrix. Integrating (B.14) by parts and fold iteration of the resulting expression yields