Diffusive Propagation of Energy in a Non-Acoustic Chain

We consider a non acoustic chain of harmonic oscillators with the dynamics perturbed by a random local exchange of momentum, such that energy and momentum are conserved. The macroscopic limits of the energy density, momentum and the curvature (or bending) of the chain satisfy a system of evolution equations}. We prove that, in a diffusive space-time scaling, the curvature and momentum evolve following a linear system that corresponds to a damped Euler-Bernoulli beam equation. The macroscopic energy density evolves following a non linear diffusive equation. In particular the energy transfer is diffusive in this dynamics. This provides a first rigorous example of a normal diffusion of energy in a one dimensional dynamics that conserves the momentum.


Introduction
Macroscopic transport in a low dimensional system, in particular the energy transport, has attracted attention in both the physics and mathematical physics literature in the latest decades. Anomalous energy transport has been observed numerically in Fermi-Pasta-Ulam (FPU) chains, with the diverging thermal conductivity [8]. Generically this anomalous superdiffusive behavior is attributed to the momentum conservation properties of the dynamics [7]. Actually one dimensional FPU-type chains have potential energy depending on the interparticle distances (i.e. the gradients of the particles displacements), and have three main locally conserved quantities: volume stretch, momentum and energy. These conserved (or balanced) quantities have different macroscopic space-time scalings, corresponding to different type of initial non-equilibrium behaviour. A mechanical non-equilibrium initial profile due to the gradients of the tension induces a macroscopic ballistic evolution, at the hyperbolic space-time scale, governed by the Euler equations (cf. [3]). When the system approaches to, or is already at a mechanical equilibrium, the temperature profile will evolve at a superdiffusive time scale.
Recent heuristic calculations based on fluctuating hydrodynamics theory [9], connect the macroscopic space-time scale of the superdiffusion of the thermal (energy) mode to the diffusive or superdiffusive fluctuations of the other conserved quantities. It turns out that this  superdiffusive behavior of the energy is governed by a fractional laplacian heat equation. This picture can be mathematically rigorously proven in the case of a harmonic chain perturbed by a local random exchange of momentum, see [6,5]. In particular, it has been shown in [6], that in the models driven by the tension, there is a separation of the time evolution scales between the long modes (that evolve on a hyperbolic time scale) and the thermal short modes that evolve in a longer superdiffusive scale. In addition, from the explicit form of the macroscopic evolution appearing in these models, it is clear that this behavior is strongly dependent on a non-vanishing speed of sound. More specifically, when the speed of sound is null, there is no macroscopic evolution either at the hyperbolic or superdiffusive time scales. This suggests that the macroscopic evolution of the system should happen at a yet longer, possibly diffusive, time scale for all modes.
In the present article we investigate the harmonic chain model with the random exchange of momenta. The interaction potential depends only on the squares of the curvature (or bending) of the chain where q x are the positions of the particles. This means that its hamiltonian is formally given by H(k, p) = x e x (k, p), (1.2) where the energy of the oscillator x is defined e x (k, p) := p 2 Here α is a positive parameter that indicates the strength of the springs. This corresponds to a special choice of attractive nearest neighbor springs and repulsive next nearest neighbor springs. It turns out that the respective speed of sound is null, even though the momentum is conserved by the dynamics. As the energy depends on the curvature and not on the volume, this system is tensionless, and the corresponding relevant conserved quantity, besides the energy and momentum, is the curvature and not the volume stretch. Our first result, see Theorem 3.1 below. asserts that these three conserved quantities (curvature, momentum and energy) evolve together in the diffusive time scale. Curvature and momentum are governed macroscopically by the damped Euler-Bernoulli beam equations: ∂ t k(t, y) = −∆ y p(t, y), ∂ t p(t, y) = α∆ y k(t, y) + 3γp(t, y) , (1.4) where γ > 0 is the intensity of the random exchange of momentum.
Defining the mechanical macroscopic energy as e mech (t, y) = 1 2 p 2 (t, y) + αk 2 (t, y) (1.5) and its thermal counterpart (or temperature profile) as e th (t, y) = e(t, y) − e mech (t, y) (1.6) the evolution of the latter is given by ∂ t e th (t, y) = ( √ 3 − 1)α 2 √ 3γ + 3γ ∆ y e th (t, y) + 3γ (∂ y p(t, y)) 2 , (1.7) see Theorem 3.2. In particular, the thermal conductivity is finite and we have a normal diffusion in this system. Notice also that because of the viscosity term, a gradient of the macroscopic velocity profile induces a local increase of the temperature. This result puts in evidence two main differences between the present and the FPU-type models: (i) the thermal conductivity is finite, even though the system is one dimensional and dynamics conserves the momentum. This suggests that the non-vanishing speed of sound is a necessary condition for the superdiffusion of the thermal energy, (ii) there is no separation of the time scales between low (mechanical) and high (thermal) energy modes: all the frequencies evolve macroscopically in the diffusive time scale. Furthermore there is a continuous transfer of energy from low modes to high modes, resulting in the rise of the temperature, due to the gradients of the momentum profile. These rigorous results on the harmonic non-acoustic chain lead us to conjecture that a similar behavior is expected for the deterministic non-linear hamiltonian dynamics corresponding to an interaction of the type V (k x ), i.e. the energy is a non linear function of the curvature of the chain.
About our proof of the hydrodynamic limit: this is a non-gradient dynamics (microscopic energy currents are not of the form of discrete space gradients of some functions). Therefore, we cannot use known techniques for such type of limits based on relative entropy methods (cf. e.g. [10], [11]) for two reasons: (i) lack of control of higher moments of the currents in terms of the relative entropy, (ii) degeneracy of the noise in the dynamics, as it acts only on the velocities.
Instead, we develop a method already used in [5], based on Wigner distributions for the energy of the acoustic chain. Thanks to the energy conservation property of the dynamics we can easily conclude, see Section 5.4, that the Wigner distributions form a compact family of elements in a weak topology of an appropriate Banach space. Our main result concerning the identification of its limit is contained in Theorem 5.1 below. The spatial energy density is a marginal of the Wigner function. We would like to highlight the fact, that in addition to proving the hydrodynamic limit of the energy functional, we are also able to identify the distribution of the macroscopic energy in the frequency mode domain, see formula (5.30). In particular the thermal energy is uniformly distributed on all modes (which is a form of local equilibrium), while the macroscopic mechanical energy is concentrated on the macroscopic low modes, see (5.30).
To show Theorem 5.1 we investigate the limit of the Laplace transforms of the Wigner distributions introduced in Section 7. The main results, dealing with the asymptotics of the Laplace-Wigner distributions, are formulated in Theorems 7.1 -7.3. Having these results we are able to finish the identification of the limit of the Wigner distributions, thus ending the proof of Theorem 5.1. The proofs of the aforementioned Theorems 7.1 -7.3, which are rather technical, are presented in Sections 8 -10, respectively.

The dynamics
2.1. Non acoustic chain of harmonic oscillators. Since in the nonacoustic chain the potential energy depends only on the bendings, see (1.1), in order to describe the configuration of the infinite chain we only need to specify (k x ) x∈Z , and the configurations of our dynamics will be denoted by ((p x , k x )) x∈Z ∈ (R × R) Z .
In case when no noise is present the dynamics of the chain of oscillators can be written formally as a Hamiltonian system of differential equationsk 2.1.1. Continuous time noise. We add to the right hand side of (2.1) a local stochastic term that conserves both p 2 x−1 +p 2 x +p 2 x+1 and p x−1 +p x + p x+1 . The respective stochastic differential equations can be written as with the parameter γ > 0 that indicates the strength of the noise in the system, and (Y x ) are vector fields given by Here (w x (t)) t≥0 , x ∈ Z are i.i.d. one dimensional, real valued, standard Brownian motions, over a probability space (Ω, F , P). Furthermore, As a result we obtain We can rewrite the system (2.2)

TOMASZ KOMOROWSKI AND STEFANO OLLA
Remark. The particular choice of the random exchange in the above dynamics is not important. The result can be extended to any other random mechanism of moment exchange, as long as total energy and momentum are conserved. Most simple dynamics would be given by exchange of momentum between nearest neighbor atoms at independent exponential times.
are stationary for the dynamics defined by (2.4). In this context τ is called the load of the chain, while as usual β −1 is the temperature and p is the average momentum. Notice that, when τ = 0, the above distribution is spatially translation invariant, only for the (k x , p x ) coordinates, but is is not translation invariant with respect to the position q x , or the stretch r x = q x − q x−1 coordinates.
2.3. Initial data. Concerning the initial data we assume that, given ǫ > 0, it is distributed according to a probability measure µ ǫ on the configuration of ((k x , p x )) x∈Z and satisfies sup ǫ∈(0,1] ǫ x e x µǫ < +∞. (2.6) Here · µǫ denotes the average with respect to µ ǫ . We denote also by E ǫ the expectation with respect to the product measure P ǫ = µ ǫ ⊗ P. The existence and uniqueness of a solution to (2.2) in ℓ 2 , with the aforementioned initial condition can be easily concluded from the standard Hilbert space theory of stochastic differential equations, see e.g. Chapter 6 of [2].
We assume furthermore that the mean of the initial configuration varies on the macroscopic spatial scale: for some functions κ, p ∈ C ∞ 0 (R). Their Fourier transformsκ andp belong to the Schwartz class S(R). As for the fluctuations around the mean we assume that their energy spectrum is uniformly L r integrable with respect to ǫ > 0 for some r > 1. We have denoted bŷ (2.8) the Fourier trasform of a given sequence f x , x ∈ Z. Here T is the unit torus, understood as the interval [−1/2, 1/2] with the identified endpoints. Let The energy spectrum is defined as: wherep(k) andk(k) are the Fourier transforms of (p x ) and (k x ), respectively. Assumption (2.6) implies in particular that The announced property of the L r integrability of the energy spectrum means that there exists r > 1 such that: Thanks to the hypothesis (2.7) we conclude that for any G ∈ C ∞ 0 (R) we have The quantities p(·), κ(·) are called the macroscopic velocity and curvature profiles. We assume furthermore that the following limits exist for any G ∈ C ∞ 0 (R). Here j(·), p 2 (·), κ 2 (·) are some functions belonging to C ∞ 0 (R).

TOMASZ KOMOROWSKI AND STEFANO OLLA
As a consequence, we conclude that the limit also exists for any G ∈ C ∞ 0 (R). Here e(y) -the macroscopic energy profile -is given by Remark. An important example of initial distributions that satisfy the above conditions is provided by local Gibbs measures, i.e. inhomogeneous product probability measures of the type Here the vector λ x = (β x , p x , τ x ) is given by λ x := λ(ǫx), where λ(x) := (β(x), p(x), τ (x)) and the functions β −1 (·), p(·), τ (·) belong to C ∞ 0 (R). The deterministic field G(λ x ), x ∈ Z, called the Gibbs potential is given by an analogue of the second equality of (2.5).
Our first result concerns the evolution of the macroscopic profiles of the velocity and curvature.
The proof of this result is fairly standard and we show it in Section 11. Define the macroscopic profile of the mechanical energy of the chain by e mech (t, y) := 1 2 p 2 (t, y) + ακ 2 (t, y) , with e mech (y) := e mech (0, y). Comparing with the energy profile at t = 0, given by (2.16), we conclude that the residual energy, called the initial thermal energy (or temperature) profile, satisfies e th (y) := e(y) − e mech (y) ≥ 0. (3.4) Concerning the evolution of the energy profile we have the following result.
(3.5) exists. In addition, we have e(t, y) = e th (t, y) + e mech (t, y) where e mech (t, y) is given by (3.3), while the thermal energy (temperature) e th (t, y) is the solution of the following Cauchy problem: ∂ t e th (t, y) =ĉ∂ 2 y e th (t, y) + 3γ(∂ y p) 2 (t, y), e th (0, y) = e th (y). (3.6) The diffusivity coefficient equalŝ Remark 3.3. Notice that the gradient of the macroscopic momentum p(t, y) appearing in (3.6) causes a local increase of the temperature. It is also straightforward to understand the appearance of this term in the aforementioned equation. Consider for simplicity the case α = 0. The dynamics is constituted then only by the random exchanges of the momentum. The conserved quantities that evolve macroscopically are the momentum p x and the kinetic energy p 2 x /2. The corresponding macroscopic equations are These can be proven easily since the microscopic dynamics is of gradient type. It follows that the macroscopic equation for the temperature field, defined by e th = e − p 2 /2, is given by The interaction α affects the thermal diffusivity, but does not influence the nonlinearity appearing in the evolution of the temperature profile.

Some basic notation
To abbreviate our notation we write Given a set A and two functions f, g : We write g(x) f (x), when only the upper bound on g is satisfied.
Denote by S the set of functions J : R × T → C that are of C ∞ class and such that for any integers l, m, n we have sup y∈R, k∈T For J ∈ S we letĴ be its Fourier transform in the first variable, i.e.
For any M > 0 let A M be the completion of S in the norm Here B M := [η : |η| < M]. We drop the subsrcipt from the notation if M = +∞. Let A ′ and A ′ M be the respective topological dual spaces of A and A M .

5.
Wigner function and its evolution 5.1. The wave function. The wave function corresponding to the configuration ((p x , k x )) x∈Z is defined as Its Fourier transform is given bŷ The energy and its spectrum (2.10) can be written as Using the decomposition into the macroscopic profile and the fluctuation part, see (2.9), we can write are the wave functions corresponding to the macroscopic profile and the fluctuation part, respectively.

Wigner functions.
By the Wigner functions corresponding to the wave function field (ψ x ) x∈Z we understand four tempered distributions W ǫ,± , Y ǫ,± that we often write together in the form of a vector and for any J ∈ A. Here W ǫ,± (η, k) and Y ǫ,± (η, k) -called the Fourier-Wigner functions -are given by Using the Cauchy-Schwartz inequality and (2.11) we get A simple calculation shows that for any function J(y, k) ≡ J(y) where Using the decomposition of the wave function into its mean, following a macroscopic profile φ(·), and the fluctuation part {ψ (ǫ) x , x ∈ Z}, see (5.4), we can correspondingly decompose the vector of the Wigner functions. Namely 12) where the Fourier-Wigner function corresponding to these wave functions shall be denoted by where, using the Poisson summation formula, we have defined The formulas for Y ǫ,± , J and Y ǫ,± , J are constructed analogously using the respective Fourier-Wigner functions. Notice that for small ǫ the expression above of W ǫ,± is well approximated by the more natural definition: As a consequence of assumption (2.14) we conclude that for functions J(y, k) = J(y): with j(·), κ 2 (·) and p 2 (·) given by (2.14). A simple calculation also shows that Thus, One can also easily check that

Evolution of the wave function.
Adjusted to the macroscopic time, we can define the wave function corresponding to the configuration at time t/ǫ 2 is the unique solution of the Itô stochastic differential equation, understood in the mild sense (see e.g. Theorem 7.4 of [2]) is a dispersion relation, A simple calculation shows thatβ(k) = 4R(k). The process B(dt, dk) is a space-time Gaussian white noise, i.e.
To characterize the limit we recall that the thermal energy density e th (t, y) is given by the solution of the Cauchy problem (3.6), while the mechanical one e mech (t, y) is defined by (3.3). The limit of the Wigner functions corresponding to the macroscopic profile wave function Our main result concerning the limit of the Wigner transform can be stated as follows.
Theorem 5.1. Suppose that the initial data satisfy the assumptions Then, (W ǫ (t)) t≥0 converge, as ǫ → 0+, * -weakly over where W (t) is a measure on R × T given by Analogously to formulas (5.10) and (5.11) we can write (5.32) Therefore, the conclusion of Theorem 3.2 is a direct consequence of Theorem 5.1.
In addition, the remainder terms r can be rewritten in the matrix form where andD ǫ is a 4 × 4 matrix that can be written in the block form where A ǫ , B ǫ , C ǫ are 2 × 2 matrices given by and B ǫ =D (ǫ) − I 2 , with I n denoting the n × n identity matrix. We have also denoted ǫ (λ, η) be the column vectors obtained by scalar multiplication of each component of w ǫ (λ, η, k) by e ι . Note that and W ǫ (η, k) is the column vector corresponding to the Fourier-Wigner transforms of the components of (7.3), and ǫ,− , r ǫ,− , (6.25) 6.2. Invertibility of matrixD ǫ . We prove that the matrixD ǫ appearing in (6.20) is invertible, thus the vector of the Laplace-Fourier transforms of Wigner functions is uniquely determined by the system. It turns out to be true, provided that λ is sufficiently large.

Hered
For the abbreviation sake we denote by d j,ǫ , j = 1, . . . , 4 the vectors corresponding to the rows of the adjugate ofD ǫ given by (6.38). Combining the above with (6.13) and (6.26) we get.

Proof of Theorem 5.1
As we have already mentioned for any sequence ǫ n → 0+ there exists a subsequence W ǫ n ′ (t) that convergences * −weakly to some W ∈ (L 1 ([0, +∞); A)) * . We prove that the element W does not depend on the choice of the sequence ǫ n ′ by showing that for any M > 0 there exists λ 0 > 0 such that the vector w ǫ n ′ (λ) made of Laplace transforms of the components of W ǫ n ′ (t) converges * -weakly over A ′ M for any λ > λ 0 . In fact one can describe the respective limit as the Laplace transform of the vector W(t) appearing in the statement of Theorem 5.1. This identifies the limit of (W ǫ (t)), as ǫ → 0+ finishing in this way the proof of Theorem 5.1.

From (6.20) we obtain
Unfortunately, the right hand side of the above system contains also terms that depend on the vector w ǫ , via the projections of its components onto the vectors e ± and f ± . To describe the behavior of w ǫ we need to determine first these projections. Using (6.20) the above system can be rewritten in the form

3)
h ǫ is given by (6.25) and the 4 × 4 matrixẼ ǫ (λ, η, k) equals Multiplying both sides of (7.2) by e ι , ι ∈ {−, +} and then integrating over T we get a system of 8 equations Here w (ι) ǫ are column vectors obtained by a scalar multiplication of the entries of w ǫ (see (6.19)) by e ι . The same concerns Matrix G ǫ (λ, η) is a 2 × 2 block matrix of the form ± are 4 × 4 matrices defined as follows: Note that vector v ǫ appearing on the right hand side of (7.6) still depends on the projections of w ǫ onto f ± , cf (6.25) and (6.16). It turns out however that the asymptotics of these projections, as ǫ → 0+, can be described by only one of them, e.g. w ǫ . This is a conclusion of our next result. Denote by δw ǫ := w ǫ . We shall also use the following convention: for a given M > 0 the constants ǫ 0 , λ 0 > 0 are selected as in the statement of Proposition 6.1 so thatδ ǫ (λ, η, k) ≈ δ The proof of the theorem is presented in Section 8.
To describe the limit of w (−) ǫ (λ, η) we can use the the system (7.6), which is "almost closed" with respect to the components of w (−) ǫ , i.e. it is closed modulo some corrections that in light of Theorem 7.1 are of lower order of magnitude.
Let us first introduce some additional notation. Given the wave function φ(t, y) we define the vector of the Laplace-Fourier transforms of the respective macroscopic Wigner functions w T φ (λ, η, h) = [w φ,+ , y φ,+ , y φ,− , w φ,− ], (7.12) where Here Hereê th (η) is the Fourier transform of e th (y) appearing in (3.4) andĉ is given by (3.7). We can show, see Section 9 below for the proof, the following result.
To obtain the asymptotics of w ǫ (λ, η, k) we use (7.1), which allows us to describe the Fourier-Laplace transforms of the Wigner functions in terms of their projections onto e ± and f ± . We obtain then the following result.
The proof of this result is contained in Section 10.
The end of the proof of Theorem 5.1. Thanks to (5.24) we know that W ǫ (t) is sequentially pre-compact , as ǫ → 0+, in the * -weak topology of (L 1 ([0, +∞), A)) * . To identify its limiting points we consider w ǫ (λ, η, k) the vector of the Laplace-Fourier transforms of W ǫ (t). Given λ > 0 this family is sequentially pre-compact in the * -weak topology opf A ′ , as ǫ → 0+. Thanks to Theorems 7.2 and 7.3 we conclude that given M > 0 one can choose λ 0 as in the statement of Proposition 6.1, such that the the components of w ǫ (λ, η, k) converge * -weakly over A ′ M to the Laplace-Fourier transforms of the respective functions appearing in the claim of Theorem 5.1 for any λ > λ 0 . To finish the proof we only need to verify that w(λ, dη, dk) -the limit of w ǫ,+ (λ, η, k) (the limit of w ǫ,− can then be trivially concluded) agrees for λ > λ 0 with the Laplace transform of W (t, dy, dk) appearing in (5.30).
We can see therefore that e ′ th (t, y) satisfies (3.6). Thus the conclusion of Theorem 5.1 follows.

Proof of Theorem 7.1
We start with the following result.
Therefore, cf (6.27), we conclude that Thus, and the first equality of (8.7) follows.
Combining this with the second formula of (8.18) we obtain, by the Lebesgue dominated convergence theorem, Using the above together with bound (7.9) we conclude that expression (8.16) can be written as Then bound (7.8) would follow, provided we can show that the right hand side of the first equation of the system (7.6), given by z , is of order of magnitude O(1), as ǫ ≪ 1. To see that we write where the terms U ǫ,i , i = 1, 2, 3 are given by The fact that z (1,−) ǫ = O(1), as ǫ ≪ 1, can be argued in a similar way as it has been done in the case of z (2,−) ǫ , see (8.14) and (8.15) above. 8.3. Proof of (7.10). From (6.20) we obtain  Using the above and the already proved estimates (7.8), (7.9) and Lemma 8.1 we obtain Here f(k) ≡ 1.

8.3.2.
Convergence of II ǫ , III ǫ and IV ǫ . Thanks to Lemma 8.1 we can write The remaining terms appearing in expressions II ǫ , III ǫ and IV ǫ can be estimated in the same manner allowing us to conclude that 8.4. Proof of (7.11). From (6.20) we obtain The analysis of the above terms is very similar to what has been done in the precious section. Using (8.10) we conclude that for any λ > λ 0 We conclude in this way that all RI ǫ , RII ǫ , RIII ǫ and RIV ǫ tend to 0 in the L 1 sense. Thus, (7.11) follows. 9. Proof of Theorem 7.2 Since functions e ± (k) are both even the fourth and eighth equation of the system (7.6) coincide with the first and the fifth ones respectively.
Adding the first and fifth equations of the system (7.6) we get Here a   .7)). The second and third equations of (7.6) read (cf (8.9)) Adding sideways these equations we get Here z The sixth and seventh equations of (7.6) yield and Summarizing, we have obtained the following system Using Theorem 7.1 we conclude that given M > 0 and λ > λ 0 the family (w Subtracting sideways from the first equation of (9.8) the sum of the remaining two and taking into account (9.10) and (9.11) we obtain Moreover, a direct calculation shows that Using Lemma 8.1 and the Lebesgue dominated convergence theorem we obtain Using the above formula and substituting (cf (7.3) and (7.7)) we can rewrite (9.12) in the form Here is the solution of the system where W ǫ (η, k) is the column vector of Fourier-Wigner functions corresponding to the initial data, see (5.5). In addition, f ǫ is given by (6.16) respectively, and r ǫ,1 := Here are the solutions of the analogues of (9.19) in which the right hand side has been replaced by W ǫ and W ǫ , respectively.
Proof. We only prove (9.39), as the argument for (9.40) is very similar. The left hand side of (9.39) for w (0) ǫ,+ can be rewritten in the form Denote by J j,ǫ , j = 1, 2, 3, 4 the respective terms arising after opening of the square bracket. Changing variables k := k/ǫ we can write (cf (5.13))  (1))dk.
Using a similar argument we infer that for any λ > λ 0 , |η| ≤ M we have Thanks to (10.11) to compute the last limit we can use the Lebesgue dominated convergence and conclude, using (9.34) and (9.38), that the right hand side of the above equality coincides with the right hand side of (10.8).