Nicolas Toussaint
email: nicolas.toussaint@sophia.inria.fr

Jean-Christophe Souplet
email: jean-christophe.souplet@sophia.inria.fr

Pierre Fillard
email: pierre.fillard@sophia.inria.fr

MedINRIA : Medical Image Navigation and Research Tool by INRIA

Processing and visualization of 3D medical data is nowadays a common problem. However, it remains challenging because the diversification and complexification of the available sources of information, as well as the specific requirements of clinicians, make it difficult to solve in a computer science point of view. Indeed, clinicians need ergonomic, efficient, intuitive and reactive softwares. Moreover, they need new solutions to fully exploit their data, but they often cannot access state-ofthe-art methods as those are mostly available in complicated softwares. The MedINRIA software was born to fill this lack and consists of a collection of tools that optimally exploit various types of data (e.g., 3D images, diffusion tensor fields, neural fibers as obtained in DT-MRI). It provides state-of-the-art algorithms while keeping a user-friendly graphical interface. For each of these tools, we first introduce its dedicated application and the processing methods it contains. Then, we focus on the features that make interactions with data even more intuitive. Med-INRIA is a free software, available on Windows, Linux and MacOSX. Other MedINRIA tools are underway to make cutting edge research in medical imaging rapidly available to clinicians. The interest clinicians have shown in MedINRIA so far indicates that the need of such simple, yet powerful softwares is real and increasing.

Introduction

The increase and diversification of sources of information in medical imaging has raised the need for adapted tools for visualization and analysis of these data. For instance, the progress made in magnetic resonance imaging (MRI) for the past two decades has prodigiously improved the quality and resolution of 3D images. Moreover, these images can take various forms: anatomical MRI, functional MRI (fMRI) and diffusion tensor MRI (DT-MRI or DTI) [START_REF] Basser | MR diffusion tensor spectroscopy and imaging[END_REF] are only a few examples. The analysis of these data is challenging too. For instance, computing neural fibers in DTI or segmenting lesions in multiple sclerosis (MS) require adapted and sophisticated methods.

Most of state-of-the-art algorithms are not easily accessible to clinicians. These methods are, however, mostly available in relatively complex softwares, Fig. 1. The MedINRIA architecture. MedINRIA is composed of three main blocks: 1. The processing block based upon ITK, 2. The visualization block using VTK and 3. The user interface block made in wxWidgets. The blue arrows indicate the data flow, from the processing block to the user interface. The user-interface can only send controls to the processing and visualization blocks (the opposite is forbidden). Thus, the interface can be re-written using another library with minimum effort.

like Slicer [START_REF] Pieper | 3d slicer[END_REF] or BrainVisa [START_REF] Cointepas | Brainvisa: Software platform for visualization and analysis of multi-modality brain data[END_REF]. Although these softwares are extremely powerful, they can be quite hard to fully master. Clinicians needs are simple yet very difficult to satisfy in a computer science point of view: they want ergonomic, reactive and intuitive softwares that offer real-time interactions with data. The MedINRIA project was born to fulfill this need.

MedINRIA is a collection of softwares targeting the clinicians. It has the native look&feel on any system (Windows, MacOSX and Linux), and effort has been made to simplify the user interface while keeping state-of-the-art algorithms. Up to now, MedINRIA contains softwares for image visualization, DT-MRI analysis, tensor visualization, and semi-automatic segmentation of MS lesions. The paper is organized as follows: In Sec. 2, some basics and concepts of MedINRIA are introduced, and a description of each application is given, from their general purpose to the analysis and interactions with medical data they introduce. We conclude in Sec. 3 with a discussion on the future work.

The MedINRIA Software

Schematically, MedINRIA is composed of three main blocks, as shown in Fig. 1. First, the processing block is based upon ITK [START_REF] Ibanez | The ITK Software Guide[END_REF]. It has all ITK features: flexibility (it can process any type of data), multithreaded and a pipelined architecture. Second, the visualization block made with VTK [START_REF]The Visualization ToolKit[END_REF], itself based on OpenGL, offers the same pipelined structure as ITK. Finally comes the user interface block written in wxWidgets [6]. This block is responsible for transmitting user-defined parameters to the two first blocks. Unlike the processing and visualization blocks, this one can only transmit data in one direction, and consequently cannot receive any information from the other blocks. The main reason is to keep the user-interface detached from the processing and visualization part. Thus, one could easily re-write the user-interface with another library, or even write command line tools instead of graphical applications.

When started, all MedINRIA's tools can be launched from a single menu. Up to now, five applications are available, each of them being dedicated to a specific task or type of data. We give below a description of each of them: we start by briefly introducing the general purpose, and continue by focusing more on the user-interactions they propose to facilitate the intuitive exploration of data.

ImageViewer: A Simple yet Powerful Image Viewer

This first application aims at displaying 3D images and interacting with them. It is also used to convert clinical data (DICOM) into volumetric images. The ImageViewer introduces interesting new features like: image tabbed browsing, a preview screen (Fig. 2 left) where all images can be compared (interactions are synchronized: slice or contrast change of one window affects all the other windows as well), which gives an ergonomic navigation through images. In the following, we will only focus on the 3D interactions with images, as this is the most innovative part w.r.t. the simplification of data manipulation.

3D Image Visualization: Images are displayed in 3D either using multiplanar reconstruction (MPR) (Fig. 2 left) or volume rendering (VR) (Fig. 2 middle). MPR consists in displaying a slice in every orthogonal direction, while VR uses 3D textures with an opacity transfer function and is hardware-accelerated. The opacity transfer function for the VR is controlled by the contrast of the 2D views. The transfer function that controls opacity in VR is given by two parameters, window and level: any voxel whose scalar value is less than level -window/2 will be transparent, and any value above level + window/2 will be opaque. The ramp in-between is linear.

Volume Of Interest (VOI) extraction: Manual extraction of a VOI in VR is a desirable feature. Clinicians may need to have an insight into the patient's brain and remove only one block from his head by visually clicking on the screen. This is done in MedINRIA using a cropping box: it consists of a white box that can be scaled and translated by clicking on white handles present on each face (Fig. 2 right). The volume inside the box is automatically removed and cropping is done in real-time (hardware-accelerated).

DTI Track: Log-Euclidean DT-MRI Processing

The DTI Track application is dedicated to DT-MRI processing and visualization. This imaging modality measures the anisotropy of water molecules diffusion within the brain [START_REF] Basser | MR diffusion tensor spectroscopy and imaging[END_REF]. From several MRI measures (at least 6), one is able to reconstruct a diffusion tensor at each voxel. This tensor, represented by a 3D ellipsoid, gives an insight on the directionality of the underlying tissue. If water molecules are free to move in every direction of space, the diffusion tensor will be a sphere. By contrast, if the diffusion is restricted by structures (myelinated axons of the brain white matter), the tensor will be cigar-shaped and aligned along the oriented tissues. This is a unique modality to assess in vivo the directionality of brain white matter, and can lead eventually to the reconstruction of nervous fibers via a process called tractography [START_REF] Weinstein | Tensorlines: Advectiondiffusion based propagation through diffusion tensor fields[END_REF].

DTI Track guides the user through the whole process, from importing the measures to the reconstruction of nervous fibers. The processing pipeline consists of three steps. First, the diffusion tensor field is estimated from the raw data: this generates an image of diffusion tensors that will be used for fiber reconstruction. Second, the diffusion tensor is smoothed using Log-Euclidean metrics [START_REF] Arsigny | Log-Euclidean metrics for fast and simple calculus on diffusion tensors[END_REF]: these metrics have proved to be very efficient practically and are among state-of-theart methods for tensor processing [START_REF] Fillard | Clinical DT-MRI estimation, smoothing and fiber tracking with Log-Euclidean metrics[END_REF]. In third comes the tractography. This last step generates thousands of nervous fibers that requires specific visualization and interaction techniques. Special attention was taken to make these two last steps as fast as possible, not to lower the overall reactivity of the software. This includes a multithreaded implementation with optimized numerical libraries. In the following, we will only focus on the manipulation and analysis of these fibers.

Fiber Manipulation: Depending on the image resolution, tractography can generate up to 100 000 fibers with an average of 100 points each. The task is not easy: one should display a maximum of information while keeping a realtime rendering. We chose to: 1. display fibers using only lines and not more complex geometric primitive like tubes, and 2. introduce a cropping box that intrinsically defines a VOI and only fibers that go through it will be displayed (Fig. 3). The cropping box can be scaled and translated with handles on each of its sides. Moving and scaling the box are real-time actions, giving clinicians a good insight of the brain connectivity.

Interactive Fiber Bundling: Navigating randomly into the reconstructed neural fibers may not be sufficient. Experts might be willing to identify and label a specific fiber bundle of interest (BOI). For this purpose, we offer two solutions. The first one is to manually define a region of interest (ROI) slice by slice on a 3D image. When this image is the patient anatomy (T1), critical structures can be easily spotted and segmented. Accessing the exact bundle that goes through these structures can be of interest. The second method relies on the cropping box. It can be used recursively, i.e., multiple VOIs can be defined by positioning it at various locations. Thus, one can easily isolate a single fiber using this method. Note that both methods can be used alternatively. Once the user is satisfied with the BOI, he can label it, color it and save it for further analysis. Quantitative values of the bundle can be computed: apparent diffusion coefficient, fractional anisotropy [START_REF] Westin | Processing and visualization for diffusion tensor MRI[END_REF], etc.

TensorViewer: Tensor Fields Visualization

Visualizing diffusion tensor fields produced by DTI Track is important for qualitycontrol. In DTI, it is frequent that tensor fields are flipped in the X, Y or Z direction. This is caused by a misalignment of the acquisition parameters with the patient acquisition frame. TensorViewer allows to visually control and correct for flipped tensors (Fig. 4 bottom). In TensorViewer, tensor fields are considered as 3D images, i.e. three slices in three orthogonal directions (namely axial, coronal and sagittal) are displayed (Fig. 4 top). However, one may question whether using ellipsoids for reprensenting tensors is the best choice in terms of speed of refresh rate.

Which tensor representation to choose?

The choice of the geometric primitive to represent a tensor is an active area of research [START_REF] Kindlmann | Superquadric tensor glyphs[END_REF]. As we said earlier, a tensor is fully represented by an ellipsoid: it can be divided into three eigenvectors (axis of the ellipsoid) and eigenvalues (axis magnitude). However, ellipsoids require lots of polygons to be rendered nicely (about 64 polygons). For instance, a dense field of 128 × 128 × 60 voxels (typical resolution of DTI) gives a total of 2 millions polygons to display three orthogonal slices, by far too many for a real-time rendering on a regular computer. To fasten the rendering, we give the possibility to downsample the field and lower the number of polygons. Moreover, we propose to represent a tensor not only by an ellipsoid, but by any other geometric primitives that can correctly represent the shape (from a sphere to a cigar-shaped ellipsoid) with fewer polygons. For instance, a cube, which has only 6 faces, gives a good approximation of the true tensor shape (the cube is rotated and scaled w.r.t. the eigenvectors/eigenvalues). To go one step further, if we get rid of the smallest eigenvalue (which is meaningless most of the time in DTI), planar shapes are possible too: a disk or a square are potential primitives.

Eventually, one can simply represent the principal eigenvector of a tensor instead of the full tensor itself: in DTI, it is admitted that the main eigenvector is aligned with the underlying oriented structure under specific conditions. These conditions are when the fractional anisotropy (FA) [START_REF] Westin | Processing and visualization for diffusion tensor MRI[END_REF] is high. FA tells how far a tensor is from the sphere: 0 means that the tensor is a sphere (absence of structure), 1 means that a tensor is cigar-shaped (presence of structure). In that case, a simple line is an acceptable shape for displaying tensors. MedINRIA proposes a listbox where all geometric shapes discussed above can be selected.

In addition to the shape chosen, we color code a tensor by its main eigenvector (absolute values of the coordinates are considered as RGB): red is for a left-right oriented tensor, green is for antero-posterior and blue is for infero-superior. The color magnitude is controlled by the FA: a dark tensor means absence of fibers. Conversely, a bright tensor means presence of fibers, and the color indicates the principal direction of the structure, according to the color code above. This provides a fast and accurate insight of the locations and directions of brain fibers.

LesionSegmentationEditor: MS Lesion Segmentation Tool

Lesion segmentation is a mandatory step in Multiple Sclerosis (MS) research. However, a precise manual segmentation is a fastidious task and is subject to inter-and intra-expert variability. LesionSegmentationEditor (LSE) aims to ease the segmentation of lesions on conventional brain MRI by experts and to compute meaningful statistics on segmentations (number of lesions, volume).

A lesion can have a different pattern and contrast depending on the MRI modality (T1, T2, Proton Density), and is evolving between two exam timepoints. Consequently, segmentation of lesions is frequently done on several images. To cover this point, the user can load several images in LSE. Each image can be visualized in a tab containing three 2D views and one 3D view. In these views, interactions are similar to ImageViewer (Sec. 2.1). Moreover, a "segmentation mask" is associated to each image: this is a binary image where a voxel is set to 1 in a presence of a lesion. This allows to overlap segmentations on images.

Practically, the segmentation is done on a window called the "Delineation Tab" (Figure 5 (a)): on the right is a 2D preview for each image (e.g. T1, T2, Proton Density). These views are synchronized to show the same part of the brain at any time. The bigger 2D view on the left is used for the actual segmentation. The image displayed can be either one of the images on the right. Lesions can be segmented either manually or automatically. In the manual mode, the user segments a lesion slice by slice with the mouse. In the automatic mode, the user only needs to spot the lesion and click inside it, and the segmentation is achieved automatically. For each lesion, a unique color can be set.

After segmenting a lesion on a specific image, the segmentation can be "transmitted" to the other images. This allows to segment a lesion on the modality where it has the best contrast and to see the result on any image. Similarly, a segmentation can be added or removed from any previous segmentation. Obtained segmentation can be visualized in 2D and 3D (Figure 5 (b)). They can be saved and statistics can be computed (Figure 5 (c)).

Conclusion

The software MedINRIA presented in this paper is a collection of applications dedicated to medical image processing and visualization. The goal was to pro-vide to clinicians an intuitive, reactive and powerful processing and visualization system. For instance, MedINRIA proposes to extract a 3D volume of interest of an image interactively in real-time using a 3D box and volume rendering techniques. In addition, it introduces an interactive method to extract of a neural fiber bundle of interest among thousands of fibers, using the same type of 3D box. MedINRIA does not only to offer user-interactions for manipulating 3D data, but it also provides state-of-the-art algorithms that are most of the time not accessible to clinicians, as these algorithms are generally available in more research software. Here, we propose to process DT-MRI using Log-Euclidean metrics, to register images using a recent diffeomorphic implementation of the demons algorithm, to automatically segment lesions in MS, and to visualize tensors using a panel of geometric shapes. The choice of the methods and their implementation was driven by the constraint of efficiency.

Future work includes the development of new applications: we are currently working on a registration tool that provides a full set of automatic image registration methods, a tumor growth simulation software, and cardiac image processing. MedINRIA is growing fast, and will remain free. A Windows, Linux and MacOSX version can be downloaded at: http://www-sop.inria.fr/asclepios/software/MedINRIA.

Fig. 2 .

 2 Fig. 2. Image visualization with MedINRIA. Left: The preview screen. Right: Image rendering in 3D. Top: Multiplanar reconstruction (MPR) of a T1 image. An axial, sagittal and coronal slices is displayed. Middle: Volume rendering (VR) of the same image. Bottom: The cropping box was used to extract a VOI. By clicking on the handles, it can be resized (face handles) and translated (center handle).

Fig. 3 .

 3 Fig. 3. Example of fiber bundling with MedINRIA. Top left: The complete set of reconstructed fibers. The white box displayed is the cropping box. Only fibers that go through that box are shown. Top right: The box was shrunk by clicking on the face handles. The red arrows show the direction of the shrinking. The user can choose to navigate only in this sub-sample of fibers. Bottom right: The box was moved to another position, thus defining a second VOI. The bundle of interest (BOI) is ready to be tagged (i.e. labeled and saved). Bottom left: The BOI is now tagged and colored in red. It can be saved for further analysis.

Fig. 4 .

 4 Fig. 4. Tensor visualization with MedINRIA. Top: An axial (left), coronal (middle) and sagittal (right) slice of a 3D tensor field. Bottom: The tensor field (left) looks flipped left-right. We applied a flipping (right) to correct this defect.

  Fig. 5. LSE Interaction, Visualization and Results


Chengyi Shen

Daniel Brito

Julien Diaz

Deyuan Zhang

Jochen Schell

Clarisse Bordes

Stéphane Garambois

High resolution spatial and temporal laboratory seismic datasets by Laser Doppler Vibrometry

High resolution spatial and temporal laboratory seismic datasets by Laser Doppler Vibrometry

Research context & objectives

We propose to perform sub-metric scale seismic measurements with innovative experimental tools in a laboratory environment which reproduces large-scale field explorations in well conditioned and controlled environment. The purpose is to develop high resolution seismic methods on various natural samples, that can be transferred later to large-scale field conditions. Seismic waves are produced in our experiments by a P-wave piezoelectric transducer stuck to the sample surface. Nanometer mechanical displacements induced by the transducer are measured by LDV around the sample. The data acquisition chain is tested in an aluminum cuboid of 280 mm side length. High-resolution (ds ≈ 1 mm, dt ≈ 10 ns) unidimensional mappings of surface displacements on various faces of the cuboid are performed ; Furthermore, multidimensional datasets are obtained on both cuboid and cylinder samples by full-field 3D Scanning Vibrometers thanks to the collaboration with Polytec GmbH. Parallel 2D/3D numerical simulations using IPDG [START_REF] Baldassari | Modélisation et simulation numérique pour la migration terrestre par équation d'ondes[END_REF] discretization method are done to match the experimental data.

LDV and geophysical applications

Laser Doppler Vibrometer

Displacement mapping on an aluminum core

Radial components of displacements are measured. Snapshots and seismogram (around the mid-section) at/until t1, t2, t3, t4 (t1 < t2 < t3 < t4) are shown.

Corresponding simulations

Corresponding 2D and 3D simulations by Finite Element method featured with IPDG [START_REF] Baldassari | Modélisation et simulation numérique pour la migration terrestre par équation d'ondes[END_REF]. The 3D simulation is performed in a reduced model with Absorbing Boundary Conditions due to memory limits.

Multidimensional measurements

Seismogram from the flat-surface dataset

Seismogram extracted from raw data. The out-of-plane components (U z ) are relatively strong on all the 3 receiver lines because the seismic source points to this direction. The relative energy of the in-plane components (U x and U y ) depends on the receiver line orientation. Along the receiver line at 0 • (X direction), U x predominates, while at 90 • (Y direction), U y predominates. U x and U y are equally distributed along the 45 • receiver line. These quantitative observations are coherent with the symmetry of wave propagation in homogeneous media. The full-field 3D scanning vibrometer measures evidently the out-of-plane components as well as the in-plane components.

On an Alu cylinder surface (Wavefield snapshots)

Snapshots of wavefields (raw data without filtering the ambient noise) at 8.90 µs, 18.50 µs, 19.78 µs and 36.42 µs, with the source at the middle of the gap. On the second snapshot we can clearly see the symmetric P wave fields in the front, followed by slower S waves. The different polarization of S waves are highlighted on the last one.

Application on carbonate core

Experimental setup and wavefield snapshots Snapshots at 21.76 µs for each source position. 8 uniformly distributed sources are deployed around the middle section of the core. Measurements close to the sources are not reliable because of instrumentation obstacles. The wavefield vectors are no longer symmetric in space due to the heterogeneities inside the natural rock. We can still recognize the P wavefield and S wavefield, however, they are strongly disturbed in the cases of source 6 and source 7, indicating probably an obstacle inside the core near these two sources.

Seismogram and X-ray tomography Conclusion :

The excellent agreement between experiments and simulations demonstrates that the full wavefield is nicely recorded by the single point LDV, reciprocally validating our numerical scheme and models. A first qualitative analysis on the multidimensional data confirmed the data quality and the capability of the full-field 3D scanning vibrometer. We can then perform basic time-lapse tomography in order to get a first insight of the elastic wave velocity models for the carbonate core. Further analyses on simulated and experimental multidimensional data will enable us to look quantitatively into the polarizations and quantify the Amplitude vs Angle (AVA [START_REF] John | Offset-dependent reflectivity -Theory and practice of AVO analysis[END_REF]) influence on the surfacic measurements. Seismic attributes such as amplitudes and frequencies contain also rich information on other rock physics properties including quality factor, anisotropy and proelasticity etc. [START_REF] Dupuy | Estimation of rock physics properties from seismic attributes ?part 2 : Applications[END_REF]. The full-field data can eventually be used to test Full Waveform Inversion schemes which would in turn yield more reliable and accurate inversion results.

C. Shen 1 , 3

 13 , D. Brito 1,3 , J. Diaz 2,3 , D. Zhang 1 , J. Schell 4 , C. Bordes 1 , S. Garambois 5 1) CNRS/ TOTAL / Univ Pau & Pays Adour/E2S UPPA, Laboratoire des Fluides Complexes et leurs Réservoirs -IPRA, UMR5150, 64000, Pau, FRANCE 2) Univ Pau & Pays Adour/CNRS, Laboratoire de Mathématiques et de leurs Applications, UMR5142, 64000, Pau, FRANCE 3) Project MAGIQUE-3D, INRIA Bordeaux-Sud-Ouest, 64013 Pau, France 4) Polytec GmbH, Polytec-Platz. 1-7, 76337 Waldbronn, Germany 5) Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, IFSTTAR, ISTerre, UMR5275, 38000 Grenoble, France

  a) : Schematic view of the LDV principle ; b) : Experimental setup of measurements on an Alu block by a single point LDV sensor. c) : Experimental setup for the Alu cylinder. Displacement mapping on an aluminum block Normal components of displacements are measured on face 1 (source face), face 2 (upper face) and face 6 (opposite face). Snapshots and seismogram (along middle lines) at/until 19.3 µs, 34.4 µs, 49.4 µs and 64.3 µs are displayed.

  On a flat surface (280 × 280 × 100 mm 3 Alu block) Snapshots of the displacements calculated by integrating the velocity measurements (raw data without filtering the ambient noise) at 19.07 µs, 22.72 µs, 36.42 µs and 48.90 µs. The first two are highlighted by P wave fields and the laser two by S waves. We can clearly see the different polarizations respectively for P and S waves.


Cedric Bernardin
email: cbernard@unice.fr

AND François Huveneers
email: huveneers@ceremade.dauphine.fr

Stefano Olla
email: olla@ceremade.dauphine.fr

Cédric Bernardin

HYDRODYNAMIC LIMIT FOR A DISORDERED HARMONIC CHAIN

come

INTRODUCTION

In this paper, we consider a one-dimensional chain of coupled harmonic oscillators with masses m x and Hamiltonian

H = x∈ p 2 x 2m x + g (q x+1 -q x ) 2 2 .
By changing units, one can assume that the stiffness coefficient g is equal to 1. The dynamics is governed by Hamilton's equations:

m x qx = p x , ṗx = (∆q) x ,
where we have used the notation ∆ = ∇ -∇ + = ∇ + ∇ -for the discrete Laplacian, with

(∇ + f ) x = f x+1 -f x and (∇ -f ) x = f x -f x-1 .
For the sake of simplicity, we consider for now the system on the infinite lattice . This is by no means necessary, and starting from the next section, we will restrict ourselves to a finite box and take the thermodynamic limit properly. This system was first analyzed in finite volume when all masses m x are equal. Putting the chain in a non-equilibrium stationary state (NESS) between two heat reservoirs at different temperatures, it was found in [START_REF] Rieder | Properties of a harmonic crystal in a stationary nonequilibrium state[END_REF] that the energy current does not decay with the size of the system, indicating that energy propagates ballistically. The situation changes if the masses are taken to be i.i.d. random variables. This case was first investigated in [START_REF] Rubin | Abnormal Lattice Thermal Conductivity of a One-Dimensional, Harmonic Isotopically Disordered Crystal[END_REF][START_REF] Casher | Heat Flow in Regular and Disordered Harmonic Chains[END_REF] and subsequently studied in [START_REF] Verheggen | Transmission Coefficient and Heat Conduction of a Harmonic Chain with Random Masses: Asymptotic Estimates on Products of Random Matrices[END_REF][START_REF] Dhar | Heat Conduction in the Disordered Harmonic Chain Revisited[END_REF][START_REF] Ajanki | Rigorous Scaling Law for the Heat Current in Disordered Harmonic Chain[END_REF]. As it turns out, the disordered harmonic chain is an Anderson insulator in disguise [START_REF] Anderson | Absence of Diffusion in Certain Random Lattices[END_REF]. However, as a consequence of the conservation of momentum, the ground state of the operator M -1 ∆, featuring in Newton's equation q = M -1 ∆q with M the diagonal matrix of the masses, is a "symmetry protected mode" [START_REF] Halperin | Quantized Hall conductance, current-carrying edge states, and the existence of extended states in a two-dimensional disordered potential[END_REF], implying a divergent localization length in the lower edge of the spectrum. This leads to a rich and unexpected phenomenology. In particular, if the chain is again in a NESS, the scaling of the energy current with the Date: January 7, 2019. Work partially supported by the grants ANR-15-CE40-0020-01 LSD and ANR-14-CE25-0011 EDNHS of the French National Research Agency (ANR).

system size happens to depend on boundary conditions and spectral factors of the reservoirs [START_REF] Dhar | Heat Conduction in the Disordered Harmonic Chain Revisited[END_REF]. This finding reveals also the complete lack of local thermal equilibrium, that results eventually from integrability (see Section 5.2).

The harmonic chain has three "obvious" conserved quantities: the total energy H, the total momentum P = x p x and the total stretch or elongation R = x r x with r x = (∇ + q) x . This gives rise to the following microscopic conservation laws:

ṙx = p x+1 m x+1 - p x m x , ṗx = r x -r x-1 , ėx = r x p x+1 m x+1 - r x-1 p x m x with e x = 1 2 p 2 x m x + r 2
x . After a hyperbolic rescaling of space and time, we ask in this paper whether the empirical densities of these conserved quantities converge to the densities r, p and e governed by the macroscopic laws Instances of rigorous derivation of Euler equations in the smooth regime rest on the ergodicity of the microscopic dynamics. In [START_REF] Braxmeier-Even | Hydrodynamic Limit for a Hamiltonian System with Boundary Conditions and Conservative Noise[END_REF][START_REF] Olla | Hydrodynamic limit for a Hamiltonian system with weak noise[END_REF][START_REF] Komorowski | Ballistic and superdiffusive scales in the macroscopic evolution of a chain of oscillators[END_REF], the Hamiltonian dynamics is perturbed by some stochastic noise acting in such a way that conserved quantities are not destroyed but that the ergodicity of the dynamics can be established rigorously.

∂ t r
One of the main motivations of this work is to show that ergodicity is not in general a necessary assumption for Euler equations to hold. Indeed, the dynamics considered here is purely Hamiltonian, non-ergodic, and possesses actually a full set of invariant quantities (see Section 5.2). In the clean case, i.e. when the masses are all equal, we show in Section 1.1 that Euler equations hold if and only if the temperature profile is constant. Instead, in Section 1.2, we argue that Euler equations hold even out of thermal equilibrium if there is disorder on the masses. We briefly discuss the fate of other conserved quantities in Section 1.3. Theorem 1 in Section 2 constitutes the main result of our paper: We show the convergence to Euler equations for the disordered harmonic chain, almost surely with respect to the masses and on average with respect to an initial local Gibbs state. The rest of the paper is devoted to the proof of this theorem.

1.1. Clean harmonic chain. Let us assume that all masses m x are equal, say m x = 1 for simplicity. In this case, the equations of motion read

qx = p x , ṗx = ∆q x . (1.1)
Let us first consider the thermal equilibrium case: Assume that the initial configuration of the chain is random and distributed according to a Gaussian law µ 0 with covariance matrix

〈 〈 〈 (∇ + q) x ; (∇ + q) y 〉 〉 〉 = 〈 〈 〈 p x ; p y 〉 〉 〉 = β -1 δ x, y , 〈 〈 〈 q x ; p y 〉 〉 〉 = 0, (1.2)
for some inverse temperature β . It is easy to prove that at time t > 0, the distribution µ t in the phase space is still given by a Gaussian law with the same covariance matrix.

To see this, just use Fourier transforms to diagonalize the dynamics:

q(k, t) = x∈ e i2πkx q x (t), p(k, t) = x∈ e i2πkx p x (t), (1.3)
and define the wave function

φ(k, t) = ω(k)q(k, t) + i p(k, t) (1.4)
where ω(k) = |2 sin(πk)| is the dispersion relation. Then, the explicit solution of (1.1) is given by φ(k, t) = e -iω(k)t φ(k, 0).

(1.5)

The correlations (1.2) imply that

〈 〈 〈 φ(k, 0) * ; φ(k ′ , 0) 〉 〉 〉 = 2β -1 δ(k -k ′ ), 〈 〈 〈 φ(k, 0); φ(k ′ , 0) 〉 〉 〉 = 0. (1.6) Consequently 〈 〈 〈 φ(k, t) * ; φ(k ′ , t) 〉 〉 〉 = e i(ω(k)-ω(k ′ ))t 〈 〈 〈 φ(k, 0) * ; φ(k ′ , 0) 〉 〉 〉 = 2β -1 δ(k -k ′ ), 〈 〈 〈 φ(k, t); φ(k ′ , t) 〉 〉 〉 = e -i(ω(k)+ω(k ′ ))t 〈 〈 〈 φ(k, 0); φ(k ′ , 0) 〉 〉 〉 = 0. (1.7)
From (1.6) and (1.7), we deduce that the covariances (1.2) are the same at any time t:

〈 〈 〈 (∇ + q) x (t); (∇ + q) y (t) 〉 〉 〉 = 〈 〈 〈 p x (t); p y (t) 〉 〉 〉 = β -1 δ x, y , 〈 〈 〈 q x (t); p y (t) 〉 〉 〉 = 0, (1.8)
which implies that the Gaussian distribution µ t differs from µ 0 only by the averages rx

(t) = 〈 〈 〈 r x (t) 〉 〉 〉 = 〈 〈 〈 (∇ + q) x (t) 〉 〉 〉 and px (t) = 〈 〈 〈 p x (t) 〉 〉
〉 that, by linearity of the dynamics, evolve following the same equation (1.1).

Assume now that the initial averages of the momentum p x and stretch r x = (∇ + q) x are slowly varying on a macroscopic scale. More precisely, let N be an integer representing a macroscopic number of sites in the chain, and let p, r : → be smooth and fast decaying initial macroscopic profiles. We let

r[N y] (0) = r( y), p[N y] (0) = p( y).
(1.9)

Let p(ξ) and r(ξ) be the Fourier transforms (in ) of p( y) and r( y). Then, as N → ∞,

1 N p( ξ N ) -→ p(ξ) and 1 N r( ξ N ) -→ r(ξ).
After a straightforward analysis we have that 1

N p ξ N , N t -→ p (ξ, t), 1 N r ξ N , N t -→ r (ξ, t) (1.10)
where

∂ t r (ξ, t) = -i2πξ p (ξ, t), ∂ t p (ξ, t) = -i2πξ r (ξ, t). (1.11)
Consequently r[N y] (N t) and p[N y] (N t) converge (as distributions) to the solution of the linear wave equation

∂ t r ( y, t) = ∂ y p ( y, t), ∂ t p ( y, t) = ∂ y r ( y, t).
(1.12)

Let us now consider the energy per particle

e x = 1 2 p 2 x + r 2 x . Its average under the distribution µ t is 〈 〈 〈 e x (t) 〉 〉 〉 = β -1 + 1 2 p2 x (t) + r2
x (t) since by (1.8), the variance of p x and r x , i.e. the temperature, remains constant in time. In the limit N → ∞ we have

〈 〈 〈 e [N y] (N t) 〉 〉 〉 -→ e ( y, t) = β -1 + 1 2 p 2 ( y, t) + r 2 ( y, t) ,
i.e. it solves the equation

∂ t e ( y, t) = ∂ y (p ( y, t) r ( y, t)) . (1.13)
We recognize that (1.12 -1.13) are the Euler equations. The above is the simplest example of propagation of local equilibrium and hydrodynamic limit in hyperbolic scaling: in a harmonic chain in thermal equilibrium at temperature β -1 , and the mechanical modes not in equilibrium, we will have a persistence of the thermal equilibrium at any time t, while the mechanical modes evolve independently from the thermal mode following the linear wave equation.

Notice that the argument above does not require the distribution µ 0 to be the thermal equilibrium measure defined by (1.2), and that it holds for any measure µ 0 with translation invariant covariance given by

〈 〈 〈 ∇q x ; ∇q y 〉 〉 〉 = 〈 〈 〈 p x ; p y 〉 〉 〉 = C(x -y), 〈 〈 〈 q x ; p y 〉 〉 〉 = 0 (1.14)
for a positive definite function C(x). The only difference is then that in (1.6-1.7) the term β -1 has to be replaced by the Fourier transform C(k). Actually, the measure µ t is not even a local equilibrium state [START_REF] Bernardin | Thermodynamics and non-equilibrium macroscopic dynamics of chains of anharmonic oscillators[END_REF], underlining that the validity of Euler equations in this example does not require the propagation of local equilibrium.

The above argument rests on the translation invariance of the distribution of the thermal modes and fails if it is space inhomogeneous, for example if the starting distribution is given by a local Gibbs state with a slowly varying temperature β -1 N ,x = β -1 (x/N ), i.e. a Gaussian measure with covariances

〈 〈 〈 ∇q x ; ∇q y 〉 〉 〉 = 〈 〈 〈 p x ; p y 〉 〉 〉 = β -1 N ,x δ x, y , 〈 〈 〈 q x ; p y 〉 〉 〉 = 0. (1.15)
In this case, even though the wave equation (1.12) still holds, generally the energy equation (1.13) is not valid. In fact the energies of each mode k evolves autonomously, as we can see studying the limit evolution of the Wigner distribution defined by

W N (ξ, k, t) := 2 N φ * k - ξ 2N , N t φ k + ξ 2N , N t W N ( y, k, t) := e -i2πξ y W N (ξ, k, t) dξ, (1.16)
(the above definitions should be understood as distributions on × Π with Π = \ ).

In the limit as N → ∞ the Wigner distribution converge to a positive distribution with an absolutely continuous part, the local distribution of the thermal modes, and a singular part concentrated on k = 0, the mechanical modes:

lim N →∞ W N (ξ, k, t) = W th (ξ, k, t) + W m (ξ, t) δ 0 (d k) (1.17)
The mechanical part W m (ξ, t) is the Fourier transform of 1 2 p 2 ( y, t) + r 2 ( y, t) . A straightforward calculation gives for the thermal part (see [START_REF] Dobrushin | One-dimensional Harmonic Lattice Caricature of Hydrodynamics[END_REF] or [START_REF] Basile | Energy Transport in Stochastically Perturbed Lattice Dynamics[END_REF] for a rigorous argument):

W th (ξ, k, t) = e -iω ′ (k)ξt W th (ξ, k, 0). (1.18) This implies that the inverse Fourier transform W th ( y, k, t) satisfies the transport equation

∂ t W th ( y, k, t) + ω ′ (k) 2π ∂ y W th ( y, k, t) = 0. (1.19)
It also follow that

W th ( y, k, t) d k = ẽ( y, t) (1.20)
where ẽ( y, t) is the limit profile of thermal energy (or temperature) defined as

1 2 〈 〈 〈 r [N y] (N t); r [N y] (N t) 〉 〉 〉 + 〈 〈 〈 p [N y] (N t); p [N y] (N t) 〉 〉 〉 ẽ( y, t). (1.21)
Consequently the thermal energy ẽ( y, t) evolves non autonomously following the equation

∂ t ẽ( y, t) + ∂ y J ( y, t) = 0, J ( y, t) = ω ′ (k)W th ( y, k, t) d k. (1.22)
We say that the system is in local equilibrium if W th ( y, k) = β -1 ( y) constant in k. This correspond to the fact that Gibbs measure gives uniform distribution on the modes. Starting in thermal equilibrium means W th ( y, k, 0) = β -1 and trivially W th ( y, k, t) = β -1 for any t > 0. But starting with local equilibrium, i.e. W ( y, k, 0) = β -1 ( y) constant in k, we have a non autonomous evolution of ẽ( y, t).

1.2. Disordered harmonic chain. The situation so far can be summarized as follows.

By linearity, the variables r x and p x admit a macroscopic limit described by (1.12) independently of the initial temperature profile. The macroscopic equation (1.13) predicts that the evolution of the energy is purely mechanical and that the temperature does not evolve with time. As it turns out, the evolution of the mechanical energy is correctly described by Euler equation (see the term N (t) in our decomposition (4.2) below), but thermal fluctuations do in general evolve with time as well, except if the temperature profile is initially flat. This picture gets strongly modified if the masses are taken to be random. On the one hand, deriving the macroscopic evolution of the fields r x and p x becomes less obvious because some homogenization over the masses is required. This difficulty can be solved by the elegant method of the "corrected empirical measure", see [START_REF] Gonçalves | Scaling Limits for Gradient Systems in Random Environment[END_REF][START_REF] Jara | Quenched non-equilibrium central limit theorem for a tagged particle in the exclusion process with bond disorder[END_REF][START_REF] Bernardin | Homogenization results for a linear dynamics in random Glauber type environment[END_REF] (though we will actually solve it another way). On the other hand, and this is the main point in considering random masses, the evolution of the energy e x is now much better approximated by Euler equation. Indeed, at a microscopic level, all thermal fluctuations are frozen thanks to Anderson localization and the evolution of the energy becomes purely mechanical.

To understand this a little bit better, it is good to realize how the disorder modifies the nature of the eigenmodes (ψ k ) 1≤k≤N of the operator M -1 ∆ for a finite chain of size N . As a consequence of Anderson localization [START_REF] Anderson | Absence of Diffusion in Certain Random Lattices[END_REF], all modes at positive energy are spatially localized. However the localization length ζ k diverges as one approaches the ground state:

ζ -1 k ∼ ω 2 k ∼ k N 2 ,
so that only the modes with k N are actually localized, while the modes k N remain comparable to the modes of the clean chain [START_REF] Verheggen | Transmission Coefficient and Heat Conduction of a Harmonic Chain with Random Masses: Asymptotic Estimates on Products of Random Matrices[END_REF][START_REF] Ajanki | Rigorous Scaling Law for the Heat Current in Disordered Harmonic Chain[END_REF]. By imposing a smooth initial profile r, p, the initial local Gibbs state attributes a weight of order 1 to a few first modes above the ground state, and a weight of order 1 to all other modes together. The first ones are responsible for the transport of mechanical energy; all modes with k ≫ N are localized and do not transport any thermal energy; all modes with 1 ≪ k ≤ o(N ) have a vanishing weight in the thermodynamic limit and can be neglected in the analysis.

Finally, we would like to mention that, while the disorder considered here and the stochastic velocity exchange noise considered in [START_REF] Braxmeier-Even | Hydrodynamic Limit for a Hamiltonian System with Boundary Conditions and Conservative Noise[END_REF][START_REF] Komorowski | Ballistic and superdiffusive scales in the macroscopic evolution of a chain of oscillators[END_REF] act in an obviously very different way, e.g. the disorder preserves integrability while the stochastic noise makes the dynamic ergodic, they do produce the same effects in some respect. Indeed the noise has only a very slow (negligeable) effect on the macroscopic modes. This bares some similarity with the fact that the disorder has very little influence on the low modes of the disordered chain, while the dynamical noise provides an active hopping mechanism among the high modes. Consequently the dynamical noise produces a superdiffusive sub-ballistic spreading of the thermal energy [START_REF] Jara | Superdiffusion of Energy in a Chain of Harmonic Oscillators with Noise[END_REF][START_REF] Komorowski | Ballistic and superdiffusive scales in the macroscopic evolution of a chain of oscillators[END_REF], that is not visible in the hyperbolic scaling. Thus, the dynamical noise plays here as well a role analogous to the disorder only in the hyperbolic scaling by freezing the temperature profile.

It is important to notice that in the disordered case the temperature profile remains frozen at any time scale, including the diffusive time scale and further, see remark 3. In particular this implies a vanishing thermal diffusivity for the disordered unpinned harmonic chain. This is not in contradiction with the divergence of the thermal conductivity observed in the NESS of the same unpinned system when connected to Langevin thermostats at different temperatures with free boundary conditions, see [START_REF] Casher | Heat Flow in Regular and Disordered Harmonic Chains[END_REF][START_REF] Ajanki | Rigorous Scaling Law for the Heat Current in Disordered Harmonic Chain[END_REF], and our result sheds actually some light on the various behaviors for the conductivity found in [START_REF] Dhar | Heat Conduction in the Disordered Harmonic Chain Revisited[END_REF] depending on the boundary conditions. In fact the thermal conductivity divergence in the NESS is due to the fluctuations of the low mechanical modes, and in our analysis there is a clear separation of the behavior of the mechanical modes (responsable for the ballistic motion) and the high thermal modes that give the temperature profile.

1.3. Other conserved quantities. Before moving on, let us briefly comment on the issue of the other conserved quantities of the system. These can also be written as a sum of local terms and lead thus to additional conservation laws. For example,

I = x d x = 1 2 x (r x -r x-1 ) 2 m x + p x+1 m x+1 - p x m x 2
is conserved (see Sections 3.1 and 5.2) and leads to the microscopic conservation law

ḋx = p x+1 m x+1 - p x m x r x+1 -r x m x+1 - p x m x - p x-1 m x-1 r x -r x-1 m x .
It is thus natural to ask whether this relation generates also some macroscopic law. In the cases where we can derive the macroscopic evolution equation (1.13) for the energy, it is easy to argue that the corresponding macroscopic density d( y, t) does not evolve with time in the hyperbolic scaling. Indeed, we can decompose d x as the sum of a mechanical and a thermal contribution, as we do in (4.2) below for the energy. In this case, contrary to what happens for the energy, the mechanical contribution vanishes in the thermodynamic limit since d x depends on r and p only through their gradients, while the contribution from the thermal modes does not evolve with time, for the same reasons as it does not for the energy.

All the other conserved quantities in this model that can be written as a sum of local terms are obtained by taking further gradients in the variables r and p (see Section 5.2), and have thus no evolution either in the hyperbolic scaling.

MODEL AND RESULTS

We define the model studied in this paper and we state our main result. For technical reasons, it is easier to work on a finite system of size N and then let N → ∞.

Hamiltonian model. The Hamiltonian H on 2N is defined by

H(q, p) = 1 2 N x=1 p 2 x m x + ((∇ + q) x ) 2 .
For concreteness, we assume free boundary conditions, i.e. q 0 = q 1 and q N +1 = q N ; other boundary conditions such as fixed or periodic could be considered just as well.

The masses (m x ) 1≤x≤N are i.i.d. random variables. In order to avoid any technical difficulty in exploiting known results from the Anderson localization literature, we assume that the law of m x admits a smooth density compactly supported in [m -, m + ] with m -> 0.

The equations of motion read M q = p and ṗ = ∆q where M is the square diagonal matrix of size N with entries defined by M x, y = δ(xy)m x (δ(z) is defined by 1 for z = 0 and 0 otherwise). It is more convenient to express the equations of motion in terms of the displacement variables

r x = (∇ + q) x (1 ≤ x ≤ N -1).
The equations of motion become

ṙx = ∇ + M -1 p x (1 ≤ x ≤ N -1), ṗx = (∇ -r) x (1 ≤ x ≤ N ) (2.1)
where we use fixed boundary conditions for r in the second equation:

r 0 = r N = 0.
2.2. Gibbs and locally Gibbs states. We consider three locally conserved quantities in the bulk:

H = N x=1 e x = N x=1 p 2 x 2m x + r 2 x 2 , P = N x=1 p x , R = N -1 x=1 r x .
The energy H and the momentum P are actually truly conserved, but the conservation of R is broken at the boundary:

Ṙ = m -1 N p N -m -1 1 p 1 .
The Gibbs states are characterized by three parameters: β > 0 and p, r ∈ . Its probability density writes

ρ G (r, p) = 1 Z G exp - β 2 N x=1 m x p x m x - p m 2 - β 2 N -1 x=1 (r x -r) 2 .
where m denotes the mean mass and Z G := Z G (β , p, r) is a normalizing constant. Local Gibbs states are obtained by replacing the constant parameters β , p, r by functions

β , p, r : [0, 1] → , with β (x) > 0 for all x ∈ [0, 1],
and by considering the measure with density

ρ loc (r, p) = 1 Z loc exp - 1 2 N x=1 β (x/N )m x p x m x - p(x/N) m 2 - 1 2 N -1 x=1 β (x/N )(r x -r(x/N )) 2 (2.2)
where Z loc := Z loc (β , p, r) is a normalizing constant. We impose the following regularity conditions on β , p, r:

β ∈ 0 ([0, 1]), r ∈ 1 ([0, 1]) with r(0) = r(1) = 0, p ∈ 1 ([0, 1]). (2.3)
We take such a local Gibbs state as initial state. Below, we denote the expectation with respect to it by

〈 〈 〈 • 〉 〉 〉: 〈 〈 〈 F 〉 〉 〉 = F (r, p)ρ loc (r, p) drdp.
Instead, expectation (resp. probability) with respect to the masses is denoted by E (resp. P).

Evolution of the locally conserved quantities.

Let us fix some maximal time T > 0. Let us define the fields , and acting on functions f ∈ 0 ([0, 1]) as

( f , t) = 1 0 r( y, t) f ( y) d y, ( f , t) = 1 0 p( y, t) f ( y) d y, ( f , t) = 1 0 e( y, t) f ( y) d y.
(2.4) for all t ∈ [0, T ]. The kernels r, p and e are defined as follows. First, at t = 0, we impose

r( y, 0) = r( y), p( y, 0) = p( y), e( y, 0) = 1 β ( y) + p 2 ( y) 2m + r 2 ( y) 2 .
Next, the evolution at all further time is governed by the following system of conservation laws:

∂ t r( y, t) = 1 m ∂ y p( y, t), r(0, t) = r(1, t) = 0, (2.5)
∂ t p( y, t) = ∂ y r( y, t), (2.6) ∂ t e( y, t) = 1 m ∂ y (r( y, t)p( y, t)). ( 2.7)
Thanks to the regularity conditions on r, p in (2.3), the solutions of these equations are classical. Since (r, p) are solution of wave equations with suitable boundary conditions, they can be obtained explicitly by expanding them in Fourier series. Then, by a time integration, e may be expressed as a function of (r, p), see (4.1). Later we will use that a classical solution for the system governing (r, p) coincides with the (unique) weak solution of this system. Because of the boundary conditions, test functions will have to be chosen appropriately (see (3.11-3.12)).

Theorem 1. Let t ∈ [0, T ] and f ∈ 0 ([0, 1]). Let us assume that the system is initially prepared in a locally Gibbs state such that β , r and p satisfy (2.3). Then, as N → ∞, almost surely (w.r.t. P),

N ( f , t) = 1 N N x=1 f (x/N ) 〈 〈 〈 r x (N t) 〉 〉 〉 → ( f , t), (2.8) N ( f , t) = 1 N N x=1 f (x/N ) 〈 〈 〈 p x (N t) 〉 〉 〉 → ( f , t), ( 2.9
)

N ( f , t) = 1 N N x=1 f (x/N ) 〈 〈 〈 e x (N t) 〉 〉 〉 → ( f , t).
(2.10) Remark 1. As pointed out in the introduction, the situation is much simpler at thermal equilibrium, i.e. for β constant, and these limits hold even for the non-disordered chain. See Section 4.3 for a derivation along the lines used to derive Theorem 1.

EVOLUTION OF N AND N

In this section, we show the limits (2.8-2.9). Moreover, in order to later deal with the field N , we show more:

The functions 〈 〈 〈 r [N y] (N t) 〉 〉 〉 and m -1 [N y] 〈 〈 〈 p [N y] (N t) 〉 〉 〉 are uniformly (in N ) Hölder regular in y ∈ [0, 1]
, with exponent at least 1/2. Hence they converge pointwise to r( y, t) and p( y, t) respectively.

3.1. A priori estimates. Given d ∈ , we denote the standard scalar product on d by 〈•, •〉 d (we will drop the subscript d when no confusion seems possible). Let us consider the two following conserved quantities:

H(r, p) = 1 2 〈p, M -1 p〉 N + 〈r, r〉 N -1 , (3.1)
I (r, p) = 1 2 〈∇ -r, M -1 ∇ -r〉 N + 〈∇ + M -1 p, ∇ + M -1 p〉 N -1 . (3.2)
The conservation of I follows from the fact that, if (r, p) solve (2.1), then (∇ + M -1 p, ∇ -r) solve the same equation, the corresponding Hamiltonian being I (since we have that

H(∇ + M -1 p, ∇ -r) = I (r, p)).
Notice also that a full set of conserved quantities can be generated by further taking gradients, see Section 5.2.

Thanks to these two conservation laws, and to the smoothness assumptions on r and p, we deduce Lemma 1. There exists a deterministic C such that, for any t ≥ 0 and any N ∈ ,

N -1 x=1 〈 〈 〈 r x (N t) 〉 〉 〉 2 ≤ CN , N x=1 〈 〈 〈 p x (N t) 〉 〉 〉 2 ≤ CN , (3.3)
N x=1 〈 〈 〈 (∇ -r) x (N t) 〉 〉 〉 2 ≤ C N , N -1 x=1 〈 〈 〈 (∇ + M -1 p) x (N t) 〉 〉 〉 2 ≤ C N . (3.4)
Proof. By linearity of the equations of motion (2.1),

(〈 〈 〈 r 〉 〉 〉, 〈 〈 〈 p 〉 〉 〉
) solve the same equations as (r, p). Therefore, the conservation of H(r, p) and I (r, p) implies the conservation of

H(〈 〈 〈 r 〉 〉 〉, 〈 〈 〈 p 〉 〉 〉) and I (〈 〈 〈 r 〉 〉 〉, 〈 〈 〈 p 〉 〉 〉).
Since the quantities to be estimated in (3.3) are bounded by

H(〈 〈 〈 r 〉 〉 〉, 〈 〈 〈 p 〉 〉 〉)
and the quantities to be estimated in (3.4) are bounded by

I (〈 〈 〈 r 〉 〉 〉, 〈 〈 〈 p 〉 〉 〉), we conclude that is it enough to establish them respectively for H(〈 〈 〈 r 〉 〉 〉, 〈 〈 〈 p 〉 〉 〉) and I (〈 〈 〈 r 〉 〉 〉, 〈 〈 〈 p 〉 〉 〉) at t = 0.
This follows from a direct computation, thanks to the product structure of the local Gibbs state (2.2) and to the hypotheses on r and p in (2.3) (in particular, this is the place where the boundary condition on r plays a role).

Remark 2. Notice that the bounds in Lemma 1 are actually valid for any time scale N α t, for any α > 0.

As a corollary, we deduce the existence of a constant C ∈ such that, for any x, y ∈ ∩

[1, N ], 〈 〈 〈 r x ′ (N t) 〉 〉 〉 -〈 〈 〈 r x (N t) 〉 〉 〉 ≤ C x ′ -x N 1/2 , m -1 x ′ 〈 〈 〈 p x ′ (N t) 〉 〉 〉 -m -1 x 〈 〈 〈 p x (N t) 〉 〉 〉 ≤ C x ′ -x N 1/2 , ( 3.5)
and therefore also such that

|〈 〈 〈 r x (N t) 〉 〉 〉| ≤ C, |〈 〈 〈 p x (N t) 〉 〉 〉| ≤ C. (3.6)
Indeed, to get e.g. (3.5) for r, we deduce from (3.4) that

〈 〈 〈 r x ′ (N t) 〉 〉 〉 -〈 〈 〈 r x (N t) 〉 〉 〉 = x ′ z=x+1 〈 〈 〈 (∇ -r) z (N t) 〉 〉 〉 ≤ N z=1 〈 〈 〈 (∇ -r) z (N t) 〉 〉 〉 2 1/2 |x -x ′ | 1/2 ≤ C x ′ -x N 1/2 .
Next (3.6) follows from (3.5) if, given N , t, there exists at least some x 0 such that the inequalities hold. This in turn follows from (3.3).

3.2.

Averaging lemma for the field N . The method of the corrected empirical measure is an elegant method to deal with the randomness on the masses in deriving the hydrodynamic limit for N and N [START_REF] Gonçalves | Scaling Limits for Gradient Systems in Random Environment[END_REF][START_REF] Jara | Quenched non-equilibrium central limit theorem for a tagged particle in the exclusion process with bond disorder[END_REF][START_REF] Bernardin | Homogenization results for a linear dynamics in random Glauber type environment[END_REF]. However, in our case, it seems more convenient to use the following lemma:

Lemma 2. Let f ∈ 0 ([0, 1]
) and t ≥ 0. Almost surely (w.r.t. the masses), for N → ∞,

1 N N x=1 f (x/N ) 〈 〈 〈 p x (N t) 〉 〉 〉 m x (m x -m) → 0, (3.7)
1

N N x=1 f (x/N ) 〈 〈 〈 p x (N t) 〉 〉 〉 m x 2 (m x -m) → 0. (3.8)
Proof. Let us start with (3.7). Let A N be the quantity in the left hand side of (3.7), let m x = m xm, and let

ϕ(x) = f (x/N ) 〈 〈 〈 p x (N t) 〉 〉 〉 m x (3.9)
where, for simplicity, we do not write explicitly the dependence of ϕ on N and t. Let 0 < τ < 1. Let

Γ τ N = {1 + k⌊N τ ⌋ : k ∈ } ∩ [1, N ] and, given x ∈ Γ τ N , let x ′ = min{ y ∈ Γ τ N ∪ {N + 1} : y > x}. We decompose A N as A N = 1 N 1-τ x∈Γ τ N 1 N τ x≤z<x ′ ϕ(z) mz = 1 N 1-τ x∈Γ τ N ϕ(x) N τ x≤z<x ′ mz + 1 N 1-τ x∈Γ τ N 1 N τ x≤z<x ′ (ϕ(z) -ϕ(x)) mz =:A (1) N + A (2)

N . To deal with A

(1) N , we observe that ϕ is bounded, see (3.6), so that by Jensen's inequality, E((A (1)

N ) 4 ) ≤ C N 1-τ x∈Γ τ N E 1 N τ x≤z<x ′ m z 4 ≤ C N 2τ
(since the masses m x are i.i.d.). Taking τ > 1/2, this shows that A

(1)

N → 0 almost surely by Borel-Cantelli's lemma. To deal with A (2) N , we start from the definition (3.9) of ϕ and we bound

|ϕ(z) -ϕ(x)| ≤ C m -1 z 〈 〈 〈 p z (N t) 〉 〉 〉 -m -1 x 〈 〈 〈 p x (N t) 〉 〉 〉 + C| f (z/N ) -f (x/N )|. (3.10)
The fact that A

(2)

N → 0 (deterministically) follows then from the bound (3.5) and from the uniform continuity of f . The proof of (3.8) is entirely analogous, with now ϕ(x) = f (x/N )(〈 〈 〈 p x (N t) 〉 〉 〉/m x ) 2 instead of (3.9). By (3.6), this function is bounded and (3.10) is still satisfied since

(m -1 z 〈 〈 〈 p z (N t) 〉 〉 〉) 2 -(m -1 x 〈 〈 〈 p x (N t) 〉 〉 〉) 2 ≤ m -1 z 〈 〈 〈 p z (N t) 〉 〉 〉 + m -1 x 〈 〈 〈 p x (N t) 〉 〉 〉 × m -1 z 〈 〈 〈 p z (N t) 〉 〉 〉 -m -1 x 〈 〈 〈 p x (N t) 〉 〉 〉 ≤ C m -1 z 〈 〈 〈 p z (N t) 〉 〉 〉 -m -1 x 〈 〈 〈 p x (N t) 〉 〉 〉 .
3.3. Proof of the convergence to the linear wave equation (2.8-2.9). For any smooth functions f , g : [0, 1] ∈ such that f (0) = f (1) = 0, the limiting fields and defined in (2.4) can be equivalently characterized as follows:

( f , t) = ( f , 0) - 1 m t 0 ( f ′ , s)ds, (3.11) (g, t) = (g, 0) - t 0 (g ′ , s)ds, (3.12)
and

( f , 0) = 1 0 f (x)r(x)dx, (g, 0) = 1 0 g(x)p(x)dx. ( 3.13)
Let us use this characterization to show that N ( f , t) → ( f , t) and N (g, t) → (g, t). The convergence at t = 0 follows from the strong law of large numbers: N ( f , 0) and N (g, 0) converge almost surely to ( f , 0) and (g, 0) given by (3.13).

Let us next consider t ≥ 0, and let us first deal with N . Integrating the equations of motion yields

N ( f , t) = N ( f , 0) + N t 0 1 N N x=1 f (x/N ) ∇ + M -1 p x (s) ds = N ( f , 0) - N t 0 1 N N x=1 ∇ -f (x/N )m -1 x 〈 〈 〈 p x (s) 〉 〉 〉 ds
where we used the boundary condition f (0) = f (1) = 0 to perform the integration by part. Since

∇ -f (x/N ) = N -1 f ′ (x/N ) + (N -2 ), we obtain N ( f , t) = N ( f , 0) - t 0 1 N N x=1 f ′ (x/N )m -1 x 〈 〈 〈 p x (N s) 〉 〉 〉ds + 1 N .
Using (3.7) in Lemma 2, as well as the dominated convergence theorem to deal with the time integral, we may replace m -1

x by (m) -1 up to an error that vanishes almost surely in the limit N → ∞. Thus

N ( f , t) = N ( f , 0) - 1 m t 0 N ( f ′ , s)ds + ǫ N , (3.14)
where ǫ N → 0 almost surely as N → ∞. Let us next deal with N . This case is simpler since no homogenization over the masses is needed. Proceeding similarly, we find

N (g, t) = N (g, 0) + N t 0 1 N N x=1 g(x/N ) ∇ -r x (s) ds = N (g, 0) - N t 0 1 N N x=1 ∇ + g(x/N )〈 〈 〈 r x (s) 〉 〉 〉 ds = N (g, 0) - t 0 N (g ′ , s)ds + ǫN (3.15)
where we used the boundary condition r 0 (s) = r N (s) = 0 for all time s ≥ 0 to perform the integration by part, and where ǫN → 0 deterministically as N → ∞.

The families N ( f , •) N and N (g, •) N are equicontinuous since a uniform bound on the time derivative of N ( f , •) and N (g, •) holds. Hence, the relations (3.14) and (3.15) implies that any limiting point must satisfy (3.11-3.12).

Pointwise convergence. Thanks to the Hölder regularity of both

〈 〈 〈 r x (N t) 〉 〉 〉 and m -1 x 〈 〈 〈 p x (N t) 〉 〉 〉 expressed by (3.5)
, we deduce a stronger result:

Proposition 1. Let y ∈]0, 1[ and let t ∈]0, 1[. As N → ∞, almost surely (w.r.t. P), 〈 〈 〈 r [N y] (N t) 〉 〉 〉 → r( y, t), 〈 〈 〈 p [N y] (N t) 〉 〉 〉 m [N y] → p( y, t) m . Proof. Let us first deal with 〈 〈 〈 r [N y] (N t) 〉 〉 〉. Let (ρ ε ) ε>0 be a regularizing family: ρ ε ∈ ∞ ( ), supp(ρ ε ) = [-ε, ε], ρ ε ≥ 0 and ρ ε ( y) d y = 1. For y ∈]ε, 1 -ε[, we decom- pose 〈 〈 〈 r [N y] (N t) 〉 〉 〉 = ρ ε y -y ′ 〈 〈 〈 r [N y] (N t) 〉 〉 〉d y ′ = ρ ε y -y ′ 〈 〈 〈 r [N y ′ ] (N t) 〉 〉 〉d y ′ + ρ ε y -y ′ 〈 〈 〈 r [N y] (N t) 〉 〉 〉 -〈 〈 〈 r [N y ′ ] (N t) 〉 〉 〉 d y ′ .
By (3.5), the second term is bounded in absolute value by

ρ ε y -y ′ 〈 〈 〈 r [N y] (N t) 〉 〉 〉 -〈 〈 〈 r [N y ′ ] (N t) 〉 〉 〉 d y ′ ≤ C ε, (3.16)
while the first term is approximated uniformly in N by 1

N N x=1 ρ ε x N -y 〈 〈 〈 r x (N t) 〉 〉 〉.
Thus, by the result shown in Section 3.3, this term converges to

1 0 ρ ε ( y -y ′ )r( y ′ , t)d y ′ (3.17)
as N → ∞. Letting next ε → 0, the continuity of r(•, t) implies that (3.17) converges to r( y, t) while (3.16) converges to 0 as N → ∞. To deal with m -1

[N y] 〈 〈 〈 p [N y] (N t) 〉 〉 〉
, we proceed similarly, using (3.7) in Lemma 2, to get the analog of (3.17).

Finally, thanks to the bound (3.6) and the pointwise convergence result in Proposition 1, and thanks to using (3.7) in Lemma 2 for the field N , we derive (2.8-2.9) by applying the dominated convergence theorem.

EVOLUTION OF THE ENERGY N

In this section we show the limit (2.10). We will assume that f ∈ 1 ([0, 1]). We can then recover the result (2.10) for f ∈ 0 ([0, 1]) by density, and using the a priori estimate x 〈 〈 〈 e x (t) 〉 〉 〉 ≤ CN at all time t ≥ 0.

4.1. Main decomposition of the energy. In order to derive the limit of N , we separate the contribution to the total energy from the temperature (that does not evolve with time) and from mechanical energy, i.e. the average kinetic and potential energy (that does evolve due to the transport of momentum and displacement).

At the macroscopic level, we deduce from (2.5-2.7) that e( y, t) = e( y, 0) + 1

m t 0 ∂ y (r( y, s)p( y, s))ds = 1 β ( y) + p 2 ( y, 0) 2m + r 2 ( y, 0) 2 + 1 m t 0 ∂ s p 2 ( y, s)/2 + mr 2 ( y, s)/2 ds = p 2 ( y, t) 2m + r 2 ( y, t) 2 + 1 β ( y) . ( 4.1)
At the microscopic level, we decompose

N ( f , t) = 1 N N x=1 f (x/N ) 〈 〈 〈 p 2 x 〉 〉 〉 2m x + 〈 〈 〈 r 2 x 〉 〉 〉 2 (N t) = 1 N N x=1 f (x/N ) 〈 〈 〈 p x 〉 〉 〉 2 2m x + 〈 〈 〈 r x 〉 〉 〉 2 2 (N t) + 1 N N x=1 f (x/N ) 〈 〈 〈 p 2 x 〉 〉 〉 2m x + 〈 〈 〈 r 2 x 〉 〉 〉 2 (N t) =: N (t) + N (t), (4.2)
with

p x = p x -〈 〈 〈 p x 〉 〉 〉, r x = r x -〈 〈 〈 r x 〉 〉 〉,
and where and stands respectively for "average" and "fluctuations". Comparing (4.1) and (4.2), we conclude that it is enough to show that, P almost surely, as N → ∞,

N (t) → 1 0 f ( y) p 2 ( y, t) 2m + r 2 ( y, t) 2 d y, ( 4.3
)

N (t) -N (0) → 0, N (0) → 1 0 f ( y) β ( y) d y. ( 4.4)
The limit (4.3) is deduced in the same ways as (2.8-2.9): We first use

(3.8) in Lemma 2 to replace 〈 〈 〈 p x 〉 〉 〉 2 /2m x by m 2 (〈 〈 〈 p x 〉 〉 〉/m x ) 2
up to an error that vanishes as N → ∞. Next, thanks to the bound (3.6) and the pointwise convergence result in Proposition 1, we derive (4.3) by applying the dominated convergence theorem.

The limit (4.4) express the fact that the profile of thermal energy (temperature) remains frozen in time. It will be established thanks to the localization of the high modes of the chain; this is the only place where localization is used. Moreover, we will show in Section 4.3 that in thermal equilibrium, the equality N (t) = N (0) holds without any assumption on the distribution of the masses (besides positivity). This shows thus that Theorem 1 holds actually also for a clean chain if β is constant. Remark 3. The proof of (4.4) holds for any larger time scale N α t, α ≥ 1, i.e.

N (N α-1 t) -N (0) → 0. ( 4.5)
4.2. Convergence of N (t). To deal with the limit (4.4), we will use the fact that any mode of the chain at positive energy is spatially localized in the thermodynamic limit. Hence, we will expand the solutions to the equations of motion into the eigenmodes of the chain. In Section 5 below, we carry this expansion in details and we deduce the needed localization estimates. For our present purposes, it suffices to know the following: There exists a basis {ψ k } 0≤k≤N -1 of N , the basis of the eigenmodes of the chain, so that the solutions to the equations of motion read

r x (t) = N -1 k=1 〈∇ + ψ k , r(0)〉 ω k cos ω k t + 〈ψ k , p(0)〉 sin ω k t (∇ + ψ k ) x ω k , (4.6
)

p x (t) = N -1 k=0 〈ψ k , p(0)〉 cos ω k t - 〈∇ + ψ k , r(0)〉 ω k sin ω k t (M ψ k ) x , (4.7)
where ω 0 = 0 and ω k > 0 for 1 ≤ k ≤ N are the corresponding eigenfrequencies of the chain (we assume that the ω k are sorted by increasing order). Observe that the first term starts from k = 1 while the second one starts from k = 0. Moreover, the orthogonality relation

〈ψ k , M ψ j 〉 = δ(k -j) holds and {ω -1 k ∇ + ψ k } 1≤k≤N -1 forms an orthonormal basis of ( N -1 , 〈•, •〉 N -1
). See Section 5.1 for more details. This representation is useful to exploit localization: all modes with k N are spatially localized. See Section 5.3 for more quantitative estimates. However, low modes with k N remain extended, and we will have to show that the contribution of these modes vanish since their proportion N 1/2 /N → 0 in the thermodynamic limit. Below, for technical reasons, we will replace 1/2 by 1α, for some α > 0 that we will need to choose small enough.

Let 0 < α ≪ 1, let

F 1 = ∩ [0, N 1-α ], F 2 = ∩]N 1-α , N -1],
and let us decompose r(t) = r (1) (t) + r (2) (t) and p(t) = p (1) (t) + p (2) (t) with

r (i) (t) = k∈F i \{0} (. . . ), p (i) (t) = k∈F i (. . . ),
for i = 1, 2 and with (...) the summand featuring in (4.6) or (4.7). We insert this decomposition in N :

N (t) = 1 N N x=1 f (x/N ) p (1)
x + p (2) x 2 2m x + r (1) x + r (2) x 2

(N t).

Let us show the two following limits:

(1)

N (t) = 1 N N x=1 f (x/N ) p (1) x 2 2m x + r (1) x 2 2 (N t) → 0, (4.8)
(2)

N (t) = 1 N N x=1 f (x/N ) p (2) x 2 2m x + r (2) x 2 2 (N t) → f ( y) β ( y) d y, (4.
N (t)| ≤ C 2N N x=1 1 m x k∈F 1 〈ψ k , p(0)〉 cos(ω k N t) - 〈∇ + ψ k , r(0)〉 ω k sin(ω k N t) m x ψ k x 2 + C 2N N x=1 k∈F 1 \{0} 〈∇ + ψ k , r(0)〉 ω k cos(ω k N t) + 〈ψ k , p(0)〉 sin(ω k N t) (∇ + ψ k ) x ω k 2 .
In both terms, one may expand the square so as to get a double sum over k, j ∈ F 1 or k, j ∈ F 1 \{0}, and insert the sum over x inside the sum over k, j. This yields

N x=1 m 2 x m x ψ j x ψ k x = 〈ψ j , M ψ k 〉 = δ(k -j), N x=1 (∇ + ψ j ) x (∇ + ψ k ) x ω j ω k = δ(k -j). Thus | (1)
N (t)| ≤ C 2N k∈F 1 〈ψ k , p(0)〉 cos(ω k N t) - 〈∇ + ψ k , r(0)〉 ω k sin(ω k N t) 2 + C 2N k∈F 1 \{0} 〈∇ + ψ k , r(0)〉 ω k cos(ω k N t) + 〈ψ k , p(0)〉 sin(ω k N t) 2 .
At this point, it suffices to show that there exists a constant C such that, for all k ∈ F 1 ,

〈 〈 〈 〈ψ k , p(0)〉 2 〉 〉 〉 ≤ C, 〈 〈 〈 〈∇ + ψ k , r(0)〉 2 〉 〉 〉 ω 2 k ≤ C, (4.10)
since, bounding sin and cos by 1, and using Cauchy-Schwarz, we obtain

| (1) N (t)| ≤ C N k∈F 1 1 = CN 1-α N → 0.
Let us deal with 〈 〈 〈 〈ψ k , p(0)〉 2 〉 〉 〉 (the other case is analogous):

〈 〈 〈 〈ψ k , p(0)〉 2 〉 〉 〉 = x ψ k x p x (0) 2 = x, y ψ k x ψ k y p x (0) p y (0) = x (ψ k x ) 2 〈 〈 〈 ( p x (0)) 2 〉 〉 〉
where we have used the fact that 〈

〈 〈 • 〉 〉 〉 is a product measure and that 〈 〈 〈 p x (0) 〉 〉 〉 = 0 for all x ∈ {1, . . . , N }. We compute 〈 〈 〈 ( p x (0)) 2 〉 〉 〉 = m x β(x/N
) . Since β is positive and continuous on [0, 1], there exists β -> 0 such that β (x/N ) ≥ β -for all x ∈ {1, . . . , N }. Hence

〈 〈 〈 〈ψ k , p(0)〉 2 〉 〉 〉 ≤ 1 β - N x=1 m x (ψ k x ) 2 = 1 β - 〈ψ k , M ψ k 〉 = 1 β - . ( 4
N (t) - (2)
N (0) → 0 as N → ∞. Let us write

(2) N (t) as a scalar product and expand it in the eigenmodes basis:

(2)

N (t) = 1 2N 〈( f • p (2) )(N t), M -1 p (2) (N t)〉 + 〈( f • r (2) )(N t), r (2) (N t)〉 = 1 2N k∈F 2 〈( f • p (2) )(N t), ψ k 〉〈ψ k , p(N t)〉 + 1 ω 2 k 〈( f • r (2) )(N t), ∇ + ψ k 〉〈∇ + ψ k , r(N t)〉 .
Here g • h denotes a function on ∩ [1, N ] obtained by the usual multiplication in real space between a function g on [0, 1] and

h on ∩ [1, N ], i.e. (g • h) x = g(x/N )h x .
By Lemma 3 stated in Section 5 below, one may associate a localization center x 0 (k) to each mode ψ k with k ∈ F 2 : x 0 (k) is the center of the interval J (k) featuring there (assuming that α is small enough so that the hypotheses of Lemma 3 are satisfied). For each k ∈ F 2 , let us decompose f as

f = f 0 (k) + f k with f 0 (k) = f (x 0 (k)/N )
(thus f 0 (k) is a constant and f k vanishes at x 0 (k)/N ). We insert this decomposition in the above expression for (2) .13) Each expression between 〈 〈 〈 . . . 〉 〉 〉 in the sum in (4.12) represents the energy of the mode ψ k and does not evolve with time, see (5.4) in Section 5 below. Hence, to show

N (t) = 1 2N k∈F 2 f 0 (k) 〈 p(N t), ψ k 〉 2 + 〈 r(N t), ∇ + ψ k 〉 2 ω 2 k (4.12) + 1 2N k∈F 2 〈( f k • p (2) )(N t), ψ k 〉〈ψ k , p(N t)〉 + 1 ω 2 k 〈( f k • r (2) )(N t), ∇ + ψ k 〉〈∇ + ψ k , r(N t)〉 . ( 4
N (t) - (2)
N (0) → 0, we only need to show that the sum in (4.13) converges to 0 as N → ∞.

Let us consider a single term in the sum (4.13), and let us focus on the term involving p (the one involving r is treated the same way). First, by Cauchy-Schwarz,

〈( f k • p (2) )(N t), ψ k 〉〈ψ k , p(N t)〉 ≤ 〈 〈 〈 〈( f k • p (2) )(N t), ψ k 〉 2 〉 〉 〉 1/2 〈 〈 〈 〈ψ k , p(N t)〉 2 〉 〉 〉 1/2 . ( 4
〈 〈 〈 〈ψ k , p(N t)〉 2 〉 〉 〉 = 〈ψ k , p(0)〉 cos ω k N t - 〈∇ + ψ k , r(0)〉 ω k sin ω k N t 2 ≤ 2 〈ψ k , p(0)〉 2 + 〈∇ + ψ k , r(0)〉 2 ω 2 k ≤ C, ( 4
〈 〈 〈 〈( f k • p (2) )(N t), ψ k 〉 2 〉 〉 〉 = x f k (x/N ) p (2) x (N t)ψ k x 2 ≤ x f 2 k (x/N )(ψ k x ) 2 x p (2) x 2 (N t) . ( 4

.16)

For f k , we have the bound

| f k (x/N )| = | f k (x/N ) -f k (x 0 (k)/N )| ≤ C |x -x 0 (k)|
N (this is the only place where we use f ∈ 1 ([0, 1])). Hence, form Lemma 3 below, we deduce that for any ε > 0, the first factor in the right hand side of (4.16) can be bounded by 1/N 2-ε by taking α > 0 small enough. From the conservation of energy (see (5.4) and the bound (4.10)) , the second factor in (4. [START_REF] Komorowski | Ballistic and superdiffusive scales in the macroscopic evolution of a chain of oscillators[END_REF]) is (N ). Hence, for α > 0 small enough, (4.16) goes to zero as N → ∞. Combining this with (4.15), we find that (4.14) goes to zero as N → ∞, and hence that (4.13) converges to 0 as N → ∞. 4.3. Thermal equilibrium case. Assume here that there exists some β > 0 such that β ( y) = β for all y ∈ [0, 1]. Then, we may relax the assumptions on the masses: requiring only that they are all strictly positive, let us show that N (t) = F N (0) for all t ≥ 0. This results from an exact computation.

Since f is arbitrary, it is necessary and sufficient to prove that, for any x,

d dt 〈 〈 〈 p2 x (t) 〉 〉 〉 2m x + 〈 〈 〈 r2 x (t) 〉 〉 〉 2 = 0. We compute p2 x (t) 2m x = 1 2 j,k 〈ψ j , p(0)〉 cos ω j t - 〈∇ + ψ j , r(0)〉 ω j sin ω j t × 〈ψ k , p(0)〉 cos ω k t - 〈∇ + ψ k , r(0)〉 ω j sin ω k t m x ψ j x ψ k x .
A similar expression holds for r2 (t)/2. In order to obtain the expectation with respect to 〈 〈 〈 • 〉 〉 〉, we compute

〈ψ k , p(0)〉〈ψ j , p(0)〉 = δ(k -j) β , (4.17)
〈ψ k , p(0)〉〈∇ + ψ j , r(0)〉 = 0, (4.18)

〈∇ + ψ k , r(0)〉〈∇ + ψ j , r(0)〉 = ω 2 k δ(k -j) β . ( 4.19)
These three properties result from the fact the product structure of ρ loc , from the fact that the variables p and r are centered, and from the the fact that β is constant for (4.17) and (4.19). E.g. (4.17) is obtained by

〈ψ k , p(0)〉〈ψ j , p(0)〉 = x, y ψ k x ψ j y 〈 〈 〈 px (0)p y (0) 〉 〉 〉 = 1 β x m x ψ k x ψ j x = δ(k -j) β .
Hence we have that

〈 〈 〈 p2 x (t) 〉 〉 〉 2m x = 1 2β k (cos 2 ω k t + sin 2 ω k t)m x (ψ k x ) 2 = 1 2β and similarly 〈 〈 〈 r2 x (t) 〉 〉 〉/2 = 1 2β .

EIGENMODES EXPANSION: INTEGRABILITY, LOCALIZATION

We describe an explicit solution to the equations of motion (2.1) in terms of the eigenmodes of the system. This representation is useful to establish the integrability of the system and to exploit the localization at all energies above the ground states (in the thermodynamic limit).

5.1. Solution to the equations of motion. From (2.1), one can deduce second order equations for r and p separately:

rx = ∇ + M -1 ∇ -r x (1 ≤ x ≤ N -1), px = ∆M -1 p x (1 ≤ x ≤ N ),
where, besides the boundary conditions r 0 = r N = 0, we have assumed free boundary conditions for M -1 p, i.e. m -1

0 p 0 = m -1 1 p 1 and m -1 N +1 p N +1 = m -1 N p N .
Notice that there are two different vector spaces: a (N -1)-dimensional space for r with fixed boundary conditions, and a N -dimensional space for M -1 p with free boundary conditions. Moreover, we observe that ∇ + = -(∇ -) † with fixed boundary conditions, and that ∆ = ∆ † with free boundary conditions.

In order to solve the equations of motion, we need to diagonalize two matrices:

∇ + M -1 ∇ - † = ∇ + M -1 ∇ -(of size N -1) and ∆M -1 † = M -1 ∆ (of size N ). Let
us start with the latter: This matrix is not symmetric but the matrix M -1/2 (-∆)M -1/2 is symmetric and non-negative. It admits thus an orthonormal set of eigenvectors, {ϕ k } 0≤k≤N -1 and corresponding eigenvalues ω 2 k that we assume to be sorted by increasing order. Moreover, the spectrum is P-almost surely non-degenerate (see e.g.

Proposition II.1 in [START_REF] Kunz | Sur le spectre des opérateurs aux différences finies aléatoires[END_REF], considering here a perturbation around the non-degenerate equal masses case). Therefore the vectors ψ k = M -1/2 ϕ k are such that M -1 (-∆)ψ k = ω 2 k ψ k , 〈ψ j , M ψ k 〉 = δ( jk).

(5.1)

Because of free boundary conditions, ω 2 0 = 0, and one may chose ψ 0 to be given by

ψ 0 = x m x -1/2
(1, . . . , 1) † .

Next, the matrix -∇ + M -1 ∇ -is symmetric and non-negative, and we denote the eigenvectors by | ψ k 〉 for 1 ≤ k ≤ N -1. It is readily checked that they may be chosen to be given by

ψ k = 1 ω k ∇ + ψ k
with the corresponding eigenvalue given by ω k for 1 ≤ k ≤ N -1. We observe that, by the free boundary conditions on ψ k , ψ k (0) = ψ k (N ) = 0. Given initial conditions r(0), p(0), we can write an explicit solution for r(t), p(t):

〈 ψ k , r〉 = -ω 2 k 〈 ψ k , r〉 (1 ≤ k ≤ N -1), 〈ψ k , p〉 = -ω 2 k 〈ψ k , p〉 (0 ≤ k ≤ N -1). Thus 〈 ψ k , r(t)〉 = 〈 ψ k , r(0)〉 cos ω k t + 〈 ψ k , ∇ + M -1 p(0)〉 ω k sin ω k t (1 ≤ k ≤ N -1),
〈ψ k , p(t)〉 = 〈ψ k , p(0)〉 cos ω k t + 〈∇ -r(0),

ψ k 〉 ω k sin ω k t (0 ≤ k ≤ N -1)
with the convention sin 0 0 = 1 in the second expression at k = 0 (notice that 〈∇ -r(0), ψ k 〉 = -〈r(0), ∇ + ψ k 〉 = 0 for k = 0). This yields therefore r(t) = Therefore, there exists almost surely N 0 so that for for all N ≥ N 0 , the event featuring in (5.5) does not occur. Hence in this case, for all k ∈ I (α) and all |x -y| > N γ , we must have either |ψ k x | ≤ 1/N 1/γ or |ψ k y | ≤ 1/N 1/γ . This means thus that for any k ∈ I (α) there exists an interval J (k) with |J (k)| ≤ 2N γ such that |ψ k x | ≤ N -1/γ for all x / ∈ J (k) ∩ .

  9) which, by Cauchy-Schwarz inequality, implies (4.4). Let us show (4.8). Let us bound | f (x/N )| ≤ C, and use the explicit solution (4

. 11 )

 11 Let us now show (4.9). A computation using the initial measure shows that (0) → f ( y) β( y) d y as N → ∞. Hence, thanks to (4.8), it holds that

( 2 )

 2 N (0) → f ( y) β( y) d y as N → ∞. Thus it suffices to show that

N - 1 k=1 1 k=0. 3 ) 5 . 2 . 2 N - 1 k=0 2 N - 1 k=1 2 N - 1 k=0 2 k 2 ω 2 k= 2 N - 1 k=1 ω 2 k 2 ω 2 k. 5 . 3 .Lemma 3 .

 1135221212122221222533 〈∇ + ψ k , r(0)〉 ω k cos ω k t + 〈ψ k , p(0)〉 sin ω k t ∇ + ψ k 〈ψ k , p(0)〉 cos ω k t -〈∇ + ψ k , r(0)〉 ω k sin ω k t M ψ k .(5Full set of invariant quantities. We observe that the dynamics has N invariant quantities, corresponding to the energy of each mode. It is thus an integrable system. Indeed, let us write the full energy asH = 1 2 〈p, M -1 p〉 + 〈r, r〉 = 1 〈p, ψ k 〉〈M ψ k , M -1 p〉 + 1 〈r, ψ k 〉〈 ψ k , r〉 = 1 〈p, ψ k 〉 2 + 〈r, ∇ + ψ k 〉 2 ωwith the convention that the second term in the last sum is 0 at k = 0. From the evolution equation of the dynamics, one gets that actuallyd dt 〈p, ψ k 〉 2 + 〈r, ∇ + ψ k 〉 0 for all 0 ≤ k ≤ N -1.(5.4)Moreover, by taking specific linear combinations of these conserved quantities, one can obtain conserved quantities that can be written as a sum of local terms. This is for instance the case of the quantity I defined in (3.2), that reads alsoI (r, p) = 1 〈p, ψ k 〉 2 + 〈r, ∇ + ψ k 〉 Localization.Localization can be expressed mathematically in the following strong sense, see[START_REF] Kunz | Sur le spectre des opérateurs aux différences finies aléatoires[END_REF][START_REF] Aizenman | Localization bounds for an electron gas[END_REF][START_REF] Aizenman | Localization at large disorder and at extreme energies: An elementary derivations[END_REF] for the general theory and[START_REF] Verheggen | Transmission Coefficient and Heat Conduction of a Harmonic Chain with Random Masses: Asymptotic Estimates on Products of Random Matrices[END_REF][START_REF] Ajanki | Rigorous Scaling Law for the Heat Current in Disordered Harmonic Chain[END_REF] for precise estimates on the localization length as one approaches the ground state. Let 0 < α <1 2 and letI (α) = ]N (1-α) , N -1] ∩ . There exist constants C, c > 0 such that E k∈I(α) |ψ k x ψ k y | ≤ Ce -c|x-y|/ζ(α) with ζ(γ) = N 2α .We will use this estimate to show that every mode in k ∈ I (α) is supported in a small subset of [1, N ] ∩ up to a small error: Let α, γ > 0 be such that 2α < γ < 1. There exists almost surely N 0 ∈ so that for all N ≥ N 0 , and for all k ∈ I (α), there exists an interval J (k) with |J (k)| ≤ 2N γ such that |ψ k x | ≤ N -1/γ for all x / ∈ J (k) ∩ .Proof. Let us first show thatP := P ∃k ∈ I (α), ∃x, y ∈ [1, N ] ∩ : |x -y| ≥ N γ , |ψ k x | ≥ N -1/γ , |ψ k y | ≥ N -1/γ ≤ C(α, γ)e -N (γ-2α)/2(5.5) Indeed we computeP ≤ k∈I(α) x, y:|x-y|≥N γ P |ψ k x | ≥ N -1/γ , |ψ k y | ≥ N -1/γ ≤ k∈I(α)x, y:|x-y|≥N γ P |ψ k x ψ k y | ≥ N -2/γ ≤ k∈I(α) x, y:|x-y|≥N γ N 2/γ E(|ψ k x ψ k y |) ≤ CN 2/γ x, y:|x-y|≥N γ e -c|x-y|/ζ(α) ≤ CN 2/γ N 2 e -N γ-2α ≤ C(α, γ)e -N (γ-2α)/2 .

  .[START_REF] Jara | Quenched non-equilibrium central limit theorem for a tagged particle in the exclusion process with bond disorder[END_REF] The second factor in (4.14) is bounded by a constant: