Gra Bu

L Miguel Kabil

Rodrigues

SPECTRAL VALIDATION OF THE WHITHAM EQUATIONS FOR PERIODIC WAVES OF LATTICE DYNAMICAL SYSTEMS

Keywords: periodic traveling waves, Whitham averaging, modulation systems, lattice dynamical systems, spectral stability, discrete Bloch transform. 2010 MSC: 34K13, 34K31, 35B10, 35B27, 35B35, 35P05, 37K60, 37L60, 39A23, 39A30

In the present contribution we investigate some features of dynamical lattice systems near periodic traveling waves. First, following the formal averaging method of Whitham, we derive modulation systems expected to drive at main order the time evolution of slowly modulated wavetrains. Then, for waves whose period is commensurable to the lattice, we prove that the formally-derived first-order averaged system must be at least weakly hyperbolic if the background waves are to be spectrally stable, and, when weak hyperbolicity is met, the characteristic velocities of the modulation system provide group velocities of the original system. Historically, for dynamical evolutions obeying partial differential equations, this has been proved, according to increasing level of algebraic complexity, first for systems of reaction-diffusion type, then for generic systems of balance laws, at last for Hamiltonian systems. Here, for their semi-discrete counterparts, we give at once simultaneous proofs for all these cases. Our main analytical tool is the discrete Bloch transform, a discrete analogous to the continuous Bloch transform. Nevertheless, we needed to overcome the absence of genuine space-translation invariance, a key ingredient of continuous analyses.

A common strategy to tackle the analysis of large time behavior of dynamical systems focuses on coherent structures, often playing the role of organizing centers for the time evolution. We Date: April 18, 2015. Research of L.M.R. was partially supported by the ANR project BoND ANR-13-BS01-0009-01. follow here this line of investigation for lattice dynamical systems and restrict our attention to periodic traveling waves.

By definition, evolution of lattice dynamical systems is continuous in time but discrete in space. We consider here particular instances of those where the lattice is one-dimensional and thus may be assumed to be Z, and the evolution obeys a differential equation d dt

U (t) = F(U (t))
where the unknown U associates to any time t ∈ R an R d -valued sequence U (t) ∈ (R d ) Z (d ∈ N * being a fixed dimension) and F is a vector-field on (R d ) Z that preserves the linear space of finitely-supported sequences and acts smoothly on it. Alternatively one may view the evolution as given by an infinite number of differential equations on R d -valued functions coupled by a finite-range interaction. We expect that our results could actually deal with nonlocalized interaction satisfying certain short-range decaying assumptions but we choose to ignore those here to keep technicalities as low as possible.

There are at least two ways in which such systems arise. They may come directly from modeling at a discrete level, as for neural networks, granular media, crystals, biological molecules, optical waveguides, chains of coupled oscillators... Alternatively they may emerge from the discretization in space -often called semi-discretization -of systems of partial differential equations, that generate fully continuous dynamical systems, In the latter case, mark that although a study of dynamical evolutions generated by fully discrete schemes including boundary conditions would be of much more direct practical use, such a study appears as a daunting task and it is a common belief that still relevant pieces of information are obtained from investigation of extended semi-discretized schemes. The reader looking for panoramas on lattice dynamical systems is referred to [START_REF] Chow | Pattern formation and spatial chaos in lattice dynamical systems. I[END_REF][START_REF] Chow | Dynamics of lattice differential equations[END_REF][START_REF] Mallet-Paret | Spatial patterns, spatial chaos, and traveling waves in lattice differential equations[END_REF][START_REF] Pankov | Travelling waves and periodic oscillations in Fermi-Pasta-Ulam lattices[END_REF][START_REF] Panayotis | Non-linear waves in lattices: past, present, future[END_REF]. Specially designed technical tools of wide application range may be found in [Rus89a, Rus89b, CMPS98, MP99a, MP99b, Ioo00, KK01, CL07, Jam03, FP99, FP02, FP04a, FP04b, HVL07, HVL08, HVL09]. Especially the latter references include key-ingredients to prove existence of traveling waves for various class of systems, a fact that is taken as an assumption in the present piece of work.

We are interested here in a particular area of the dynamics, namely what occurs in neighborhoods of periodic traveling waves. Compared to other coherent structures of lattice dynamical systems -such as fronts, kinks, pulses, shocks, solitary waves or breathers -, periodic waves have received relatively less attention. In particular, while the stability analysis of some of the former patterns seems to have reached some maturity -see for instance [MP99b, FP99, FP02, FP04a, FP04b, BGHR03, BHSZ10, HS13] -, the authors of the present paper are not aware of even a single example of a comprehensive stability study for periodic waves1 of some lattice dynamical system. Nevertheless, the interested reader may benefit from looking at [Kee80, TT95, BR97, WZ97, FV99, PP00, IK00, FR02, Pan05, Car05, FR05, Bak07, GLR09, Mak11, VZ11, Her12, Jam12, BP13, LHC14], that are at least partly dedicated to periodic waves of lattice dynamical systems, and mostly focusing on proving their existence.

To set our precise framework, let us specify that we say that, for systems as above, a traveling wave, sometimes called a uniformly sliding solution, -that is, a solution U with special form U (t) = (u(j -ct)) j∈Z for some speed c ∈ R and some profile u : R → R d -is said to be periodic if the corresponding profile u is itself periodic. In this case we rescale phase variable of the profile to ensure that the profile period is one. Explicitly any periodic traveling wave U is written as U (t) = (u(k j + ω t)) j∈Z with a one-periodic profile u. This brings out spatial wavenumber k ∈ R * and time frequency ω = -k c. Mark that up to now we are not assuming that the spatial wavenumber k is rational. When its is not the case the spatial periodicity of the profile is not a priori easily observed on the solution itself, though when standing waves are discarded, that is when restricting to c = 0 this periodicity is apparent in time, the solution being 1/ω-periodic in time. We refer the reader to [START_REF] Herrmann | Oscillatory waves in discrete scalar conservation laws[END_REF] for further elaboration on this comment including fascinating illustrating pictures. This already stresses that, in contrast with analogous issues for continuous dynamical systems, analysis of the behavior near standing waves follows from different kind of arguments. Indeed, while profile equations for genuinely-traveling waves are differential functional equations, profile equations for standing waves are algebraic functional equations. This results in the fact that even when traveling and standing waves coexist the standing-wave limit c → 0 is very singular. For this reason, our analysis focused on on genuinely-traveling waves do not extend in a straightforward way neither to breathers nor to steady spatially periodic solutions.

Our goals are, first, to derive -on a formal basis -a system of partial differential equations that is expected to describe the time evolution of slowly modulated wavetrains of lattice dynamical systems, second, to elucidate -on mathematical grounds -what are, at the spectral level, the connections between the linearization about a constant state of the formally-derived modulation systems and the linearization about a periodic wave of our initial lattice system of equations. Our formal derivation is in many ways very classical and it is well-known as one of the possibilities to recover a continuous description from a discrete model. The reader is referred to [HLM94, HLM95, DHM06, DHR06, GHM06, DH08] for various perspectives on this point of view. Our main motivation here is to provide a rigorous spectral validation, that should be compared with some contributions of [DSSS09, Ser05, NR13, BGNR14] on partial differential equations, essentially as [FP99, FP02, FP04a, FP04b] may be thought as a discrete counterpart of [START_REF] Pego | Eigenvalues, and instabilities of solitary waves[END_REF][START_REF] Pego | Asymptotic stability of solitary waves[END_REF].

As a direct consequence of our results stems the proof that characteristic speeds of averaged systems do provide group velocities for the original one. Mark that even if we illustrate our strategy mostly on systems obtained from semi-discretization of continuous models we provide a proof at the discretized level so that, in these particular cases, group velocities implicitly depend on the mesh size and of the particular choice of discretization. Given the fundamental role of a precise knowledge of group velocities in a quantitative description of the dynamics, and since accurate reproduction of dispersion relations obtained by linearizing about constant states has served for a long time now to discriminate between numerical schemes with, from other respects, similar performances, we do not exclude that in a near future comparison between discrete group velocities and continuous ones could serve similar purposes.

Our contribution may seem relatively modest in that the proved connection is purely at the spectral level. Yet on the other hand our strategy is very robust as we emphasize by treating here in a single place various classes systems that were considered at the continuous level in separate contributions. Moreover, as our proofs provide relations for eigenmodes -including eigenvectors and not restricting to eigenvalues -, we expect that it could be directly useful for a dynamical validation of averaged systems in the large-time limit, as, in continuous cases, [DSSS09, Section 4.2] for [SSSU12, JNRZ13], [START_REF] Noble | Whitham's modulation equations and stability of periodic wave solutions of the Korteweg-de Vries-Kuramoto-Sivashinsky equation[END_REF] for [START_REF] Johnson | Behavior of periodic solutions of viscous conservation laws under localized and nonlocalized perturbations[END_REF] or [START_REF] Benzoni-Gavage | Slow modulations of periodic waves in Hamiltonian PDEs, with application to capillary fluids[END_REF] for [Rod].

From our point of view, our main restriction is that in the spectral validation we only consider periodic waves whose period is a multiple of the generator of the lattice. This stringent restriction does not fit well in the spirit of slowly modulated wavetrains that involve a continuum of such periodic waves but it enables us to use techniques originating in the study of differential operators with periodic coefficients -the classical Floquet theory -rather than those for quasi-periodic coefficients. For a glimpse at the fascinating but technically daunting quasi-periodic theory the reader is referred to [START_REF] Håkan Eliasson | Discrete one-dimensional quasi-periodic Schrödinger operators with pure point spectrum[END_REF][START_REF] Håkan Eliasson | Perturbations of linear quasi-periodic system[END_REF][START_REF] Håkan Eliasson | One-dimensional quasi-periodic Schrödinger operators-dynamical systems and spectral theory[END_REF][START_REF] Håkan Eliasson | On the discrete one-dimensional quasi-periodic Schrödinger equation and other smooth quasi-periodic skew products[END_REF] where is analyzed the spectrum of operators similar to those obtained by linearizing discrete equations about a general spatially-periodic stationary solution.

From a technical point of view, in our analysis of lattice systems, departures from the continuous strategy originate in the loss of a genuine space-translation invariance of the original system -see Remark 3.1 though. To balance this, we make an extensive use of the remaining time-translation invariance and of the space-translation invariance of the profile equation. In particular, instead of moving to a mobile frame in which the wave is at rest and analyzing the spectrum of the generator of the linearized dynamics, we study directly Floquet multipliers of the linearized time evolution, that is, we study the spectrum of the map that encodes the evolution over a time 1/ω according to the linearized dynamics. Likewise, when dealing with Hamiltonian systems, we replace in our arguments the missing momentum conservation law with local conservation of energy. On the other hand, we may still use a normal, Floquet exponent by Floquet exponent, Bloch-wave decomposition by relying upon the discrete Bloch transform. Then to examine the spectrum of resulting discrete Bloch symbols we perform here a direct spectral perturbation analysis "à la Kato" [START_REF] Kato | Perturbation theory for linear operators[END_REF] as in [START_REF] Noble | Whitham's modulation equations and stability of periodic wave solutions of the Korteweg-de Vries-Kuramoto-Sivashinsky equation[END_REF][START_REF] Benzoni-Gavage | Slow modulations of periodic waves in Hamiltonian PDEs, with application to capillary fluids[END_REF] for continuous systems. Alternatively we could also have introduced a suitable Evans function, as done in [START_REF] Serre | Spectral stability of periodic solutions of viscous conservation laws: large wavelength analysis[END_REF] for continuous dynamics and in [KK01, [START_REF] Benzoni-Gavage | Stability of semi-discrete shock profiles by means of an Evans function in infinite dimensions[END_REF][START_REF] Cramer | Fredholm determinants and the Evans function for difference equations[END_REF] for other classes of coherent structures in lattice systems. For a detailed discussion of respective advantages of proofs by direct spectral perturbation or by Evans function expansions see [START_REF] Rodrigues | Asymptotic stability and modulation of periodic wavetrains, general theory & applications to thin film flows[END_REF].

The rest of the present paper is organized as follows. In the second section, examining formal expansions, we derive modulated equations. In the third section we introduce tools from functional analysis -including the discrete Bloch transform -that we shall subsequently use to analyze the spectrum of the linear one-period evolution operator. In Section 4 we state and prove our main results concerning the spectral validation of the previously obtained averaged systems. We close this paper with further comments and open questions. In Sections 2 & 4, we analyze three classes of systems in increasing order of algebraic complexity. First, we consider a lattice system of reaction-diffusion type, in the sense that the system does not support any local conservation law. In this case, the averaged evolution obeys a scalar equation and both derivation and spectral validation are considerably simpler so that we also provide high-order versions of those. It is actually those higher-order versions that are needed in the dynamical large-time validation of the slow modulation scenario. After that, we examine a general system that includes built-in conservation laws. The averaged dynamics is now prescribed by a genuine system and the spectral validation includes spectral perturbation of a Jordan block, hence requires a preliminary desingularization. At last, we consider a lattice system of Hamiltonian type that comes with an "extra" conservation law encoding conservation of energy. Again, all the corresponding theorems, Theorem 4.1, Theorem 4.2 and Theorem 4.3, are stated and proved in Section 4.

Notation. Elements U ∈ (C d ) Z of (C d ) Z are implicitly assumed to have coordinates U j = U (j) ∈ C d , j ∈ Z. Likewise for U ∈ (C d ) Z×R , U j (t) = U (j, t), for (j, t) ∈ Z × R, and U (t) = U (•, t) ∈ (C d ) Z , for t ∈ R. All along T denotes the left shift operator on (C d ) Z , (T U ) (j) = U (j + 1) , for U ∈ (C d ) Z .
The set N is the set of nonnegative integers, and, for (M,

N ) ∈ Z, [[M, N ]] = [M, N ] ∩ Z.
The resolvent set and the spectrum set of an operator are denoted respectively by ρ(•) and σ(•).

Formal derivation of averaged equations

In this section we derive, from formal considerations, averaged modulation systems for three classes of nonlinear lattice systems, starting with the simplest case of lattice systems of reactiondiffusion type. To do so, we will follow the two-timing method introduced by Whitham [START_REF] Whitham | Linear and nonlinear waves[END_REF]. Only afterwards, in Section 4 , shall we prove that the original dynamical evolution is indeed related to these modulation equations.

Reaction-diffusion case. We consider a discretization of the scalar reaction-diffusion equation for u

: R × R + → R d ∂ t u = ∆u + f (u) (2.1)
where f : R d → R d is a given (in general nonlinear) smooth function. Using a centered difference scheme to discretize the spatial Laplacian operator, we obtain the following system of coupled ordinary differential equations

d dt U j (t) = µ[U j+1 (t) -2U j (t) + U j-1 (t)] + f (U j (t)), j ∈ Z, t > 0, (2.2)
where µ ∈ R is a given fixed constant. To consider it in the abstract form of dynamical lattice systems, we introduce F the nonlinear operator on U ∈ (C d ) Z associated with f and given by F (U ) = (f (U j )) j∈Z . Then, we write system (2.2) for

U (•) = (U j (•)) j∈Z as d dt U (t) = µ (T -2Id + T -1 ) U (t) + F (U (t))), t > 0 . (2.3)
A profile u generates through U (t) = (u(kj + ω t)) j∈Z , a periodic traveling wave to (2.3) with spatial wavelength k = 0 and time frequency ω if and only if

ωu (ζ) = µ[ u(ζ + k) -2 u(ζ) + u(ζ -k) ] + f (u(ζ)), ζ ∈ R , (2.4) and u is one-periodic (2.5) u(ζ + 1) = u(ζ), ζ ∈ R .
To obtain (even formal) pieces of information about a given periodic wave, we need a precise knowledge of neighboring periodic waves. With this respect we make the following "non degeneracy" assumption:

(A)
Periodic traveling waves of (2.3) -identified when coinciding up to time translation -form a smooth manifold of dimension one, regularly parametrized by k.

For latter use, we denote by u k (•), ω(k) such a parametrization by k. We refer the reader to references given in the introduction for proofs of existence results. We also stress that our assumption is consistent with both the continuous limit leading to (2.1) and similar considerations in [HLM94, HLM95, DHM06, DHR06, GHM06, DH08].

Our formal ansatz is of two-scale type, one scale slow -but otherwise arbitrary -, involving macroscopic variables, another one fast but oscillatory, hence described through an oscillation phase. Explicitly we want to gain some insight about the asymptotic behavior -when the slow frequency goes to zero -of families of solution (U ( ) (•)) expanding as U ( )

j (t) = u ( )     εj X , εt T , φ ( ) ( j, t) Θ     + O( 2 ), (2.6) with phases φ ( ) φ ( ) (X, T ) = φ 0 (X, T ) + φ 1 (X, T ) + 2 φ 2 (X, T ) + O( 3 ) , (2.7) and profiles u ( ) one-periodic in the third variable Θ u ( ) (X, T, Θ) = u 0 (X, T, Θ) + u 1 (X, T, Θ) + O( 2 ) .
This yields the following relation

∂ T u ( ) X, T, φ ( ) (X, T ) + (∂ T φ ( ) (X, T )) • ∂ Θ u ( ) X, T, φ ( ) (X, T ) = µ u ( ) X + , T, φ ( ) (X + , T ) -2u ( ) X, T, φ ( ) (X, T ) + u ( ) X -, T, φ ( ) (X -, T ) + f u ( ) X, T, φ ( ) (X, T ) + O( 2 ) .
Usual identification of powers of and replacement of the phase φ 0 (X, T )/ + φ 1 (X, T ) with Θ then provide at order 0

(∂ T φ 0 (X, T )) ∂ Θ u 0 (X, T, Θ) = µ [ u 0 (X, T, Θ + ∂ X φ 0 (X, T )) -2u 0 (X, T, Θ) + u 0 (X, T, Θ -∂ X φ 0 (X, T )) ] + f (u 0 (X, T, Θ)) ,
which is, in the Θ variable, the functional profile equation with wavenumber and time frequency

k 0 (X, T ) := ∂ X φ 0 (X, T ) and ω 0 (X, T ) := ∂ T φ 0 (X, T ) . (2.8)
Hence, under a suitable normalization of parametrization, without loss of generality (2.9)

u 0 (X, T, Θ) = u k 0 (X,T ) (Θ)
where evolution of the local wavenumber obeys

(2.10)

∂ T k 0 -∂ X (ω(k 0 )) = 0 , stemming from Schwarz' identity ∂ T ∂ X φ 0 = ∂ X ∂ T φ 0 .
As a result of our heuristics, we obtain the conjecture that slow/oscillatory solutions evolve at main order according to a slow modulation scenario U ( )

j (t) = u k 0 ( j, t) φ 0 ( j, t) + φ 1 ( j, t) + O( ),
that is, locally, at scale 1, the solutions look like one member of the periodic traveling wave family, but with local parameters -phase shift φ 1 and wavenumber k 0 -evolving on a slow scale 1/ . In this description local phase and local spatial wavenumber are tied by k 0 = ∂ X φ 0 and the evolution of the local wavenumber k 0 obeys at main order (2.10). In particular, unsurprisingly, in this simple context, linear group velocity around a given wave U generated by a profile u k is expected to be ω (k).

In this simple case where the family of periodic waves is one-dimensional and hence modulation equation (2.10) is scalar, algebra required to raise the order of description remains relatively simple so that we explain it now. It requires a more precise account of assumption (A). In assuming a regular parametrization by k, we include that for any wave under considerationu with wavenumber k and frequency ω -the operator L acting on L 2 per ([0, 1]) with domain

H 1 per ([0, 1]) through (2.11) Lu = -ω u + µ[ u(• + k) -2u + u(• -k) ] + df (u(•)) u
has one-dimensional kernel. By time translation invariance, this kernel is thus C u . Moreover since, as a (relatively) compact perturbation of -ω( ) , L is Fredholm of index 0 -when considered as usual as a bounded operator on its domain endowed with its graph norm -, we conclude the existence of a u ad k ∈ H 1 per ([0, 1]) spanning the kernel of the adjoint L * of L and such that u ad k , u = 1 (in the L 2 per ([0, 1]) sense). In particular by differentiating the profile (2.4) and taking a scalar product with u ad k , one receives (2.12)

ω (k) = u ad k , µ[u (• + k) -u (• -k)] .
Moreover, u ad k allows us to give a precise account of what is the appropriate normalization of parametrization that we have alluded to above. Indeed we explicitly require (N) u ad k , ∂ k u k = 0 in order to enforce the distinct respective roles of wavenumber and phase shift. The above condition may be achieved by shifting the original u k in a suitable k-dependent way.

After these preliminaries we may go back to the identification of terms of order 1 . By dropping all (X, T ) dependencies, we receive

-Lu 1 + ∂ T φ 1 ∂ Θ u 0 -∂ X φ 1 µ [∂ Θ u 0 (• + k 0 ) -∂ Θ u 0 (• -k 0 )] = -∂ T k 0 ∂ k u k 0 + µ∂ X k 0 [∂ k u k 0 (• + k 0 ) -∂ k u k 0 (• -k 0 )] + 1 2 µ∂ X k 0 [∂ Θ u 0 (• + k 0 ) + ∂ Θ u 0 (• -k 0 )]
where L is here an (X, T )-parametrized version of original L, associated with u 0 and acting in the Θ-variable. Averaging against u ad k 0 yields

∂ T φ 1 -ω (k 0 ) k 1 = d(k 0 ) ∂ X k 0 .
where, for concision's sake, we have introduced a diffusion coefficient

(2.13) d(k) = u ad k , µ [∂ k u k (• + k) -∂ k u k (• -k)] + 1 2 µ [∂ ζ u k (• + k) + ∂ ζ u k (• -k)] .
As a refined conclusion, we obtain that slow/oscillatory solutions evolve at main order according to U ( )

j (t) = u k( j, t) φ( j, t) + O( ) ,
with local phase and local spatial wavenumber tied by k = ∂ X φ and the evolution of the local wavenumber k obeys at second order

(2.14) ∂ T k -∂ X (ω(k)) = ∂ X (d(k)∂ X k) .
Mark that in this second-order description, there is no more undetermined local phase shift and that, when present, large-time diffusive decay or anti-diffusive growth may be captured by equation (2.14).

General mixed case.

The main point about discrete reaction-diffusion systems considered in the previous subsection is that they do not support any built-in conservation law. Assumption A then requires that indeed no hidden conservation is present. This leads to a scalar modulation behavior encoded by (2.10) or (2.14). We now relax this stringent assumption to consider a system mixing conservative equations with non conservative ones. At the continuous level such structure is ubiquitous. For instance it naturally emerges from the modeling of the evolution of isentropic compressible flows undergoing external forces by the Navier-Stokes or Euler systems. Typically, in the former case one receives for the evolution of mass density ρ and velocity field u a system in the form (2.15)

∂ t ρ + ∂ x (ρ u) = 0 ∂ t (ρu) + ∂ x (ρu 2 + P (ρ)) = ∂ x (ν(ρ)∂ x u) + g(ρ, u),
where P , ν and g provide respectively the pressure law, Lamé viscosity coefficient and external forces.

Here we consider for

U = (R, W ) T -valued in (R d ) Z = (R d 1 ) Z × (R d 2 ) Z -systems of the following general form (2.16) d dt R(t) + D 1 (f r (R(t), W (t))) = 0 d dt W (t) + D 2 (f w (R(t), W (t))) = D 3 (B(R(t), W (t))D 4 W (t)) + g(R(t), W (t))
,

where f = (f r , f w ) T : R d → R d , g : R d → R d 2 and B : R d → L(R d 2 ) are identified with their component-wise counterparts actions on sequences, of respective type (R d ) Z → (R d ) Z , (R d ) Z → (R d 2 ) Z and (R d ) Z → L((R d 2 ) Z
), and D j , j = 1, 2, 3, 4, are constant-coefficient discrete differential operators, D 1 being conservative in the sense that the kernel of D * 1 (the adjoint of D 1 on 2 (Z; R d 1 )) contains constant sequences. From now on, for definiteness' sake, we will restrict our attention to the case where for some η > 0

(2.17)

D 1 = D 4 = η (T -Id) , D 2 = D 3 = η (Id -T -1 ) .
As in the reaction-diffusion case, we make a "non degeneracy" assumption on the set of periodic traveling waves:

(A)
Periodic traveling waves of (2.16) -identified when coinciding up to time translation -form a smooth manifold of dimension d 1 + 1, regularly parametrized by (k, M), their wavenumber and the average-values of their d 1 th first components.

Explicitly M = (M 1 , • • • , M d 1 ) with M j =
1 0 e j • u with e j the jth element of the canonical basis of R d . For latter use, we denote by u k,M (•), ω(k, M) a parametrization provided by the previous assumption. Again we refer the reader to references in the introduction for proofs of existence results.

Again, a formal ansatz U ( )

j (t) = u ( )     εj X , εt T , φ ( ) ( j, t) Θ     + O( 2 ),
(2.18)

with phases φ ( ) φ ( ) (X, T ) = φ 0 (X, T ) + φ 1 (X, T ) + 2 φ 2 (X, T ) + O( 3 ) , (2.19)
and profiles u ( ) one-periodic in the third variable Θ

u ( ) (X, T, Θ) = u 0 (X, T, Θ) + u 1 (X, T, Θ) + O( 2 )
leads, under a suitable normalization of parametrization, to

u 0 (X, T, Θ) = u (k 0 ,M 0 )(X,T ) (Θ) with k 0 (X, T ) := ∂ X φ 0 (X, T ) , ω 0 (X, T ) := ∂ T φ 0 (X, T )
and

M 0 (X, T ) := 1 0 I d 1 0 d 1 ×d 2 u 0 (X, T, Θ) dΘ .
But now the law of conservation of waves

∂ T k 0 -∂ X (ω(k 0 , M 0 )) = 0
fails to describe completely the time-evolution.

To proceed, one collects terms of power 1 , replacing φ 0 (X, T )/ + φ 1 (X, T ) with Θ, and, with usual implicit notation for coordinates splitting u = (r, w) T , receives

-Lu 1 + ∂ T (u 0 ) + ηT k 0 (df (u 0 ) • ∂ X u 0 ) = ∂ Θ (• • • ) + 0 • • • ,
where L is the operator associated with the linearized profile equation and T k 0 acts on functions of Θ through T k 0 (v) = v(• + k 0 ). Motivated by the fact that constant functions with values e j , j = 1, . . . , d 1 , lie in the kernel of L * , we average over Θ ∈ (0, 1) the scalar product of the previous system with those e j and obtain

∂ T M 0 + ∂ X F (k 0 , M 0 ) = 0
where

F (k, M) = η 1 0 f r (u k,M (Θ)) dΘ .
Altogether we have derived the following averaged system

   ∂ T k -∂ X ω(k, M) = 0 ∂ T M + ∂ X F (k, M) = 0 (2.20)
for the evolution of local parameters involved in a slow modulation description of the dynamics.

2.3. Hamiltonian case. Up to now we have implicitly assumed that our systems do not contain any "hidden" conservation law. In doing so we were motivated by our will to restrict to nondegenerate cases. However for special algebraic classes of systems 'hidden' conservation laws are indeed generic and we need to accommodate them. As a typical example we analyze systems that are (at least formally) of Hamiltonian type. Incidentally, although we do not follow this path here, we stress that the presence of an extra structure also offers alternative ways of deriving the same modulation equations; see [START_REF] Whitham | Linear and nonlinear waves[END_REF].

To emphasize analogy with results on continuous dynamical systems in [START_REF] Benzoni-Gavage | Slow modulations of periodic waves in Hamiltonian PDEs, with application to capillary fluids[END_REF] we focus on a class of Hamiltonian lattice dynamical systems that include discrete counterparts to the Korteweg-de Vries equation and the Euler-Korteweg system. Explicitly we fix some Hamiltonian H : (R d ) 2 → R and some skew-symmetric constant-coefficient discrete differential operator J in a conservative form, say for definiteness' sake J = D B with B a symmetric matrix of size d and D = η 2 T -T -1 . Moreover we choose some other constant-coefficient discrete differential operator D (not necessarily skew-symmetric), say D = η (T -Id), for some η > 0.

Then we consider the following lattice dynamical system

(2.21) d dt U (t) = J δH[U (t)],
where δ denotes a discrete Euler operator providing variational derivatives for functionals stemming from local functions of (U, DU ),

(2.22) δH[U ] = ∇ U H(U, D U ) + D * ∇ D U H(U, DU ) .
By definition the previous system comes with a local conservation law for the discrete Hamiltonian

(2.23) d dt [H(U, D U )] = D 1 2 T -1 (δ H[U ]) B δ H[U ] + T -1 (∇ D U H(U, D U )) • J δ H[U ] .
Notice that to derive (2.23) from (2.21) we have used that

A DB -( D * A) B = D((T -1 A)B) for any A, B ∈ (C d ) Z .
As we have already done in previous subsections, to the shift operator T acting on sequences we may naturally associate a family of shifts on functions T k (v) = v(• + k). For concision's sake, from now on we will likewise denote by (P (T)) k the operator P (T k ). In particular, we will freely use notation D k , D k , δ k ... Now, our "non-degeneracy" assumption on the set of periodic traveling waves takes the following form (A) Periodic traveling waves of (2.21) -identified when coinciding up to time translation -form a smooth manifold of dimension d + 2, regularly parametrized by (k, M, E), their wavenumber, the average-values of their components and their Hamiltonian.

More explicitely

E = 1 0 H(u (k,M,E) (ζ), D k u (k,M,E) (ζ)) dζ .
With notational convention similar to the one of previous subsections, a two-scale ansatz leads, under a suitable normalization of parametrization, to

u 0 (X, T, Θ) = u (k 0 ,M 0 ,E 0 )(X,T ) (Θ) with k 0 (X, T ) := ∂ X φ 0 (X, T ) , ω 0 (X, T ) := ∂ T φ 0 (X, T ) , M 0 (X, T ) := 1 0 u 0 (X, T, Θ) dΘ , and
E 0 (X, T ) := 1 0 H(u 0 (X, T, Θ), D k 0 (X,T ) u 0 (X, T, Θ)) dΘ .
Moreover we still have a law of conservation of waves

∂ T k 0 -∂ X (ω(k 0 , M 0 )) = 0
and a conservation law for the averaged-values

∂ T M 0 = B ∂ X F (k 0 , M 0 , E 0 ),
where

F (k, M, E) = η 1 0 δ k (H)[u k,M,E ](ζ) dζ = η 1 0 ∇ U H(u k,M,E (ζ), D k u k,M,E (ζ)) dζ .
The latter fact is easily deduced from the fact that by collecting terms of power 1 stemming from the insertion of our ansatz in our system we receive an equation in the form

-Lu 1 + ∂ T (u 0 ) -η 2 (T k 0 + T -k 0 )B∂ X (δ k 0 (H)[u 0 ]) = ∂ Θ (• • • ) + D k 0 (• • • ) .
The derivation of an equation for the time evolution of E 0 requires more algebra. To write computations in a compact way we introduce some more pieces of notation. First we set L ξ := e -iξ• L e iξ• then we denote L (1) and L (2) the operators involved in the expansion

L ξ = L + (iξ k) L (1) + (iξ k) 2 L (2) + O(|ξ| 3 ) .
With these preliminaries the equation obtained from the collection of coefficients of 1 is explicitly written as

-Lu 1 + ∂ T (u 0 ) + (∂ T φ 1 + c(k 0 , M 0 , E 0 )∂ X φ 1 ) ∂ Θ u 0 = ∂ X φ 1 L (1) ∂ Θ u 0 + ∂ 2 X φ 0 L (2) ∂ Θ u 0 + (L (1) -c(k 0 , M 0 , E 0 )) ∂ X u 0 . Since δ k 0 [H](u 0 )
belongs to the kernel of L * , we shall average over Θ ∈ (0, 1) the scalar product of the previous system with δ k 0 [H](u 0 ). Our claim is that we receive

∂ T E 0 -∂ X (S(k 0 , M 0 , E 0 )) = 0 where S(k, M, E) = η 1 0 1 2 T -k (δ k (H)[u (k,M,E) ]) B δ k (H)[u (,M,E) ] + η 1 0 T -k (∇ D U H(u (k,M,E) , D k u (k,M,E) )) • J k δ k (H)[u (k,M,E) ]
Indeed it does follow from direct computations. One one hand, when a

∈ {ζ, k, M 1 , • • • , M d , E}, one computes that (2.24) δ k [H](u) • ∂ a u = ∂ a (H(u, D k u)) -ηT -k (∇ D U H(u, D k u)) • ∂ ζ u ∂ a k + (• • • )
where (• • • ) denotes mean-free terms and dependence of u on parameters has been omitted. On the other hand, for a ∈ {ζ, k, M 1 , • • • , M d , E}, we also have that (2.25)

δ k [H](u) • (L (1) -c)∂ a u + L (2) (∂ ζ u)∂ a k = ∂ a η 1 2 T -k (δ k (H)[u]) B δ k (H)[u] + ηT -k (∇ D U H(u, D k u)) • J k δ k (H)[u] -ηT -k (∇ D U H(u, D k u)) • ∂ ζ u ∂ a ω + (• • • )
with the same convention. This yields our claim by integration using that

∂ T k 0 = ∂ X (ω(k 0 , M 0 , E 0 )).
Alternatively one could have derived this last equation by directly expanding and averaging equation (2.23).

Altogether we have derived the following averaged equations for the slow modulation of local parameters (2.26)

     ∂ T k = ∂ X ω(k, M, E) ∂ T M = B∂ X F (k, M, E) ∂ T E = ∂ X S(k, M, E) .

Analytic framework

We collect in this section some preliminaries from functional analysis that will allow us to formulate in which sense the formally derived averaged equations provide some valuable pieces of information at the spectral level.

Integral transform.

To analyze linearized lattice systems with periodic coefficients we shall make use of an adapted integral transform, usually called discrete Bloch transform. Though probably less known that its continuous counterpart it is a natural functional tool that have already received various applications in related contexts, see for instance [START_REF] Chirilus-Bruckner | Rigorous description of macroscopic wave packets in infinite periodic chains of coupled oscillators by modulation equations[END_REF][START_REF] Chirilus | Detection of standing pulses in periodic media by pulse interaction[END_REF].

3.1.1. Discrete Bloch Transform. Let N ∈ N * be given. We define the following transform

(3.1) B : 2 (Z; C d ) → L 2 [[0, N -1]] × -π N , π N ; C d , f → f where for j ∈ [[0, N -1]] (3.2) f (j, •) := lim k 0 →∞ k∈[[-k 0 ,k 0 ]]
e -i(kN +j) • f kN +j , the limit being taken in the L 2 -π N , π N -sense. Up to a multiplicative constant this is a total isometry whose reciprocal is given as

(3.3) B -1 : L 2 [[0, N -1]] × -π N , π N ; C d → 2 (Z; C d ) , f → f where for (k, j) ∈ Z × [[0, N -1]] (3.4) f kN +j = N 2π π N - π N
e i (kN +j) ξ f (j, ξ) dξ .

From now on we'll constantly identify functions on Z/(N Z) with functions on [[0, N -1]] through usual restriction/periodic extension processes. To this purpose we stress that formula (3.2) already provides the needed periodic extension. With this identification, inverse formula (3.4) simply reads for any j ∈ Z (3.5)

f j = N 2π π N - π N
e i j ξ f (j, ξ) dξ .

Discrete Bloch symbols. To any element

ξ → A ξ of L ∞ ( -π N , π N ; L( 2 (Z/(N Z); C d ))) one may associate a bounded operator A ∈ L( 2 (Z; C d )) such that, for any f ∈ 2 (Z; C d ), (Af )ˇ(•, ξ) = A ξ f (•, ξ) . Explicitly (3.6) (Af ) j = N 2π π N - π N e i j ξ (A ξ f (•, ξ)) j dξ .
Moreover up to a multiplicative constant the map

L ∞ ( -π N , π N ; L( 2 (Z/(N Z); C d ))) → L( 2 (Z; C d ))
is a partial isometry. From this follows the classical identity between the spectrum of the synthetized operator A and the essential range of ξ → σ(A ξ ). In particular, in the case where

ξ → A ξ belongs to C 0 ( -π N , π N ; L( 2 (Z/(N Z); C d ))), this yields σ(A) = ξ∈[- π N , π N ] σ(A ξ )
so that the spectrum of A may be analyzed by looking at the family of spectra of the finitedimensional symbols A ξ .

Reciprocally for a large class of A ∈ L( 2 (Z; C d )) one may define corresponding symbols ξ → A ξ formally given by

A ξ = e -i • ξ A e i • ξ that is, for f ∈ 2 (Z/(N Z); C d ) (A ξ f ) j = e -i j ξ (A (e i • ξ f )) j
with an extended sense of A. To make things more concrete, we focus now on the case where A is an N -periodic discrete differential operator in the sense that it is given as

A = a(•, T) with a(•, X) = j∈Z a j (•) X j
where T is again the left shift operator and a : Z → L(C d ) is N -periodic. The Bloch symbols of A are then given by

A ξ = a(•, e i ξ T)
where T is now the left shift operator acting on (C d ) Z/(N Z) and a is considered as a function on Z/(N Z).

3.2. Linear evolution. By linearizing original equations we receive equations in the form d dt

V (t) = A(t) V (t)
where A(•) is a 1/ω-periodic family of bounded operators generating S(•, •) an evolution system on 2 (Z; C d ) (in the sense defined for instance in [Paz83, Chapter 5]). Explicitly we have (∀s ∈ R, S(s, s) = Id) and (∀(s, r, t) ∈ R 3 , S(t, s) = S(t, r)S(r, s)) and ∀(s, t) ∈ R 2 , (∂ t S)(t, s) -A(t) S(t, s) = 0 and (∂ s S)(t, s) + S(t, s) A(s) = 0 .

Mark that periodicity implies that linear stability in the sense that

sup t∈R + S(t, 0) < ∞ is equivalent to sup k∈N S(1/ω, 0) k < ∞ .
Likewise linear exponential instability in the sense that

∃γ > 0 , inf t∈R + e -γ t S(t, 0) > 0 is equivalent to ∃κ > 1 , inf k∈N κ -k S(1/ω, 0) k > 0 ,
that is to ρ(S(1/ω, 0)) > 1.

Remark 3.1. Here we use time-translation invariance and periodicity in time as a substitute to space-translation invariance and stationarity in a moving frame that play a crucial role in similar analyses of continuous systems. However as pointed out for instance in [START_REF] Friesecke | Solitary waves on FPU lattices. I. Qualitative properties, renormalization and continuum limit[END_REF][START_REF] Friesecke | Solitary waves on FPU lattices. II. Linear implies nonlinear stability[END_REF][START_REF] Friesecke | Solitary waves on Fermi-Pasta-Ulam lattices. III. Howland-type Floquet theory[END_REF][START_REF] Friesecke | Solitary waves on Fermi-Pasta-Ulam lattices. IV. Proof of stability at low energy[END_REF] the loss of an exact invariance in space does not preclude arguments based on a pseudo-invariance. Notably, mark that for all t ∈ R +

A(t + 1/c) = T -1 A(t)T .
This implies that for all (t, s) ∈ R

S(t + 1/c, s + 1/c) = T -1 S(t, s)T .
As a consequence, linear stability is also equivalent to

sup k∈N (TS(1/c, 0)) k < ∞ .
The fact that such an alternative exist is obviously crucial when dealing with fronts, shocks or solitary waves rather than with periodic waves.

From now on, in order to make the most of the discrete Bloch transform we focus on cases where, for some N ∈ N * , A(•) is valued in N -periodic discrete differential operators that is we restrict to cases where linearization is about a periodic traveling wave with an integer period, N . Then A ξ (•) the value at ξ of the corresponding family of Bloch symbols generates the evolution system given by the value at ξ, S ξ (•, •), of the family of Bloch symbols associated with S(•, •). In our context some continuity in the ξ-variable is available and linear stability reads

sup ξ∈[- π N , π N ] sup k∈N S ξ (1/ω, 0) k < ∞ ,
while linear exponential instability may be written as ∃ξ ∈ -π N , π N , ρ(S ξ (1/ω, 0)) > 1 . This implies that a necessary condition for linear stability is that for any ξ ∈ -π N , π N σ (S ξ (1/ω, 0)) ⊂ B(0, 1) .

The fact that each periodic wave is embedded in a family of neighboring waves implies that 1 always belongs to the spectrum of S 0 (1/ω, 0). It is then crucial in the analysis of linear stability to determine how the eigenvalue 1 evolves when ξ is varied. Our goal is to prove that this piece of information is precisely given by the averaged equations.

Algebraic interplay.

To analyze the evolution of the spectrum of S ξ (1/ω, 0) when ξ is small, we first stress the following classical identity, for any ξ

∈ -π N , π N (3.7) ∀(s, t) ∈ R 2 , S ξ (t, s) = S 0 (t, s) + t s S 0 (t, r) (A ξ (r) -A 0 (r)) S ξ (r, s) dr .
In particular, expanding A ξ (•) as

A ξ (•) ξ→0 = A 0 (•) + i ξ A (1) (•) + (i ξ) 2 A (2) (•) + O(|ξ| 3 ) ,
uniformly in time, for some A (1) (•), A (2) (•), we receive locally uniformly in time

S ξ (t, s) ξ→0 = S 0 (t, s) + i ξ S (1) (t, s) + (i ξ) 2 S (2) (t, s) + O(|ξ| 3 ) , with (3.8) S (1) (t, s) = t s S 0 (t, r) A (1) (r) S 0 (r, s) dr S (2) (t, s) = t s S 0 (t, r) A (2) (r) S 0 (r, s) dr + t s
S 0 (t, r) A (1) (r) S (1) (r, s) dr .

We have already encountered a similar expansion for L ξ := e -iξ• L e iξ• . Operators L ξ acting on functions of period 1 are actually continuous Bloch symbols of the operator obtained by extending L to functions defined on the full line. Since we won't need this relation here, we mostly use the fact that by expanding

L ξ = L + (iξ k) L (1) + (iξ k) 2 L (2) + O(|ξ| 3 ) .
we obtain operators L (j) that are practical shorthands to denote relations obtained by differentiating profile equation with respect to parameters. For instance, when a denotes a parameter which is not k (but possibly ζ, to account for phase shifts), one receives

(3.9) L∂ a u = ∂ a ω∂ ζ u while (3.10) L∂ k u = -k∂ k c ∂ ζ u -L (1) ∂ ζ u .
Proofs of our main results strongly rely on the interplay between L ξ and A ξ that allows to use the previous relations. Explicitly to any function v : R → C d of period 1, we associate

V v : R → (C d ) Z by V v (j, t) = v(k (j -ct)) .
Observe that when v and w are one-periodic

1/ω 0 V w (s), V v (s) ds = 1 kω w, v .
Note also that when v is one-periodic V v (•, t) is 1/ω-periodic in time and valued in N -periodic sequences. In this case, (3.11)

d dt V v (t) -A ξ (t) V v (t) = -V L ξ v (t) -iξωV v (t) d dt V v (t) + A ξ (t) * V v (t) = V L * ξ v (t) -iξωV v (t) .
In particular

(3.12) d dt ( V w , V v )(t) = V L * ξ w (t), V v (t) -V w (t), V L ξ v (t) ,
which turns to be crucial to transport duality relations. Obviously by expanding in ξ the foregoing relations we obtain corresponding relations between operators A (j) and operators L (j) , for instance

(3.13) d dt V v (t) -A 0 (t) V v (t) = -V L v (t) -A (1) (t) V v (t) = -k V (L (1) -c) v (t) -A (2) (t) V v (t) = -k 2 V L (2) v (t) d dt V v (t) + A 0 (t) * V v (t) = V L * v (t)
. Mark that from (3.9) and (3.11)

∀(t, s) ∈ R 2 , V ∂ ζ u (t) = S 0 (t, s) V ∂ ζ u (s) .
In particular, V ∂ ζ u (0) ∈ ker(S 0 (1/ω, 0) -Id). More generally (3.9) and (3.11) provide as many members of the generalized kernel of S 0 (1/ω, 0) -Id as the expected dimension of the family of periodic traveling waves (counted up to translation). Regardless of the version of the equations that we are considering, assumption A may be reformulated as the fact that these elements form a basis of this generalized kernel. Moreover (3.10), (3.11) and (3.8) yield

(3.14) ∀(t, s) ∈ R 2 , V ∂ k u (t) = S 0 (t, s) V ∂ k u (s) -(t -s)∂ k ω V ∂ ζ u (s) + 1 k S (1) (t, s)V ∂ ζ u (s) .

Spectral validation of averaged equations

4.1. Reaction-diffusion case. Firstly we analyze the linearization of (2.3) about U (•) the periodic traveling wave associated to u k , with k = 1/N . The linearized evolution obeys (∀t , d dt V (t) = A(t) V (t)) where

A(•) = µ (T -2Id + T -1 )) + df (U (•))
with usual identification of functions with their local action on sequences, that is, here,

∀(t, j) , [df (U (t))V ] j = df (U j (t)) V j .
Explicitly here

A (1) (•) ≡ µ (T -T -1 ) and A (2) (•) ≡ 1 2 µ (T + T -1 ) .
Now our main goal is to prove the following result that validates (2.10) and (2.14) at the spectral level.

Theorem 4.1. Assuming that condition A holds in a neighborhood of u k , there exist positive 0 and ξ 0 ∈ (0, π/N ) and an analytic curve λ :

[-ξ 0 , ξ 0 ] → C such that ∀ξ ∈ [-ξ 0 , ξ 0 ] , σ(S ξ (1/ω k , 0)) ∩ B(1, 0 ) = {λ(ξ)} . Moreover, for ξ ∈ [-ξ 0 , ξ 0 ], eigenvalue λ(ξ) of S ξ (1/ω k , 0) is simple and (4.1) λ(ξ) ξ→0 = e 1 ω k (ikξ∂kω k + (ikξ) 2 d(k)) + O(|ξ| 3 ) ,
where d(k) is given by (2.13) under normalization N.

The above theorem has many direct implications. For instance, we readily see that d(k) > 0 yields a linear exponential instability caused by side-band perturbations, that is, by perturbations with arbitrary small -but non zero -spatial Floquet exponents. It also proves as expected that the single linear group velocity is indeed ∂ k ω k .

Proof. Assumption A implies that 1 is a simple eigenvalue of S 0 (1/ω, 0) so that the only point that does not follow from classical results on regular perturbations of simple eigenvalues is expansion (4.1). To identify coefficients we expand

λ(ξ) ξ→0 = 1 + i ξ λ (1) + (i ξ) 2 λ (2) + O(|ξ| 3 ) ,
and a corresponding eigenvector ϕ ξ , which also depends analytically on ξ, as

ϕ ξ ξ→0 = V ∂ ζ u (0) + iξ ϕ (1) + (iξ) 2 ϕ (2) + O(|ξ| 3 ) .
Setting p = 1/ω k for writing convenience, we now expand relation S ξ (p, 0)ϕ ξ = λ(ξ)ϕ ξ .

At first order in iξ we receive

(4.2) (S 0 (p, 0) -Id) ϕ (1) = λ (1) V ∂ ζ u (0) -S (1) (p, 0)V ∂ ζ u (0) . Hence (S 0 (p, 0) -Id) (ϕ (1) -kV ∂ k u (0)) = λ (1) -kp∂ k ω k V ∂ ζ u (0) .
Since 1 is a simple eigenvalue of S 0 (p, 0), V ∂ ζ u (0) does not belong to the image of S 0 (p, 0) -Id. From this we infer the desired

λ (1) = kp∂ k ω k and ϕ (1) -kV ∂ k u (0) ∈ ker(S 0 (p, 0) -Id) = CV ∂ ζ u (0
V ∂ ζ u , V ∂ M 1 u , • • • , V ∂ M d 1 u
) and (kV u ad , kV e 1 , • • • , kV e d 1 ) form dual bases of the generalized kernels of S 0 (1/ω k,M , 0) -Id and S 0 (1/ω k,M , 0) * -Id. By Kato's perturbation method [Kat76, pp. 99-100] we may continue these bases as dual bases (q 0 (ξ), • • • , q d 1 (ξ)) and (q 0 (ξ), • • • , qd 1 (ξ)) of the generalized eigenspaces of S ξ (1/ω k,M , 0) and S ξ (1/ω k,M , 0) * associated with their spectra in B(1, 0 ) as long as |ξ| ≤ ξ 0 provided ξ 0 and 0 are small enough. Then, for |ξ| ≤ ξ 0 , the spectrum of S ξ (1/ω k,M , 0) -Id in B(0, 0 ) is the spectrum of the matrix

Ω ξ = qα (ξ), S ξ (1/ω k,M , 0) -Id q β (ξ) 0≤α,β≤d 1 .
Observe that

Ω 0 =    0 p∂ M 1 ω • • • p∂ M d 1 ω . . . 0 d 1 ×d 1 0    where p = 1/ω k,M . Expanding Ω ξ = Ω 0 + ikξ Ω (1) + (ikξ) 2 Ω (2) + O(|ξ| 3 ) ,
we also obtain since

1 k S (1) (p, 0) V ∂ ζ u (0) = -(S 0 (p, 0) -Id) V ∂ k u (0) + p ∂ k ω k V ∂ ζ u (0) that Ω (1) =      * • • • * 0 . . . * 0      .
The upshot is that the matrix Ωξ defined by

Ωξ := 1 ikξ Σ -1 ξ Ω ξ Σ ξ with Σ ξ = (ikξ) -1 O 1×d 1 0 d 1 ×1 I d 1 ×d 1
depends analytically on ξ. The heuristics behind this transformation is that the change of basis transforms a phase-like coordinate to a wavenumber-like coordinate as in the derivation of the Whitham system (2.20) while the division by ikξ changes eigenvalues in velocities. Now the point is to identify Ω0 . We already know that

Ω0 = (Ω (1) ) α,β (Ω 0 ) α,β (Ω (2) ) α,β (Ω (1) ) α,β 0≤α,β≤d 1 . and, for β = 1, • • • , d 1 , (Ω 0 ) 1,β = p∂ M β ω. First (Ω (1) ) 0,0 = 1 ik q0 (0), S 0 (p, 0) -Id ∂ ξ q 0 (0) + 1 k q0 (0), S (1) (p, 0)q 0 (0) = p∂ k ω + 1 ik q0 (0), S 0 (p, 0) -Id (∂ ξ q 0 (0) -ikV ∂ k u (0)
) . Our claim is that the last term vanishes. Indeed since (S 0 (p, 0) -Id)q 0 (0) = 0 and, for |ξ| ≤ ξ 0 , (S ξ (p, 0) -Id)q 0 (ξ) ∈ Span({q 0 (ξ),

• • • , q d 1 (ξ)}) we conclude 1 ik (S 0 (p, 0) -Id)∂ ξ q 0 (0) + 1 k S (1) (p, 0)q 0 (0) ∈ Span({q 0 (0), • • • , q d 1 (0)}) hence S 0 (p, 0) -Id (∂ ξ q 0 (0) -ikV ∂ k u (0)) ∈ Span({q 0 (0), • • • , q d 1 (0)}) . This implies ∂ ξ q 0 (0) -ikV ∂ k u (0) ∈ Span({q 0 (0), • • • , q d 1 (0)})
which in turn yields

S 0 (p, 0) -Id (∂ ξ q 0 (0) -ikV ∂ k u (0)) ∈ Span({q 1 (0), • • • , q d 1 (0)}) ,
proving the cancellation. To make the remaining computations easier we pick (z

1 , • • • , z d 1 ) such that ∂ ξ q 0 (0) -ikV ∂ k u (0) = 1≤α≤d 1 z α q α (0)
and replace (q 0 (ξ), q1 (ξ),

• • • , qd 1 (ξ)) with (q 0 (ξ) -ξ d 1 α=1 z α q α (ξ), q1 (ξ) + ξz 1 q0 (ξ), • • • , qd 1 (ξ) + ξz d 1 q0 (ξ)). This ensures ∂ ξ q 0 (0) = ikV ∂ k u (0).
With this in hands we may also compute, for 1

≤ α, β ≤ d 1 , (Ω (1) ) α,β = 1 ik ∂ ξ qα (0), S 0 (p, 0) -Id q β (0) + 1 k qα (0), S (1) (p, 0)q β (0) = p∂ M β ω ik ∂ ξ qα (0), q 0 (0) + 1 k qα (0), S (1) (p, 0)q β (0) = - p∂ M β ω ik qα (0), ∂ ξ q 0 (0) + 1 k qα (0), S (1) (p, 0)q β (0) = -p∂ M β ω qα (0), V ∂ k u (0) + 1 k qα (0), S (1)
(p, 0)q β (0) . Now using that, for any (s, t),

S 0 (t, s)V ∂ M β u (s) = V ∂ M β u (t) + (t -s)∂ M β ωV ∂ ζ u (t) ,
we go further to obtain

1 k qα (0), S (1) (p, 0)q β (0) = p 0 V eα (t), A (1) (t)[V ∂ M β u (t) + t∂ M β ωV ∂ ζ u (t)] dt with p 0 V eα (t), A (1) (t)V ∂ M β u (t) dt = p e α , (L (1) -c)∂ M β u = p 1 0 ηe α • df (u(ζ))∂ M β u(ζ)dζ = -p ∂ M β F α (k, M )
and, by using (L

(1) -c)∂ ζ u = ∂ k ω∂ ζ u -L∂ k u and integrating by part, p 0 t V eα (t), A (1) (t)V ∂ ζ u (t)] dt = -k p 0 t V eα (t), V L∂ k u (t)] dt = k p 0 t V eα (t), d dt -A 0 (t)) V ∂ k u (t)] dt = k p V eα (0), V ∂ k u (0) . Hence, for 1 ≤ α, β ≤ d 1 , (Ω (1) ) α,β = -p ∂ M β F α (k, M ) .
Lastly we compute, for 1

≤ α ≤ d 1 , (Ω (2) ) α,0 = 1 (ik) 2 ∂ ξ qα (0), S 0 (p, 0) -Id ∂ ξ q 0 (0) + 1 ik 2 qα (0), S (1) (p, 0)∂ ξ q 0 (0) + 1 ik 2 ∂ ξ qα (0), S (1)
(p, 0) q 0 (0) + 1 k 2 qα (0), S (2) (p, 0) q 0 (0) .

On one hand, integrating by parts,

1 ik 2 ∂ ξ qα (0), S (1) (p, 0) q 0 (0) = 1 ik 2 p 0 ∂ ξ qα (0), S 0 (p, t)A (1) (t)V ∂ ζ u (t) dt = 1 ik p 0 ∂ ξ qα (0), S 0 (p, t)∂ k ω V ∂ ζ u (t) dt + 1 ik p 0 ∂ ξ qα (0), S 0 (p, t) d dt -A 0 (t) V ∂ k u (t) dt = 1 ik p∂ k ω ∂ ξ qα (0), V ∂ ζ u (0) -1 ik ∂ ξ qα (0), S 0 (p, 0) -Id ∂ ξ q 0 (0) .
On the other hand, using (3.14),

1 k 2 qα (0), [ 1 i S (1) (p, 0)∂ ξ q 0 (0) + S (2) (p, 0) q 0 (0)] = p 0 V eα (t), A (1) (t) V ∂ k u (t) + t∂ k ω V ∂ ζ u (t) dt + 1 k p 0 V eα (t), A (2) (t) V ∂ ζ u (t) dt = p e α , (L (1) -c)∂ k u + L (2) ∂ ζ u + p ∂ k ω 1 ik qα (0), ∂ ξ q 0 (0) . Hence (Ω (2) ) α,0 = p e α , (L (1) -c)∂ k u + L (2) ∂ ζ u = -p ∂ k F α (k, M ) so that as expected Ω0 = p   ∂ k ω ∂ M 1 ω • • • ∂ M d 1 ω -∂ k F α (k, M ) -∂ M 1 F α (k, M ) • • • -∂ M d 1 F α (k, M )   .
4.3. Hamiltonian case. Now we turn our attention to the role of (2.26) in the spectral analysis of the linearization of (2.21). Our analysis hereafter follows closely the proof of the foregoing subsection. The only difference stems from the presence of one extra "hidden" conservation law.

Theorem 4.3. Assuming that condition A holds in a neighborhood of u k,M,E , there exist positive 0 and ξ 0 ∈ (0, π/N ) and continuous curves λ

α : [-ξ 0 , ξ 0 ] → C, α = 1, • • • , d + 2, such that ∀ξ ∈ [-ξ 0 , ξ 0 ] , σ(S ξ (1/ω k,M , 0)) ∩ B(1, 0 ) = {λ 1 (ξ), • • • , λ d+2 (ξ)} . Moreover, for ξ ∈ [-ξ 0 , ξ 0 ] and α = 1, • • • , d + 2 (4.5) λ α (ξ) ξ→0 = e 1 ω k,M ikξ aα + O(|ξ| 2 ) ,
where a 1 , • • • , a d+2 are the characteristic speeds of the Whitham system (2.26) linearized about parameters (k, M, E).

Proof. As in the previous subsection, perturbatively we may obtain dual bases (q 0 (ξ), • • • , q d+1 (ξ)) and (q 0 (ξ), • • • , qd+1 (ξ)) of the generalized eigenspaces of S ξ (1/ω k,M,E , 0) and S ξ (1/ω k,M,E , 0) * associated with their spectra in B(1, 0 ) as long as |ξ| ≤ ξ 0 provided ξ 0 and 0 are small enough.

Here however we start from (q 0 (0), .

• • • , q d+1 (0)) = (V ∂ ζ u , V ∂ M 1 u , • • • , V ∂ M d u , V ∂ E u )
The cancellation discussed in the previous subsection shows that the matrix Ωξ defined by

Ωξ := 1 ikξ Σ -1 ξ Ω ξ Σ ξ with Σ ξ = (ikξ) -1 O 1×d+1 0 d+1×1 I d+1×d+1
depends analytically on ξ. The point is again to connect Ω0 with our averaged system, here (2.26). Manipulations as above already yield (Ω (2) ) d+1,0 = 1 (ik) 2 ∂ ξ qd+1 (0), S 0 (p, 0) -Id ∂ ξ q 0 (0) + 1 ik 2 qd+1 (0), S (1) (p, 0)∂ ξ q 0 (0) + 1 ik 2 ∂ ξ qd+1 (0), S (1) (p, 0) q 0 (0) + 1 k 2 qd+1 (0), S (2) (p, 0) q 0 (0) , and, for β = 1, • • • , d + 1, (Ω (1) ) d+1,β = 1 ik ∂ ξ qd+1 (0), S 0 (p, 0) -Id q β (0) + 1 k qd+1 (0), S (1) (p, 0)q β (0) . Moreover we may enforce ∂ ξ q 0 (0) = ik∂ k u.

Ω0 = p       ∂ k ω ∂ M 1 ω • • • ∂ M d ω ∂ E ω -B∂ k F -B∂ M 1 F • • • -B∂ M d F -B∂ E F 1 p (Ω (2) ) d+1,0 1 p (Ω (1) ) d+1,1 • • • 1 p (Ω (
For β = 1, • • • , d + 1, setting

a β = M β if 1 ≤ β ≤ d E if β = d + 1
and proceeding as in the previous subsection we derive From this stems as expected (Ω (1) ) d+1,β = p ik ∂ a β ω [ ∂ ξ qd+1 (0), q 0 (0) + qd+1 (0), ∂ ξ q 0 (0) ] + p ∂ a β S = p ∂ a β S . Lastly, on one hand, we still have 1 ik 2 ∂ ξ qd+1 (0) , S (1) (p, 0) q 0 (0) = 1 ik p ∂ k ω ∂ ξ qd+1 (0), q 0 (0) -1 ik ∂ ξ qd+1 (0), S 0 (p, 0) -Id ∂ ξ q 0 (0) . On the other hand, Hence as expected (Ω (2) ) d+1,0 = p ∂ k S. This achieves the proof of

Ω0 = p       ∂ k ω ∂ M 1 ω • • • ∂ M d ω ∂ E ω -B∂ k F -B∂ M 1 F • • • -B∂ M d F -B∂ E F ∂ k S ∂ M 1 S • • • ∂ M d S ∂ E S      
that yields the theorem.

Conclusions and Remarks

We have proved that the slow modulation ansatz captures accurately the essential spectral features of the space-time low 2 Floquet-multiplier evolution linearized about a periodic wave whose period belongs to the lattice. As already mentioned this opens at least three classes of questions.

In many cases one also expects that in the large-time the dynamics is effectively reduced to a slow modulation evolution. Our analysis provides the spectral background -including, as appears from our proofs, expansion of critical eigenfunctions -for a proof of such claims, either in the linear regime as, for continuous systems, in [Rod] for the Korteweg-de Vries equation, or in the complete nonlinear regime, as in [START_REF] Johnson | Behavior of periodic solutions of viscous conservation laws under localized and nonlocalized perturbations[END_REF] for general parabolic systems of partial differential equations.

With this in mind, one expects to be able to identify on averaged systems the key features of the discrete dynamics. It is then natural to ask if one may design discrete systems leading to desired modulation properties, either with applications to the analysis of numerical schemes where the goal is to preserve the slow modulation asymptotics at the discrete level, or in the effective conception of smart materials, a fast-growing field of investigation.

At last, even for the direct problem at the spectral level, remains the technical challenge of dealing with waves of general period, a seemingly quasi-periodic problem.

and (q 0

 0 (0), • • • , qd+1 (0)) = (kV u ad , kV e 1 , • • • , kV e d , kV δ k H[u] ) ,where u ad belongs to the generalized kernel of L * , is orthogonal to∂ M 1 u, • • • , ∂ M d u, ∂ E u and satisfies u ad , ∂ ζ u = 1. Then, for |ξ| ≤ ξ 0 , the spectrum of S ξ (1/ω k,M,E , 0) -Id in B(0, 0 ) is the spectrum of the matrix Ω ξ = qα (ξ), S ξ (1/ω k,M,E , 0) -Id q β (ξ) 0≤α,β≤d+1

,

  where p = 1/ω k,M,E , for α = 1, • • • , d + 1,

(Ω ( 1 ) 0 V

 10 ) d+1,β = p ik ∂ a β ω ∂ ξ qd+1 (0), q 0 (0) + p δ k H[u], (L (1) -c)∂ a β u + k ∂ a β ω p 0 t V δ k H[u] (t), d dt -A 0 (t) V ∂ k u (t)] dt .Now using (2.24) and (2.25) we also obtainp δ k H[u], (L (1) -c)∂ a β u = p ∂ a β S + k ∂ a β ω p δ k H[u] (t), V ∂ k u (t)] dt .

0 V 0 V

 00 1) (p, 0)∂ ξ q 0 (0) + S (2) (p, 0) q 0 (0)] = p eα (t), A (1) (t) V ∂ k u (t) + t∂ k ω V ∂ ζ u (t) dt + 1 k p eα (t), A (2) (t) V ∂ ζ u (t) dt = p e α , (L (1) -c)∂ k u + L (2) ∂ ζ u + k ∂ k ω p 0 t V δ k H[u] (t), d dt -A 0 (t) V ∂ k u (t)] dt = p ∂ k S + p ∂ k ω 1ik qα (0), ∂ ξ q 0 (0) .

  ) .Proof. Let us denote by e 1 , • • • , e d 1 both the d 1 first members of the canonical basis of C d and the constant functions with corresponding values. As a consequence of assumption A, there exists u ad belonging to the generalized kernel of L * , orthogonal to ∂ M 1 u, • • • , ∂ M d 1 u and such that u ad , ∂ ζ u = 1. Moreover, by appealing to Subsection 3.3 as in the reaction-diffusion case, one obtains that (

In a sense analogous to [Sch96, Sch98, MSU01, JNRZ14, Rod13] for dissipative continuous systems.

In the sense of close to 1.

Acknowledgement:

The main part of this work was carried out during the 3-months stay of B.K. in Lyon on 2013. B.K. would like to thank Sylvie Benzoni-Gavage for the kind invitation that was at the origin of this stay, and acknowledge the hospitality of Institut Camille Jordan at the Université Lyon 1 and the financial support of the CMIRA ACCUEIL DOC program of the region Rhône-Alpes.

By expanding now at the second order in iξ we obtain (4.3) (S 0 (p, 0) -Id) ϕ (2) = λ (2) V ∂ ζ u (0) + kp∂ k ω k ϕ (1) -S (1) (p, 0)ϕ (1) -S (2) (p, 0)V ∂ ζ u (0) .

To proceed we also need to involve u ad k defined in Subsection 2.1, that is u ad k spans ker L * and u ad k , ∂ ζ u k = 1 (in the L 2 per ([0, 1]) sense). From (3.11) we deduce that

thus that V u ad k (0) ∈ ker(S 0 (p, 0) * -Id) and from (3.12) that

, by using (3.13) and normalization N and integrating by part, we observe that

4.2. General mixed case. Now our main goal is to unravel the role of (2.20) in the spectral analysis of the linearization of (2.16).

Theorem 4.2. Assuming that condition A holds in a neighborhood of u k,M , there exist positive 0 and ξ 0 ∈ (0, π/N ) and continuous curves

are the characteristic speeds of the Whitham system (2.20) linearized about parameters (k, M).

In contrast with what happens in the reaction-diffusion case, that is averaged to a scalar equation, the above first-order expansion already provides an instability criterion. If the Whitham system (2.20) is not weakly hyperbolic at (k, M), that is, if some of the a α are not real, then the corresponding wave is linearly exponentially instable to side-band perturbations. In contrast, when weak hyperbolicity is met, then the above theorem proves that those a α do give the d 1 + 1 linear group velocities, a non trivial fact.


Tony Flores Pulgar

Jean-Luc Scharbarg

Katia Jaffrès-Runser

Christian Fraboul

Extending real-time networks over WiFi: related issues and first developments

Keywords: Hybrid real-time networks, interconnection, IEEE802.11, DCF, EDCA

The flexibility of wireless connectivity is appealing in the context of real-time industrial networks. This paper discusses the use of wireless communication protocols to interconnect remotely located fieldbuses. The focus of this paper is to analyze the feasibility and design issues related to this type of hybrid real-time network architecture. Investigations are presented by addressing an interconnection through the well-mastered WiFi technology. On an example architecture, we discuss the impact of the different distributed medium access protocols available (DCF, EDCA) on the real-time flows. We outline the main design issues related to these choices and illustrate them on a first case study where remotely located CAN buses are interconnected through an IEEE802.11g network in DCF mode. Using this very simple and cost-effective architecture, we show as a first result that transmitting soft real-time data over such an architecture is feasible.

I. INTRODUCTION

Industrial fieldbus technologies are widely rolled out to offer real-time communication capabilities on the factory floor. A large set of protocols offer deterministic and timely bounded transmissions using tailored medium access schemes and architectures (e.g. PROFIBUS, PROFINET, TTEthernet, etc.).

Recent developments for industrial communications consider introducing wireless transmissions into the global network architecture [START_REF] Pellegrini | On the use of wireless networks at low level of factory automation systems[END_REF] [START_REF] Willig | Wireless technology in industrial networks[END_REF]. First studies have assessed the capabilities of mainstream wireless technologies such as WiFi (IEEE802.11 [3]), Bluetooth or ZigBee (IEEE802.15.4) [START_REF] Willig | Wireless technology in industrial networks[END_REF] for real-time communications. In parallel, new real-time wireless protocols have been designed [START_REF] Ng | Performance evaluation of wireless controller area network (wcan) using token frame scheme[END_REF][5] [START_REF] Petersen | WirelessHART versus ISA100.11a: The format war hits the factory floor[END_REF]. Recently, two TDMA-oriented solutions (Wire-lessHART and ISA100.11a) have been commercialized for factory automation applications [START_REF] Petersen | WirelessHART versus ISA100.11a: The format war hits the factory floor[END_REF]. The main pitfall of wireless communications is of course the increased unreliability the medium suffers from due to interference and pathloss compared to shielded wires.

Among others, one of the interesting benefits of wireless transmissions is to provide a cost-effective network to interconnect distant heterogeneous or homogeneous legacy fieldbuses. The focus of this paper is to discuss this use case of wireless communications.

A wireless interconnection will benefit architectures where several fieldbuses, located far from each other, need a backhaul network to exchange data. Depending on the application, this data may be time-sensitive in the hard or in the soft sense. Either legacy wireless technologies such as WiFi (IEEE802.11 technology) or dedicated wireless protocols such as WirelessHART may be chosen, depending on the nature of the traffic exchanged between the remote buses.

For hard real-time data, a dedicated reliable wireless solution has to be picked, while for soft real-time data, a cheaper and probably less reliable wireless technology can be chosen. This paper focuses on this last case, where ofthe-shelf WiFi controllers are used to interconnect remote real-time buses. Several Medium Access Control (MAC) protocols are available with the IEEE802.11 technology. A first discussion presents the benefits and issues related to these MAC protocols to carry soft real-time data on an example hybrid network architecture.

This discussion is then illustrated with a first case study where remotely located CAN buses [START_REF] Bosch | CAN Specification version 2.0[END_REF] are interconnected through an IEEE802.11g network using a decentralized MAC protocol relying on CSMA/CA (Carrier Sense Multiple Access with Collision Avoidance). Using this very simple and cost-effective architecture, we show that transmitting soft real-time data over such an architecture is feasible, provided some design issues are correctly accounted for (e.g. gateway policies, wireless network load, interference, . . . ).

Section II illustrates the hybrid architecture we consider with an example network. Then, Section III discusses the problem of using WiFi to interconnect the remote real-time networks using different distributed IEEE802.11 MAC protocols. Next, Section IV presents the first results of extending CAN networks over IEEE802.11 DCF. Finally, Section V concludes this paper.

II. TARGET INTERCONNECTION ARCHITECTURE

An example of the type of hybrid architecture of interest in this paper is described in Figure 1. In this example, four real-time buses are interconnected through a wireless local area network that follows the mainstream IEEE802.11 standard [3].

To interconnect the real-time buses of Figure 1, four dedicated wireless gateways are implemented composed of a real-time controller and a wireless IEEE802.11 interface. Additional wireless transmitters are represented in this architecture, emitting pure wireless traffic with non real-time guarantees. All wireless transmitters may be connected in ad hoc mode (distributed medium access) or using a central controller localized in an Access Point (AP). Such an AP is not represented in Figure 1 but may be necessary Unique to pure and wireless RT flows, all frames have a critical delay which is the maximum allowed duration between their generation and reception times at their destination controllers.

III. INTERCONNECTION VIA WIFI

A. IEEE 802.11 MAC protocol overview IEEE 802.11-2012 [3] defines several standards to offer a wireless connectivity at transmission rates ranging from 11Mbps (e.g. legacy versions such as IEEE802.11b) up to 600Mbps (IEEE802.11n). Bandwidth increase is due to improvements at the physical layer (OFDM, larger bandwidth, MIMO transmissions, etc.). Next generation physical layers are expected to increase data rates beyond 600Mbps (cf. IEEE802.11ac).

The fundamental medium access in IEEE 802.11 is a Distributed Coordination Function (DCF) which is a distributed random access scheme based on CSMA/CA. Additional protocols are defined to meet specific requirements but all use the service provided by DCF as shown on the MAC architecture in Fig. 2. For instance, the Point Coordination Function is a centralized protocol where an Access Point (AP) provides a contention-free medium access. Quality of Service (QoS), where channel access is differentiated for different classes of service, is introduced with the Hybrid Coordination Function (HCF) either in a distributed manner (EDCA protocol) or in a centralized manner (HCCA protocol). Finally, wireless meshed networks architectures can be handled in a controlled fashion with the MCCA protocol. This paper discusses the use of the distributed medium access protocols (DCF and EDCA) since they are widely deployed and available for of-the-shelf designs. Thus, since these access are random based, only soft-real time data can be envisioned on this type of architecture. For more real-time constrained data, it is clear that further studies using centralized HCCA or even MCCA protocols should be performed.

B. Distributed medium access with IEEE802.11

This section introduces the main characteristics of the distributed MAC protocols of IEEE802.11, namely DCF and EDCA protocols. It discusses their implementation aspects in the context of real-time bus interconnection.

In DCF mode: a station performs carrier sensing to detect ongoing communications. If the channel is free for a period of time called Distributed InterFrame Space (DIF S), it transmits its frame immediately. If the channel is sensed busy, it defers its transmission until the end of the current transmission. Then, the station selects a random backoff b following an exponential backoff scheme. If the medium is idle for a DIF S period of time, the backoff is decremented every aSlotTime duration. The backoff interval time is decremented as long as the channel is idle and is frozen as the node detects a transmission. At the end of this transmission, when the channel remains idle during DIF S, decrementation resumes. As b reaches zero, the transmission is attempted immediately by the station. After a successful transmission, the receiver sends an ACK after a duration called Short Inter Frame Space (SIF S) If no ACK is received by the emitter after an Extended Interframe Space (EIF S), the transmission is tried again.

A new backoff b is then uniformly chosen in the range [0, w -1] where w is the contention window. This window depends on the number of failed attempts experienced by the current transmission. At the first attempt w is equal to the minimum contention window CW min . Each unsuccessful transmission involves the multiplication of w by 2 until a maximum value of CW max is reached.

DCF mode doesn't guarantee any prioritized access for real-time wireless frames. Thus, the real time traffic created by the gateways is directly competing with the pure wireless flows of non-real time nodes. In order to achieve a soft real-time guaranty, the interconnection policies implemented at the wireless gateways have to be designed with the aim of reducing the number of collisions of RT flows that occur on the wireless medium with DCF mode. This can be arranged by selecting, in a timely manner, the RT frames to be encapsulated in the wireless RT frames as we will investigate in Section IV.

In EDCA mode: QoS is being introduced on top of DCF so as to provide different levels of priority to wireless transmissions. A wireless node requests a transmission opportunity (TXOP) which is an interval of time when a particular quality-of-service (QoS) station (STA) has the right to initiate frame exchange sequences onto the wireless medium. A TXOP is defined by a starting time and a maximum duration. The TXOP is obtained by a station by successfully contending for the channel in EDCA. During an EDCA TXOP, a STA may initiate multiple frame exchange sequences.

EDCA defines four access (or traffic) categories ranging from background (lowest priority), best effort, video to voice traffic (highest priority). All nodes emitting flows in the same access category AC contend for TXOPs together using their own set of EDCA parameters which are given by the triplet (AIFS[AC], CW min [AC], CW max [AC]). The AC with the highest priority has the lowest parameter values among all ACs so as to have a higher probability to access channel first if contending with lower priority AC data frames. If there are frames with different ACs waiting in the output queue of at the same STA, collisions between them are resolved within the STA so that the data frames of higher priority compete for the channel. Remaining lower priority frames behave as if there was an external collision on the medium.

Interconnecting real-time buses using EDCA will be more challenging than using DCF because several options are available. First, besides defining encapsulation policies, the wireless RT flows have to be properly allocated to ACs. This can be done statically by the gateways. Another more subtile adaptation resides in the development of a dynamic adjustment of the EDCA parameters for each AC based on end-to-end delay statistics and missed-deadline rates, using monitoring. This last issue is an ongoing work which is not further addressed in this paper.

IV. A FIRST STUDY: CAN OVER DCF

Interconnection of a real-time bus with IEEE802.11 in DCF mode is investigated in the following. The real-time bus of interest is the well known Controller Area Network (CAN) [START_REF] Bosch | CAN Specification version 2.0[END_REF] standard.

A. Investigated architecture

CAN is one of the mainstream standards for embedded communications. Despite the fact that it has been originally developed for automotive communications, CAN has found its place in factory automation applications to handle sensor-actuator communications because of its ease of use and the low cost of its controllers. We recall that CAN medium access is based on CSMA/CR (with Collision Resolution): the start of frame transmissions on the bus are synchronous. When two or more stations start a transmission simultaneously, the one with the smallest frame identifier wins and the others stop their transmission. This mechanism guarantees strict priority order on identifiers. It implies limitations on the bandwidth and the maximal length of the bus (e.g. 1 Mbps for 40 meters). Bitstuffing is used to avoid the transmission of long sequences of bits with identical value. The computation of the frame length has to take into account these additional bits. In this paper, we use the upper bound given in [START_REF] Davis | Controller Area Network (CAN schedulability analysis: refuted, revisited and revised[END_REF]. The length C of a frame carrying x bytes of data is:

C = (55 + 10 × x) (1)
In this example, all wireless nodes (gateways, pure wireless emitters) function in ad hoc mode using DCF medium access protocol following the specification of the OFDM-PHY layer of 802.11g (20 MHz channel spacing). Since WiFi access points are static, we can consider that they are located at a distance where they can operate at the highest rate of 54 Mbps. We assume proper channel assignment has been performed so as to mitigate internode interference. In this ideal case study, transmissions are error-free and the transmission duration (in µs) at 54 Mbps of a MPDU of x bytes is derived according to [3]:

d(x) = 20 + 4 22 + 8(34 + x) 216 (2)
Following timing parameters are assumed: SIF S = 16µs, DIF S = 34µs, EIF S = 78µs, CW min = 16 and CW max = 1024.

B. Wireless gateway policies

As stated earlier, for DCF mode, interconnection policies at the gateway are crucial to ensure a timely distribution of wireless RT frames. Interconnection is done by encapsulating CAN frames of the wireless RT flows into IEEE802.11 frames. Encapsulation is chosen to facilitate the addressing of flows in the global network, and thus a basic static switching table is stored at the gateways to route wireless RT flows.

A gateway encapsulation strategy is characterized by the following parameters:

• The maximum number N ,m of CAN frames which can be encapsulated by gateway for destination gateway m in one IEEE802.11 frame, • The maximum time W max j a frame of a wireless RT flow f H j can wait in its source gateway.

C. Results

Table I summarizes the strategies analyzed in this study on the illustrative configuration of Figure 1. Each CAN bus carries 4 pure RT flows. Each bus 1, 2 and 3 are respectively the source of 8 wireless RT flows each which are all transmitted to bus 4. These wireless RT flows have a period and a critical delay of 10ms each. The size in bytes for each of these 8 flows are equally distributed between 2, 4, 6 and 8 data bytes.

The encapsulation strategies of Table I are considered, where uniform timers are assigned to the flows. T0 denotes a basic strategy where one CAN frame is encapsulated in exactly one wireless frame. T2, T3 and T4 are purely timed strategies where CAN frames wait for at most 2, 3 and 4ms respectively.

Table I SIMULATED STRATEGIES

A quantitative analysis of the proposed bridging strategies has been conducted. This analysis is based on simulations. Therefore, a home-made simulation tool has been developed using QNAP2 [START_REF] Simulog | Qnap2[END_REF]. Table II gives results concerning the utilization of the wireless medium. The study is conducted with different numbers of identical pure wireless flows f W (between 4 and 7), all of them with an average period of 2ms and for packets of 200 data bytes. The first line of the table shows the relative number of IEEE802.11 frames generated by wireless RT flows. The number of IEEE802.11 frames decreases when the value of the timer increases. Indeed, the average number of wireless RT frames encapsulated in an IEEE802.11 frame increases when the timer increases. Table II also shows the percentage of collisions (% Col) and the average delays (AvgD) of pure wireless flows for each strategy and each simulated scenario (i.e. number of pure wireless flows). Empty values mean that the wireless network is saturated and transmission delays diverge. It can be noticed that these percentages of collisions and delays are deeply linked to the number of contending IEEE802.11 frames.

Table III presents the rates of wireless RT frames which miss their deadlines. More precisely, each value in Table III is the average number of wireless RT frames that miss their deadlines when 10 6 of such frames are generated and transmitted. The main result of Table III is that it is still possible to reach a reasonably low value of missed deadlines (up to 200 every 10 6 frames) with the use of DCF when a limited number of pure wireless flows (up to 7) is present in the network.

It can be inferred from these results as well that there is a trade-off between the time spent waiting at the gateway due to W max j and the time spent contending for the wireless medium. In this very simple example, it can be seen that the best trade-off is obtained with T3

V. CONCLUSION

The motivation of this paper is to highlight that IEEE802.11 technology is a credible solution for interconnecting remote field buses when soft real-time data are exchanged between these buses. Future work needs to discuss the improvements provided by enhanced medium access protocols of IEEE802.11, including QoS capabilities and/or centralized mechanisms. Even though IEEE802.11 is the mainstream wireless technology, it is topical as well to evaluate the benefits of a dedicated wireless real-time protocol such as WirelessHART or ISA100.11a.

Figure 1 .

 1 Figure 1. Interconnection architecture overview

Figure 2 .

 2 Figure 2. IEEE802.11 MAC (source: [3])

T0 N ,m = 1 W

 1 max j = 0 ms ∀j T2 N ,m = 100 W max j = 2 ms ∀j T3 N ,m = 100 W max j = 3 ms ∀j T4 N ,m = 100 W max j = 4 ms ∀j


Jonathan Jung
email: jonathan.jung@ljll.math.upmc.fr

A two-fluid finite-volume solver based on OpenCL

In this paper, we propose a new very simple numerical method for solving liquid-gas compressible flows on two dimensional cartesian meshes. For achieving high performance, the scheme is tested on recent multicore processors and Graphics Processing Units (GPU), using the OpenCL environment. We describe how to install and to run the code for computing a shock-bubble interaction on your GPU.

1 Mathematical model and numerical strategy

Model

We are studying the numerical resolution of a two-fluid compressible fluid flow. The model is the Euler equations with an additional transport equation of a color function ϕ. The function ϕ is equal to 1 in the gas and 0 in the liquid. It allows to locate the two-fluid interface. We consider the system

∂ t W + ∂ x F (W ) + ∂ y G(W ) = 0, (1.1)
where W = (ρ, ρu, ρv, ρE, ρϕ) T , F (W ) = (ρu, ρu 2 + p, ρuv, (ρE + p)u, ρuϕ) T , G(W ) = (ρu, ρuv, ρv 2 + p, (ρE + p)v, ρvϕ) T .

The pressure law p is a stiffened gas pressure law p(ρ, e, ϕ) = (γ (ϕ) -1) ρe -γ (ϕ) p ∞ (ϕ) ,

where e = E -u 2 +v 2

2

, and

(γ, p ∞ )(ϕ) = (γ gas , p ∞,gas ), if ϕ = 1, (γ liq , p ∞,liq ), if ϕ = 0.

(1.3)

Mesh

We want solve (1.1) on the study domain [a; b] × [c; d] between time t = 0 and time t = T > 0. We consider regular subdivisions

(x i-1 2 ) 0≤i≤Nx+1 of [a; b] and (y j-1 2 ) 0≤j≤Ny+1 of [c; d] x i-1 2 = a + ih x , i = 0 • • • N x + 1, y j-1 2 = c + ih y , j = 0 • • • N y + 1,
where h x = b-a N +1 and h y = d-c M +1 are space steps in the x and y direction. The cell C i,j is the volume centrer on (x i , y j ) with

C i,j = x i-1 2 ; x i+ 1 2 × y j-
x i = x i-1 2 + x i+ 1 2 2 , i = 1 • • • N x , y j = y j-1 2 + y j+ 1 2 2 , j = 1 • • • N y .
The cell C i,j for i = 0, i = N x , j = 0 and j = N y are used to apply the boundary conditions. We consider also a sequence of times t n , n ∈ N such that the time step ∆t n := t n+1 -t n > 0.

Numerical strategy

We solve the two dimensional equations (1.1) on cartesian meshes with a dimensional splitting. Each step of the splitting is solved with a finite volume scheme. We use a Arbitrary-Lagrangian-Eulerian (ALE) method coupled with a remap step. It allows us to schwith between a Lagrangian approach at the liquid-gas interface and an Eulerian approach in the pure phases. The ALE step used a relaxation solver described in [START_REF] Helluy | Opencl simulations of two-fluid compressible flows with a random choice method[END_REF] and the remap step is done randomly. For more details, we refer to [START_REF] Bachmann | Random sampling remap for compressible two-phase flows[END_REF][START_REF] Helluy | Opencl simulations of two-fluid compressible flows with a random choice method[END_REF][START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluides[END_REF].

Initial conditions

We present a two-dimensional test that consists in simulating the interaction of a shock wave propagating in a liquid and interacting with a gas cylinder (bubble). It is a difficult problem, both numerically and physically. The initial conditions are depicted in Figure 1.1: the bubble of gas is surrounded by liquid at rest in atmospheric conditions. The computational domain is 2 meters long and 1 meter high. The gas is inside a disk (bubble) whose center is is located at (0.5, 0.5). The initial radius of the bubble is 0.4 m. A piston hits the left side at the velocity of 300 m.s -1 yielding a shock pressure of about 3 × 10 9 Pa. The initial data of this problem are given in the table of the Figure 1.1.

Boundary conditions

Top and bottom boundary conditions are set to solid walls while we use constant state boundary conditions for the left and right boundaries.

The code CLBubble

For performance reasons, we decided to implement the 2D scheme on recent multicore processor architectures, such as a Graphic Processing Unit (GPU).

What is a GPU?

A modern GPU is made of a global memory (≈ 1 GB) and compute units (≈ 28). Each compute unit is made of processing elements (≈ 10) and a local memory (≈64 kB). The same program (a kernel) can be executed on all the processing elements at the same time. There are some rules to respect. All the processing elements have access to the global memory but have only access to the local memory of their compute unit. The access to the global memory is slow while the access to the local memory is fast. The access to global memory is much faster if two neighboring processing elements read (or write) into two neighboring memory locations, in this case we speak about "coalescent memory access".

OpenCL and GPU implementation

The OpenCL [START_REF]Opencl online documentation[END_REF] implementation is described in [START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluides[END_REF][START_REF] Helluy | Opencl simulations of two-fluid compressible flows with a random choice method[END_REF]. We recall only the main lines.

• We compute the time step ∆t n . We compute a local time step (∆t n ) i,j on each cell and we use a reduction algorithm (see [START_REF] Blelloch | Scans as primitive parallel operations[END_REF]) to compute ∆t n = min i,j

(∆t n ) i,j .

• We perform the ALE-projection update in x-direction. We compute the fluxes balance in the x-direction for each cell of each row of the grid: a row or a part of a row is associated to one compute unit and one cell to one virtual processor. As of October 2012, the OpenCL implementations generally imposes a limit (typically 1024) for the number of work-items inside a work-group [START_REF]Opencl online documentation[END_REF]. This forces us to split the rows for some large computations. The values in the cells are then loaded into the local cache memory of the compute unit. It is possible to perform the ALE-projection algorithm with all the data into the cache memory in order to achieve the highest performance. The memory access are coalescent for reading and writing.

• We transpose the data matrix (exchange x and y) with an optimized memory transfer algorithm [START_REF] Ruetsch | Optimizing matrix transpose in cuda[END_REF]. The optimized algorithm includes four steps:

the data matrix are splitted into smaller tables of size 32 × 32. • We perform the ALE-projection update in y-direction. The memory access are coalescent because of the transposition,

• We transpose again the data matrix for the next time step.

The repartition of the computational time on each kernel is the following: the ALE-projection update represents 80%, the transposition 11% and the time step computation 9% of the global time computation. We observe high speedups (see [START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluides[END_REF][START_REF] Helluy | Opencl simulations of two-fluid compressible flows with a random choice method[END_REF]) with the GPU implementation. The efficiency is explained by two important points. We used an optimized transposition algorithm to have coalescent access in x and y directions and we also used a relaxation solver. With this solver, fluxes have a simpler expression than the exact Godunov's flux.

Documentation

We perform a documentation of the code with Doxygen [START_REF] Doxygen | [END_REF]. To obtain an HTML and a LATEX documentation, go in the folder Doxygen.

cd Doxygen/ and run:

doxygen D o x y f i l e
The HTML documentation is build in html/index.html and the LATEX documentation is in latex/refman.tex

Run the shock bubble interaction on GPU

In this section, we explain how to install the code. We choose CMake [START_REF] Cmake | Cmake online documentation[END_REF] to build the executable file. Then, you need CMake on your computer.

Environnement and GPU

The code was tested on different hardware. It works on LINUX and MACOS environnement. There are two important points for the code, the first one is OpenCL and the second one is the drivers for your GPU. The code can be execute on Graphics Processing Units from different brand (Nvidia, AMD, Intel).

How to compile and run the code?

The code run on LINUX or MACOS, the way to compile it is the same. There are two steps. The first one is to detect all OpenCL devices and the second one is the run the code on the chosen OpenCL device.

Choose the OpenCL device

To found all OpenCL devices, go in the folder "Found_OpenCL_Devices" cd Found_OpenCL_Devices/Make/ and you run: In this example, the GPU device is the number 1. Then, to run the code on the GPU, you choose DEVICEID = 1.

Run the code on the chosen device

To compile the code, go in the folder "Make" cd . . / . . / Make/ and run:

How to run the code with an interactive visualisation (OPENGL)?

As the computation is done on the GPU, it is also possible to visualize the data during the computation. For an interactive visualization, OpenGL, GLEW and GLUT have to be installed in your computer. Moreover, you should also have at least the version 2.8.10 of cmake to found the GLEW libraries. To compile the code, go in the folder "Make" and run:

cmake . -DOPENGL:BOOL=TRUE -DDEVICEID :

INT=1
where 1 is the number of the chosen OpenCL device.

Remark 1. It is possible that Cmake returns an error. Sometimes, it does not found the direction for the header (for example GLEW_INCLUDE_DIRS) or a library (for example GLEW_LIBRARIES). In this case, go in the file CMakeCache.txt and put the good link at the right place. You could found the good link with a "locate". When the CMakeCache.txt is completed, run again "cmake" to be sur that everything is good.

After that, you compile the code with make You can run the code on the device number "DEVICEID" with . / . . / Binary / c l b u b b l e

To stop the computation and build an output file, press "Escape".

Output

The code construct a Gmsh [START_REF] Gmsh | Gmsh online documentation[END_REF] file in the folder "Output". The file contains the value of ρ, u, v, p and ϕ on the domain at the final time.

How to modify the mesh?

You can increase or decrease the number of cell of the computation in each direction. To increase (resp. decrease) the number of cell in the x direction, increase (resp. decrease) the variable _NBLOCKSX at line 41 of file "CLFunctions.hpp". For the y direction, we can modify the variable _NBLOCKSY at line 47 of file "CLFunctions.hpp". You could also increase the number of cells in the x (resp. y) direction by increasing the number of work-item _NBWORKSX (resp. _NBWORKSY). If you modify this number, you have to be careful because there is a maximal size that is different for each GPU. The maximal size is given in the information at the beginning of the execution. It is the number referenced after

Nb Works Max=

The code should be more efficient if you put this maximal value in _NBWORKSX (resp. _NBWORKSY).

How to modify the initial condition?

You can change the initial condition by modifying the density ρ (_RHO), the x-velocity u (_U), the y-velocity v (_V), the pressure p (_P) and the gas mass fraction ϕ (_PHI) in the file "CLFunctions.hpp". The initial configuration is the following The parameters γ gas , γ liq , p ∞,gaz and p ∞,liq of the stiffened gas pressure law (1.2)-(1.3) in each phase are fixed by _GAMGAZ, _GAMLIQ, _PINFGAZ and _PINFLIQ. The variable ϕ (_PHI) allows to fix the phase (liquid or gas) of each states at the initial time. ϕ is equal to 1 in the gas and to 0 in the liquid phase.

How to change the pressure law?

It is possible to change the pressure law. You have to change three functions. The first one is the functions W2Y where we use the expression of the pressure p as a function of the density ρ, the internal energy e and the gas mass fraction ϕ. The second one is Y2W where we use the expression of the internal energy e as a function of the density ρ, the pressure p and the gas mass fraction ϕ. The third one is the expression of the sound_speed where we describe the sound speed c as a function of the density ρ, the internal energy e and gas mass fraction ϕ.

Conclusion

In this paper, we describe how you can install and run on your GPU the code CLBUBBLE to compute liquid-gas compressible flows on a two dimensional cartesian meshe.

Figure 1 . 1 :

 11 Figure 1.1: Shock bubble interaction. Description of the initial conditions on the left and initial data on the right.

  OpenCL is install and construct all the link for the OpenCL libraries. After that, run make It construct an executable file "cldevices" to list all OpenCL devices. You can show the list with . / . . / Binary / c l d e v i c e s The answer is something like this CL_PLATFORM_NAME = Apple CL_PLATFORM_VERSION = OpenCL 1 . 2 ( Jun 3 2014 1 2 : 4 3 : 4 1 ) 2 d e v i c e s found D ev ic e #0 name = I n t e l (R) Core (TM) i 5 -4250U CPU @ 1 . 3 0GHz D r i v e r v e r s i o n = 1 . 1 G l o b a l Memory (MB) : 8192 G l o b a l Memory Cache (MB) : 0 L o c a l Memory (KB) : 32 Max c l o c k (MHz) : 1300 Max Work Group S i z e : 1024 Number o f p a r a l l e l compute c o r e s : 4 D ev ic e #1 name = HD G ra ph ic s 5000 D r i v e r v e r s i o n = 1 . 2 ( Jun 9 2014 1 3 : 2 4 : 1 9 ) G l o b a l Memory (MB) : 1536 G l o b a l Memory Cache (MB) : 0 L o c a l Memory (KB) : 64 Max c l o c k (MHz) : 1200 Max Work Group S i z e : 512 Number o f p a r a l l e l compute c o r e s : 280

cmake . -DDEVICEID : INT=1 where 1

 1 is the number of the chosen OpenCL device. If the number of your GPU device is 0, replace 1 by 0. After that, you can compile the code with make You can run the code on the device number "DEVICEID" with . / . . / Binary / c l b u b b l e

Acknowledgements

The author wish to thank Philippe Helluy for many fruitful discussions.