A A. Proof of Equations (5)–(7)

# Spreading of Perturbations in Long-Range Interacting Classical Lattice Models

## Abstract

Lieb-Robinson-type 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 Lieb-Robinson bounds, is required to correctly capture the shape of the causal region, as confirmed by numerical results for classical long-range chains. We discuss the relevance of our findings for the relaxation to equilibrium of long-range 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 Lieb-Robinson bounds Lieb and Robinson (1972); Nachtergaele and Sims (2010); ?, 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 equal-time 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 Lieb-Schultz-Mattis theorem in higher dimensions Hastings (2004). Experimental observations related to Lieb-Robinson bounds have also been reported Cheneau et al. (2012); ?.

The original proof by Lieb and Robinson Lieb and Robinson (1972) requires interactions of finite range. An extension to power-law-decaying long-range 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 long-range 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 Lieb-Robinson-type results are known. This fact nicely fits into the larger picture that, for , the behavior of long-range interacting systems often differs substantially from that of short-range interacting systems. Examples of such differences include nonequivalent equilibrium statistical ensembles and negative response functions Lynden-Bell and Wood (1968); ?; ?; ?; ?; ?, or the occurrence of quasistationary states whose lifetimes diverge with the system size Antoni and Ruffo (1995); ?; ?; 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 Lieb-Robinson bounds.

In most cases the peculiarities of long-range interacting systems have been investigated in the framework of classical Hamiltonian systems Campa et al. (2009), but little is known about Lieb-Robinson bounds in classical mechanics. Exceptions are restricted to specific models with nearest-neighbor interactions Marchioro et al. (1978); Buttà et al. (2007); ?. In the context of classical Hamiltonian mechanics, a Lieb-Robinson 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),

 |{fA(0),gB(t)}|≤|A||B|∥∇f∥∞∥∇g∥∞uAB(t) (1)

where and are the (bounded) maxima of all the partial derivatives of and with respect to the phase space coordinates, and

 Missing or unrecognized delimiter for \right (2)

The partial derivatives on the right-hand side of (2) quantify the effect that a variation of the initial momentum or position , at the lattice site have on the time-evolved momentum or position , at the lattice site . A classical Lieb-Robinson 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 long-range interacting lattice models. We have proved a Lieb-Robinson-type result, providing an upper bound on in (2) [and hence on the Poisson bracket (1)] for a broad class of classical long-range interacting lattice models in arbitrary spatial dimension. Dipolar interactions in condensed matter systems are the prime example of such long-range 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, 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 (1-norm) distance between lattice sites and . The value of the exponent determines the range of the interaction, from mean-field-type (distance-independent) interactions at to nearest-neighbor 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 finite-range 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 short-range case, such a bound can be used to rigorously control finite-size 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 long-range interacting chain. Besides confirming the validity of the bound, the numerical results reveal that the refined version of our bound, sharper than conventional Lieb-Robinson-type 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

 H=∑i∈Λp2i2−JΛ2∑i,j∈Λi≠jcos(qi−qj)|i−j|α. (3)

For , the second sum on the right-hand 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 Lieb-Robinson bound requires the Hamiltonian to be extensive. We enforce extensivity also for by allowing the coupling constant to depend explicitly on the lattice,

 JΛ=J/supi∈Λ∑j∈Λ∖{i}1|i−j|α, (4)

where is a real constant 2.

Classical long-range Lieb-Robinson bound.—Upper bounds on the partial derivatives on the right-hand side of (2) are given by

 ∣∣∣∂qj(t)∂qi(0)∣∣∣ ≤∞∑n=1Uijnt2n(2n)!|i−j|α≤cosh(vt)−1|i−j|α=:Bqqij(t), (5a) ∣∣∣∂qj(t)∂pi(0)∣∣∣ ≤∞∑n=1Uijn|t|2n+1(2n+1)!|i−j|α≤sinh|vt|−|vt|v|i−j|α=:Bqpij(t), (5b) ∣∣∣∂pj(t)∂qi(0)∣∣∣ ≤∞∑n=1Uijn|t|2n−1(2n−1)!|i−j|α≤vsinh|vt||i−j|α=:Bpqij(t), (5c) ∣∣∣∂pj(t)∂pi(0)∣∣∣ ≤∞∑n=1Uijnt2n(2n)!|i−j|α≤cosh(vt)−1|i−j|α=:Bppij(t). (5d)

The positive coefficients are defined recursively by

 Uijn+1={|J|Uijn+|JΛ|Uiin+CijUmaxnfor i≠j,|J|Uiin+CiiUmaxnfor i=j, (6)

with and , and we use the constants

 Umaxn=supi,j∈ΛUijn,v=supi,j∈Λ√|J+JΛ|+Cij, (7) Cij=⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩|JΛ|∑k∈Λ∖{i,j}|i−j|α|i−k|α|j−k|αfor i≠j,|JΛ|∑k∈Λ∖{i}1|i−k|2αfor i=j. (8)

The proof of the bounds combines techniques for classical lattices with nearest-neighbor interactions Marchioro et al. (1978) with those used for proving Lieb-Robinson bounds for long-range quantum systems Hastings and Koma (2006); Nachtergaele et al. (2006), with the additional refinement of allowing lattice-dependent 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 right-hand 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

 v|t|>ln(2ϵ)+αln|i−j|. (9)

This effective causal region, in which the effect of an initial perturbation is non-negligible, has a logarithmic shape, and differs in this respect from the linear (cone-shaped) region derived in Marchioro et al. (1978) for short-range interactions. As a consequence, no finite group velocity limits the spreading of perturbations in long-range interacting lattices, and supersonic propagation can occur. For a refined understanding of the spatiotemporal behavior, we go beyond the usual Lieb-Robinson-type 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 single-particle phase spaces , many-particle 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 one-dimensional chain of sites with periodic boundary conditions. The partial derivatives in (5a)–(5d) are approximated by difference quotients,

 ∂pj(t)∂qi(0)≈~pj(t)−pj(t)δqi=:Qpqij(t), (10)

and similarly for the other derivatives. Here, is the time-evolved 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 time-evolved momenta and are obtained by numerically integrating Hamilton’s equations using a sixth-order 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 distance-dependence is shown in Fig. 1. The plots illustrate the supersonic propagation of perturbations in the presence of long-range 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

 ~Bpqij(t)=~c|i−j|α∞∑n=1Uijn|t/~z|2n−1(2n−1)! (11)

based on the sharper bound in (5c), with fit parameters and . This fit is of excellent quality, indicating that the distance-dependence 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 system-size 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, well-captured 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.— Lieb-Robinson bounds were previously known for long-range 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 short-time 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 short-time 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 Lieb-Robinson-type inequalities, bounding the speed at which a perturbation can travel across the lattice, for a broad class of long-range interacting classical lattice models (including models on arbitrary graphs , with multidimensional single-particle phase spaces , many-particle interactions, and for rather general forms of the interaction potential). By a suitable rescaling, we extended the bounds to arbitrary non-negative long-range exponents , deep into the regime of strong long-range interactions. The weaker bounds on the right-hand side of (5a)–(5d) are direct analogs of the quantum mechanical version of Lieb-Robinson bounds for long-range 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 distance-dependence introduced through the coefficients , reproduce the functional form of the propagation front. These findings are in contrast to the short-range case, where already the weaker “conventional” form of the Lieb-Robinson bound yields the correct, cone-shaped 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 long-range interacting classical lattice models, from heat conduction to information transmission, energy transfer, and the approach to equilibrium. In the latter context, different finite-size 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 Lieb-Robinson bounds, similar in spirit to (5a)–(5d), can also be derived for quantum mechanical lattice models with long-range 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 (NAPOF-USP). 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

 ddtpj(t) =−JΛ∑i∈Λi≠jsin[qj(t)−qi(t)]|i−j|α, (A.1a) ddtqj(t) =pj(t). (A.1b)

Integrating (A.1a) and (A.1b) over time and taking partial derivatives with respect to or , one obtains

 ∂qj(t)∂qi(0) =δij+∫t0∂pj(s)∂qi(0)ds, (A.2a) ∂qj(t)∂pi(0) =∫t0∂pj(s)∂pi(0)ds, (A.2b) ∂pj(t)∂qi(0) =−∫t0JΛ (A.2c) × ∑k∈Λk≠j[cos[qj(s)−qk(s)]|j−k|α(∂qj(s)∂qi(0)−∂qk(s)∂qi(0))]ds, ∂pj(t)∂pi(0) =δij−∫t0JΛ (A.2d) × ∑k∈Λk≠j[cos[qj(s)−qk(s)]|j−k|α(∂qj(s)∂pi(0)−∂qk(s)∂pi(0))]ds.

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

 ψi =(ψij)j∈Λwithψij(t)=∂qj(t)∂qi(0), (A.3a) Ajj(t) =−JΛ∑k∈Λk≠j(cos[qj(t)−qk(t)]|j−k|α), (A.3b) Ajk(t) Missing or unrecognized delimiter for \right (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

 ψij(t)=δij+∫t0∫t10[Aψi]j(t2)dt2dt1, (A.4)

where denotes the th component of the vector . Integrating by parts we obtain

 ψij(t)=δij+∫t0(t−t1)[Aψi]j(t1)dt1. (A.5)

-fold iteration of this formula gives

 ψij(t)=δij+M∑m=1(∫t0∫t10⋯∫tm−10dtm⋯dt1×(t−t1)⋯(tm−1−tm)[A(t1)⋯A(tm)δi]j)+∫t0∫t10⋯∫tM0(t−t1)⋯(tM−tM+1)×[A(t1)⋯A(tM+1)ψi(tM+1)]jdtM+1⋯dt1, (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),

 ∣∣∣∑k1,…,kn−1∈ΛAjk1(t1)⋯Akn−1i(tn)∣∣∣≤⎧⎪⎨⎪⎩Uijn|i−j|α\!\!\!for i≠j,Uiin\!\!\!for i=j, (A.7)

with as defined in (6). For the proof of (A.7) we require that

 Cij<∞∀i,j∈Λ, i≠j (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 long-range 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

 |Aii| Missing or unrecognized delimiter for \right (A.9a) |Aij| ≤|JΛ||i−j|α=Uij1|i−j|α% for i≠j. (A.9b)

Induction hypothesis: Assume (A.7) holds for some .
Inductive step: For and , the left-hand side of (A.7) can be bounded by

 ∣∣∣∑k1,…,kn∈ΛAjk1(t1)⋯Akn−1kn(tn)Akni(tn+1)∣∣∣ Missing or unrecognized delimiter for \right ≤Umaxn∑kn∈Λkn≠i,j1|j−kn|α|JΛ||kn−i|α Missing or unrecognized delimiter for \right CijUmaxn+|J|Ujin+|JΛ|Ujjn|i−j|α=Ujin+1|i−j|α. (A.10)

For and , the matrix product is bounded by

 ∣∣∣∑k1,…,kn∈ΛAik1(t1)⋯Akn−1kn(tn)Akni(tn+1)∣∣∣ (A.11) Umaxn∑kn∈Λkn≠i1|i−kn|α|Akni(tn+1)|+Uiin|Aii(tn+1)| ≤ CiiUmaxn+|J|Uiin=Uiin+1.

This completes the proof of (A.7) for all .

Making use of the time-independent bound (A.7), the time-dependence of the integrand in (A.6) becomes trivial. Hence the integration can be performed by elementary means, yielding

 Missing or unrecognized delimiter for \right (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),

 Umaxn+1 ≤supi,j∈Λ(|J|Uijn+|JΛ|Uiin+CijUmaxn) ≤Umaxn(supi,j∈ΛCij+|J+JΛ|)=Umaxnv2 (A.13)

which is solved by . Inserting the latter inequality into (A.10) yields

 ∣∣∣∑k1,…,kn−1,kn∈ΛAjk1(t1)⋯Akn−1kn(tn)Akni(tn+1)∣∣∣≤v2n+2. (A.14)

Inserting this time-independent bound into (A.6) and performing the integration we obtain

 Missing or unrecognized delimiter for \right (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

 Missing or unrecognized delimiter for \left (A.16)

which proves the weaker (rightmost) bound in (5c).

## Appendix B B. Bounds for general classical long-range 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,

 HΛ=∑i∈Λp2i2mi+NΛ∑X⊂ΛΦXΛ({qk}k∈X), (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.

For proving a Lieb-Robinson-type bound, we require the interactions to satisfy

 supi,j∈Λ∑X∋i,jsupq∈\mathbbmYX∣∣ ∣∣∂2ΦXΛ∂qi∂qj({qk}k∈X)∣∣ ∣∣F(d(i,j)) <∞for i≠j, (B.2a) NΛsupi∈Λ∑X∋isupq∈\mathbbmYX∣∣ ∣∣∂2ΦXΛ∂q2i({qk}k∈X)∣∣ ∣∣F(d(i,i)) <∞. (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 non-increasing function has the following properties:

1. is uniformly summable over , i.e.,

 ∥FΛ∥:=supi∈Λ∑j∈ΛNΛF(d(i,j))<∞, (B.3)
2. satisfies

 CΛ=NΛsupi,j∈Λ∑k∈Λk≠iF(d(i,k))F(d(k,j))F(d(i,j))<∞. (B.4)

Within this setting, we obtain the following bounds on the partial derivatives on the right-hand side of (2),

 ∣∣∣∂qj(t)∂qi(0)∣∣∣ ≤F(d(i,j))μmj∞∑n=1μnVijnt2n(2n)! ≤F(d(i,j))μmj(cosh(vt)−1), (B.5a) ∣∣∣∂qj(t)∂pi(0)∣∣∣ ≤F(d(i,j))μmj∞∑n=1μnVijnt2n+1(2n+1)! ≤F(d(i,j))μvmj(sinh(vt)−vt), (B.5b) ∣∣∣∂pj(t)∂qi(0)∣∣∣ ≤F(d(i,j))μmj∞∑n=1μnVijnt2n(2n)! ≤vF(d(i,j))μmj(sinh(vt)), (B.5c) ∣∣∣∂pj(t)∂pi(0)∣∣∣ ≤F(d(i,j))μmj∞∑n=1μnVijnt2n(2n)! ≤F(d(i,j))μmj(cosh(vt)−1), (B.5d)

with

 v=√μ[CΛ+F(0)]. (B.6)

The coefficients are defined recursively by

 Vijn+1=CijVmaxn+F(0)Vijn (B.7)

with and for ,

 Vmaxn =supi,j∈ΛVijn, (B.8) Cij =NΛ∑k∈Λk≠iF(d(i,k))F(d(k,j))F(d(i,j)). (B.9)

Proof.—The equations of motion corresponding to the Hamiltonian (B.1) are

 ddtpj(t) =−NΛ∑X∋j∂ΦXΛ∂qj({qk(t)}k∈X), (B.10a) ddtqj(t) =pj(t)mj. (B.10b)

Integrating over and take the partial derivative with respect to or , we obtain

 ∂qj(t)∂qi(0)= δij\mathbbm1μ+1mj∫t0∂pj(s)∂qi(0)ds, (B.11a) ∂qj(t)∂pi(0)= 1mj∫t0∂pj(s)∂pi(0)ds, (B.11b) ∂pj(t)∂qi(0)= −NΛ∫t0∑X∋j∂∂qi(0)∂ΦXΛ({qk(s)}k∈X)∂qjds = −NΛ∫t0∑X∋j∑k∈X∂2ΦXΛ∂qj∂qk∂qk(s)∂qi(0)ds = ∫t0∑k∈ΛAjk∂qk(s)∂qi(0)ds, (B.11c) ∂pj(t)∂pi(0)= δij\mathbbm1μ+∫t0∑k∈ΛAjk∂qk(s)∂pi(0)ds, (B.11d)

where denotes the identity matrix and

 Ajk=−NΛ∑X∋j,k∂2ΦXΛ∂qj∂qk. (B.12)

Inserting (B.11c) into (B.11a) and introducing the definition

 ψi=(ψij)j∈Λwithψij(t)=∂qj(t)∂qi(0) (B.13)

we obtain

 ψij(t)=δij\mathbbm1μ+∫t0∫t10[Aψi]j(t2)dt2dt1. (B.14)

This equation is formally identical to (A.4), the only difference being that each entry of the