Jean-François Coulombel
email: jean-francois.coulombel@math.univ-toulouse.fr

The Leray-Gårding method for finite difference schemes. II. Smooth crossing modes

Keywords: AMS classification: 65M06, 65M12, 35L03, 35L04 hyperbolic equations, difference approximations, stability, boundary conditions, semigroup estimates

In [Cou15] a multiplier technique, going back to Leray and Gårding for scalar hyperbolic partial differential equations, has been extended to the context of finite difference schemes for evolutionary problems. The key point of the analysis in [Cou15] was to obtain a discrete energy-dissipation balance law when the initial difference operator is multiplied by a suitable quantity. The construction of the energy and dissipation functionals was achieved in [Cou15] under the assumption that all modes were separated. We relax this assumption here and construct, for the same multiplier as in [Cou15], the energy and dissipation functionals when some modes cross. Semigroup estimates for fully discrete hyperbolic initial boundary value problems are deduced in this broader context by following the arguments of [Cou15].

Throughout this article, we keep the same notation as in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF]. We introduce the subsets of the complex plane:

U := {ζ ∈ C, |ζ| > 1} , U := {ζ ∈ C, |ζ| ≥ 1} , D := {ζ ∈ C, |ζ| < 1} , S 1 := {ζ ∈ C, |ζ| = 1} , D := D ∪ S 1 .
We let M n (K) denote the set of n × n matrices with entries in K = R or C. If M ∈ M n (C), M * denotes the conjugate transpose of M . We let I denote the identity matrix or the identity operator when it acts on an infinite dimensional space. We use the same notation x * y for the Hermitian product of two vectors x, y ∈ C n and for the Euclidean product of two vectors x, y ∈ R n . The norm of a vector x ∈ C n is |x| := (x * x) 1/2 . The induced matrix norm on M n (C) is denoted • .

The letter C denotes a constant that may vary from line to line or within the same line. The dependence of the constants on the various parameters is made precise throughout the text.

Introduction

This article is a sequel of our previous work [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] where we have developed a multiplier technique for finite difference schemes. The theory in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] encompasses the well-known example of the leapfrog scheme for the transport equation. Our main motivation was to derive stability estimates for finite difference schemes with a method that bypasses as much as possible Fourier analysis. This was a first step towards later considering multistep time integration techniques with finite volume space discretizations on unstructured meshes. We extend the results of [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] by dropping a simplicity assumption that was made in this work, which now allows us to consider crossing eigenmodes. Namely, the situation we consider here is the one where the latter crossing occurs in a smooth way. We also deal completely with the case of multistep schemes with two time levels for which the eigenmode crossing need not be smooth. In order to avoid repeating many arguments from [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF], we shall refer to this work whenever possible. We warn the reader that the introduction below is mostly the same as in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] since the considered problem is the same and we have found it easier for the reader to recall all the assumptions needed in the proof of our main result (which is Theorem 1 below). The main difference lies in the statement of Assumption 2 below.

We now set some more notation. With d ∈ N * being the considered space dimension, we let ∆x 1 , . . . , ∆x d , ∆t > 0 denote the space and time steps where the ratios, also known as the so-called Courant-Friedrichs-Lewy parameters, λ k := ∆t/∆x k , k = 1, . . . , d, are fixed positive constants. We keep ∆t ∈ (0, 1] as the only free small parameter and let the space steps ∆x 1 , . . . , ∆x d vary accordingly. The ℓ 2 -norms with respect to the space variables have been previously defined and thus depend on ∆t and the CFL parameters through the cell volume (either ∆x 2 • • • ∆x d on Z d-1 or ∆x 1 • • • ∆x d on Z d ). We always identify a sequence w indexed by either N (for time), Z d-1 or Z d (for space), with the corresponding step function. In particular, we shall feel free to take Fourier or Laplace transforms of such sequences.

For all j ∈ Z d , we set j = (j 1 , j ′ ) with j ′ := (j 2 , . . . , j d ) ∈ Z d-1 . We let p, q, r ∈ N d denote some fixed multi-integers, and define p 1 , q 1 , r 1 , p ′ , q ′ , r ′ according to the above notation. We also let s ∈ N denote some fixed integer. This article is devoted to recurrence relations of the form:

                 s+1 σ=0 Q σ u n+σ j = ∆t F n+s+1 j , j ′ ∈ Z d-1 , j 1 ≥ 1 , n ≥ 0 , u n+s+1 j + s+1 σ=0
B j 1 ,σ u n+σ 1,j ′ = g n+s+1 j , j ′ ∈ Z d-1 , j 1 = 1 -r 1 , . . . , 0 , n ≥ 0 ,

u n j = f n j , j ′ ∈ Z d-1 , j 1 ≥ 1 -r 1 , n = 0, . . . , s , (1)
where the operators Q σ and B j 1 ,σ are given by:

Q σ := p 1 ℓ 1 =-r 1 p ′ ℓ ′ =-r ′ a ℓ,σ S ℓ , B j 1 ,σ := q 1 ℓ 1 =0 q ′ ℓ ′ =-q ′ b ℓ,j 1 ,σ S ℓ . (2)
In (2), the a ℓ,σ , b ℓ,j 1 ,σ are real numbers and are independent of the small parameter ∆t (they may depend on the CFL parameters though), while S denotes the shift operator on the space grid: (S ℓ v) j := v j+ℓ for j, ℓ ∈ Z d . We have also used the short notation

p ′ ℓ ′ =-r ′ := d ν=2 pν ℓν =-rν , q ′ ℓ ′ =-q ′ := d ν=2 qν ℓν =-qν .
Namely, the operators Q σ and B j 1 ,σ only act on the spatial variable j ∈ Z d , and the index σ in (1) keeps track of the dependence of (1) on the s + 2 time levels involved at each time iteration. The numerical scheme (1) is understood as follows: one starts with ℓ 2 initial data (f 0 j ), ..., (f s j ) defined on [1 -r 1 , +∞) × Z d-1 . The source terms (F n j ) and (g n j ) in (1) are given. Assuming that the solution u has been defined up to some time index n + s, n ≥ 0, then the first and second equations in (1) should uniquely determine u n+s+1 j for j 1 ≥ 1 -r 1 , j ′ ∈ Z d-1 . The mesh cells associated with j 1 ≥ 1 correspond to the interior domain while those associated with j 1 = 1 -r 1 , . . . , 0 represent the discrete boundary. Recurrence relations of the form (1) arise when considering finite difference approximations of hyperbolic initial boundary value problems, see [START_REF] Gustafsson | Time dependent problems and difference methods[END_REF], which is our main motivation (the Dirichlet and extrapolation boundary conditions considered in [START_REF] Coulombel | The neumann numerical boundary condition for transport equations[END_REF] are typical examples). We wish to deal here simultaneously with explicit and implicit schemes and therefore make the following solvability assumption.

Assumption 1 (Solvability of (1)). The operator Q s+1 is an isomorphism on ℓ 2 (Z d ). Moreover, for all

F ∈ ℓ 2 (N * ×Z d-1 ) and for all g ∈ ℓ 2 ([1-r 1 , 0]×Z d-1 ), there exists a unique solution u ∈ ℓ 2 ([1-r 1 , +∞)× Z d-1 ) to the system Q s+1 u j = F j , j ′ ∈ Z d-1 , j 1 ≥ 1 , u j + B j 1 ,s+1 u 1,j ′ = g j , j ′ ∈ Z d-1 , j 1 = 1 -r 1 , . . . , 0 .
The first and second equations in (1) therefore uniquely determine u n+s+1 j for j 1 ≥ 1 -r 1 and j ′ ∈ Z d-1 ; one then proceeds to the following time index n + s + 2. Existence and uniqueness of a solution (u n j ) in ℓ 2 ([1-r 1 , +∞)×Z d-1 ) N to (1) follows from Assumption 1 as long as the source terms lie in the appropriate functional spaces, so the last requirement for well-posedness is continuous dependence of the solution on the three possible source terms (F n j ), (g n j ), (f n j ). This is a stability problem for which several definitions can be chosen according to the functional framework. The following one dates back to [START_REF] Gustafsson | Stability theory of difference approximations for mixed initial boundary value problems[END_REF] in one space dimension and to [START_REF] Michelson | Stability theory of difference approximations for multidimensional initialboundary value problems[END_REF] in several space dimensions.

Definition 1 (Strong stability). The finite difference approximation (1) is said to be "strongly stable" if there exists a constant C such that for all γ > 0 and all ∆t ∈ (0, 1], the solution (u n j ) to (1) with zero initial data (that is, (f

0 j ) = • • • = (f s j ) = 0 in (1)
) satisfies the estimate:

γ γ ∆t + 1 n≥s+1 ∆t e -2 γ n ∆t |||u n ||| 2 1-r 1 ,+∞ + n≥s+1 ∆t e -2 γ n ∆t p 1 j=1-r 1 u n j 1 ,• 2 ℓ 2 (Z d-1 ) ≤ C    γ ∆t + 1 γ n≥s+1 ∆t e -2 γ n ∆t |||F n ||| 2 1,+∞ + n≥s+1 ∆t e -2 γ n ∆t 0 j 1 =1-r 1 g n j 1 ,• 2 ℓ 2 (Z d-1 )    . (3)
The main contributions in [START_REF] Gustafsson | Stability theory of difference approximations for mixed initial boundary value problems[END_REF][START_REF] Michelson | Stability theory of difference approximations for multidimensional initialboundary value problems[END_REF] are to show that strong stability can be characterized by an algebraic condition which is usually referred to as the Uniform Kreiss-Lopatinskii Condition. We shall assume here from the start that (1) is strongly stable. We can thus control, for zero initial data, ℓ2 type norms of the solution to (1). Our goal, as in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF], is to understand which kind of stability estimate holds for the solution to (1) when one considers nonzero initial data (f 0 j ), . . . , (f s j ) in ℓ 2 . We are specifically interested in showing semigroup estimates for (1), that is in controlling the ℓ ∞ n (ℓ 2 j ) norm of the solution to (1) (which is stronger than the ℓ 2 n (ℓ 2 j ) control encoded in (3)). Our main assumption is the following. It is a relaxed version of the corresponding assumption in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] where the roots of the dispersion relation (4) below were assumed to be always simple.

Assumption 2 (Stability for the discrete Cauchy problem). For κ ∈ (C \ {0}) d , let us set :

Q σ (κ) := p ℓ=-r κ ℓ a ℓ,σ ,
where the coefficients a ℓ,σ are the same as in (2) and we use the classical notation

κ ℓ := κ ℓ 1 1 • • • κ ℓ d d for κ ∈ (C \ {0}) d and ℓ ∈ Z d .
Then there exists a finite number of points κ (1) , . . . , κ (K) in (S 1 ) d such that the following properties hold: (1) , . . . , κ (K) }, the roots to the dispersion relation 1 :

• if κ ∈ (S 1 ) d \ {κ
s+1 σ=0 Q σ (κ) z σ = 0 , (4)
are simple and located in D.

• if κ equals one of the κ (k) 's, the dispersion relation (4) has one multiple root z (k) ∈ D (its multiplicity is denoted m k ) and all other roots are simple.

• for all k = 1, . . . , K, there exists a neighborhood V k of κ (k) in C d and there exist holomorphic functions z 1 , . . . , z m k on V k such that

z 1 (κ (k) ) = • • • = z m k (κ (k) ) = z (k) ,
and for all κ ∈ V k , z 1 (κ), . . . , z m k (κ) are the m k roots to (4) that are close to z (k) .

Assumption 2 means that the dispersion relation (4) can have multiple roots (for stability reasons, multiple roots may only belong to D and not to S 1 ). When multiple roots occur, we only ask that the splitting of the multiple eigenvalue around each such point be smooth (analytic). The fact that we only consider one multiple root at a time is only a matter of clarity and notation. There is no doubt that more elaborate crossings (e.g., with one root remaining double along a submanifold of (S 1 ) d ) could be considered by further refining the techniques developed below. Eventually, we observe that multiple roots of the dispersion relation (4) occur for instance when one uses the Adams-Bashforth or Adams-Moulton time integration technique of order 3 or higher, see [HNW93, Chapter III] (which is the reason why extending the result of [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] was necessary). We now make the following assumption, which already appeared in [START_REF] Gustafsson | Stability theory of difference approximations for mixed initial boundary value problems[END_REF][START_REF] Michelson | Stability theory of difference approximations for multidimensional initialboundary value problems[END_REF] and several other works on the same topic.

Assumption 3 (Noncharacteristic discrete boundary). For ℓ 1 = -r 1 , . . . , p 1 , z ∈ C and η ∈ R d-1 , let us define

a ℓ 1 (z, η) := s+1 σ=0 z σ p ′ ℓ ′ =-r ′ a (ℓ 1 ,ℓ ′ ),σ e i ℓ ′ •η . ( 5
)
Then a -r 1 and a p 1 do not vanish on U × R d-1 , and they have nonzero degree with respect to z for all η ∈ R d-1 .

Our main result is comparable with [Wu95, Theorem 3.3], [CG11, Theorems 2.4 and 3.5] and [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF]. It shows that strong stability (or "GKS stability") in the sense of Definition 1 is a sufficient condition for incorporating ℓ 2 initial conditions in (1) and proving optimal semigroup estimates. Our result reads just as in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] but it now holds in the broader context of Assumption 2.

Theorem 1. Let Assumptions 1, 2 and 3 be satisfied, and assume that the scheme (1) is strongly stable in the sense of Definition 1. Then there exists a constant C such that for all γ > 0 and all ∆t ∈ (0, 1], the solution to (1) satisfies the estimate:

sup n≥0 e -2 γ n ∆t |||u n ||| 2 1-r 1 ,+∞ + γ γ ∆t + 1 n≥0 ∆t e -2 γ n ∆t |||u n ||| 2 1-r 1 ,+∞ + n≥0 ∆t e -2 γ n ∆t p 1 j 1 =1-r 1 u n j 1 ,• 2 ℓ 2 (Z d-1 ) ≤ C    s σ=0 |||f σ ||| 2 1-r 1 ,+∞ + γ ∆t + 1 γ n≥s+1 ∆t e -2 γ n ∆t |||F n ||| 2 1,+∞ + n≥s+1 ∆t e -2 γ n ∆t 0 j 1 =1-r 1 g n j 1 ,• 2 ℓ 2 (Z d-1 )    . (6)
In particular, the scheme (1) is "semigroup stable" in the sense that there exists a constant C such that for all ∆t ∈ (0, 1], the solution (u n j ) to (1) with (F n j ) = (g n j ) = 0 satisfies the estimate

sup n≥0 |||u n ||| 2 1-r 1 ,+∞ ≤ C s σ=0 |||f σ ||| 2 1-r 1 ,+∞ . (7)
The scheme (1) is also ℓ 2 -stable with respect to boundary data, see [START_REF] Trefethen | Instability of difference models for hyperbolic initial boundary value problems[END_REF]Definition 4.5], in the sense that there exists a constant C such that for all ∆t ∈ (0, 1], the solution (u n j ) to (1) with (F n j ) = (f n j ) = 0 satisfies the estimate

sup n≥0 |||u n ||| 2 1-r 1 ,+∞ ≤ C n≥s+1 ∆t 0 j 1 =1-r 1 g n j 1 ,• 2
ℓ 2 (Z d-1 ) .
Sections 2 and 3 below are devoted to the proof of Theorem 1. We follow the lines of [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] and first explain why the same multiplier as in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] yields an energy-dissipation balance law for the Cauchy problem (in the whole space) in the broader framework of Assumption 2. The analysis relies on a suitable construction of the energy and dissipation functionals, which are more involved than in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF]. The end of the proof of Theorem 1 follows [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] almost word for word. We explain where the specificity of the broader framework of Assumption 2 comes into play. In an Appendix, we deal with the specific case s = 1 (recurrence relations with two time levels) for which energy and dissipation functionals with local densities can be constructed. This gives hope to later deal with finite volume space discretization techniques on unstructured meshes.

The Leray-Gårding method for fully discrete Cauchy problems

This section is devoted to proving stability estimates for discretized Cauchy problems in the whole space Z d , which is the first step before considering the discretized initial boundary value problem (1). More precisely, we consider the simpler case of the whole space j ∈ Z d , and the recurrence relation in ℓ 2 (Z d ):

     s+1 σ=0 Q σ u n+σ j = 0 , j ∈ Z d , n ≥ 0 , u n j = f n j , j ∈ Z d , n = 0, . . . , s , (8)
where the operators Q σ are given by (2). We recall that in (2), the a ℓ,σ are real numbers and are independent of the small parameter ∆t (they may depend on the CFL parameters λ 1 , . . . , λ d ), while S denotes the shift operator on the space grid: (S ℓ v) j := v j+ℓ for j, ℓ ∈ Z d . Stability of ( 8) is defined as follows.

Definition 2 (Stability for the discrete Cauchy problem). The numerical scheme defined by (8) is (ℓ 2 -) stable if Q s+1 is an isomorphism from ℓ 2 (Z d ) onto itself, and if furthermore there exists a constant C 0 > 0 such that for all ∆t ∈ (0, 1], for all initial conditions f 0 , . . . , f s ∈ ℓ 2 (Z d ), there holds

sup n∈N |||u n ||| 2 -∞,+∞ ≤ C 0 s σ=0 |||f σ ||| 2 -∞,+∞ . (9)
Let us quickly recall, see e.g. [START_REF] Gustafsson | Time dependent problems and difference methods[END_REF], that stability in the sense of Definition 2 is in fact independent of ∆t ∈ (0, 1] (because (8) nowhere involves ∆t and the norms in (9) can be simplified on either side by the cell volume k ∆x k ), and can be characterized in terms of the uniform power boundedness of the so-called amplification matrix

A (κ) :=      -Q s (κ)/ Q s+1 (κ) . . . . . . -Q 0 (κ)/ Q s+1 (κ) 1 0 . . . 0 0 . . . . . . . . . 0 0 1 0      ∈ M s+1 (C) , (10)
where the Q σ (κ)'s are defined in (4) and where it is understood that A is defined on the largest open set of C d on which Q s+1 does not vanish. Let us also recall that if Q s+1 is an isomorphism from ℓ 2 (Z d ) onto itself, then Q s+1 does not vanish on (S 1 ) d , and therefore does not vanish on an open neighborhood of (S 1 ) d . With the above definition (10) for A , the following well-known result holds, see e.g. [START_REF] Gustafsson | Time dependent problems and difference methods[END_REF]:

Proposition 1 (Characterization of stability for the fully discrete Cauchy problem). Assume that Q s+1 is an isomorphism from ℓ 2 (Z d ) onto itself. Then the scheme (8) is stable in the sense of Definition 2 if and only if there exists a constant C 1 > 0 such that the amplification matrix A in (10) satisfies

∀ n ∈ N , ∀ κ ∈ (S 1 ) d , A (κ) n ≤ C 1 . (11)
In particular, the spectral radius of A (κ) should not be larger than 1 (the so-called von Neumann condition).

The eigenvalues of A (κ) are the roots to the dispersion relation (4). When these roots are simple for all κ ∈ (S 1 ) d , the von Neumann condition is both necessary and sufficient for stability of (8), see, e.g., [Cou13, Proposition 3]. However, Assumption 2 is more general than the situation considered in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] where the roots always remain simple. Nevertheless, since the occurence of a multiple root only occurs in the interior D and not on the boundary S 1 , we easily deduce from Assumption 2 that the matrix A (κ) in ( 10) is geometrically regular in the sense of [Cou13, Definition 3]. Hence we can still apply [Cou13, Proposition 3] and conclude that Assumption 2 implies stability for the Cauchy problem (8) (in the sense of Definition 2). As in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF], our goal now is to derive the semigroup estimate (9) not by applying Fourier transform to (8) and using uniform power boundedness of A , but rather by multiplying the first equation in (8) by a suitable local multiplier. As a warm-up, and to make things as clear as possible, we first deal with the simpler case where one only considers the time evolution and no additional space variable (the standard recurrence relations in C).

Stable recurrence relations

In this Paragraph, we consider sequences (v n ) n∈N with values in C. The index n should be thought of as the discrete time variable, which is the reason why we always write n as an exponent in order to be consistent with the notation used for discretized partial differential equations. Let then ν ≥ 1 and let a ν , . . . , a 0 be some complex numbers with a ν = 0 (in the next Paragraphs, we choose ν = s + 1). It is well known that all solutions (v n ) n∈N to the recurrence relation

∀ n ∈ N , a ν v n+ν + • • • + a 0 v n = 0 ,
are bounded if and only if the polynomial:

P(X) := a ν X ν + • • • + a 1 X + a 0 , (12)
has all its roots in D and the roots on S 1 are simple, see [START_REF] Hairer | Solving ordinary differential equations[END_REF]chapter III.3]. This is equivalent to requiring that the companion matrix (compare with (10)):

     -a ν-1 /a ν . . . . . . -a 0 /a ν 1 0 . . . 0 0 . . . . . . . . . 0 0 1 0      ∈ M ν (C) ,
be power bounded. In that case, the Kreiss matrix Theorem [START_REF] Strikwerda | A survey of the Kreiss matrix theorem for power bounded families of matrices and its extensions[END_REF] implies that the latter matrix is a contraction (it has a norm ≤ 1) for some Hermitian norm on C ν . In [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF], we have obtained some explicit construction of such a Hermitian norm and an associated dissipation functional in the case where all the roots of P in (12) are simple and located in D. The construction is based on a multiplier technique which is the discrete analogue of [Går56, Lemme 1.1]. The inconvenience of the result in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] is that even the roots in D, which are associated with an exponentially decaying behavior in time, are assumed to be simple. We suppress this technical assumption here and explain why the multiplier technique developed in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] allows to deal with the general case with multiple roots in D.

As in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF], we introduce the notation T for the shift operator in time, that is, for any sequence (v n ) n∈N , we define: (T m v) n := v n+m for all m, n ∈ N. The following Lemma is an extension of [Cou15, Lemma 1].

Lemma 1 (The energy-dissipation balance law for recurrence relations). Let P ∈ C[X] be a polynomial of degree ν, ν ≥ 1, that satisfies the following two properties:

• If P (z) = 0, then z ∈ D.
• If P (z) = 0 and z ∈ S 1 , then z is a simple root of P .

Then there exists a positive definite Hermitian form q e on C ν , and a nonnegative Hermitian form q d on C ν such that for any sequence (v n ) n∈N with values in C, there holds:

∀ n ∈ N , 2 Re T (P ′ (T) v n ) P (T) v n = ν |P (T) v n | 2 + (T -I) q e (v n , . . . , v n+ν-1 ) + q d (v n , . . . , v n+ν-1 ) .
In particular, for any sequence (v n ) n∈N that satisfies the recurrence relation

∀ n ∈ N , P (T) v n = 0 ,
the sequence (q e (v n , . . . , v n+ν-1 )) n∈N is nonincreasing.

The multiplier T P ′ (T) v n used in Lemma 1 is the same as in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF]. We shall see below in the proof why the expressions provided in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] for the energy and dissipation functions q e , q d can not cover the case of multiple roots and how they should be modified.

Proof. Let us first recall the proof in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] in the case of simple roots because this is the starting point for the general case we consider here. We therefore assume for now that P has degree ν and only has simple roots z 1 , . . . , z ν located in D. We write

P (X) = a ν j=1 (X -z j ) ,
with a = 0, and introduce the Lagrange polynomials:

∀ k = 1, . . . , ν , P k (X) := a ν j=1 j =k (X -z j ) .
Since the z j 's are pairwise distinct, the P k 's form a basis of C ν-1 [X]. Moreover, the following relation was obtained in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF]:

2 Re T (P ′ (T) v n ) P (T) v n -ν |P (T) v n | 2 = (T-I) ν k=1 |P k (T) v n | 2 + ν k=1 1 -|z k | 2 |P k (T) v n | 2 . ( 13
)
The conclusion of Lemma 1 is then obtained by introducing the energy (q e ) and dissipation (q d ) forms: ∀ (w 0 , . . . , w ν-1 ) ∈ C ν , q e (w 0 , . . . , w ν-1 ) :=

ν k=1 |P k (T) w 0 | 2 , ( 14
)
q d (w 0 , . . . , w ν-1 ) := ν k=1 (1 -|z k | 2 ) |P k (T) w 0 | 2 . ( 15
)
When the roots of P are located in D, q d is obviously nonnegative (this property does not depend on the fact that the roots are simple). When furthermore the roots of P are simple, the P k 's form a basis of C ν-1 [X] and q e is positive definite. The conclusion follows.

We now turn to the general case and therefore no longer assume that the roots of P in D are simple. For the sake of clarity, we label the pairwise distinct roots of P as z 1 , . . . , z m and let µ 1 , . . . , µ m denote the corresponding multiplicities. We thus have:

P (X) = a m j=1 (X -z j ) µ j ,
for some a = 0, and we introduce the polynomials:

∀ k = 1, . . . , m , P k (X) := a (X -z k ) µ k -1 m j=1 j =k (X -z j ) µ j .
We thus get the relation:

P ′ = m k=1 µ k P k ,
and it is a simple exercise to adapt the computation in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] to obtain the relation (compare with (13)):

2 Re T (P ′ (T) v n ) P (T) v n -ν |P (T) v n | 2 = (T -I) m k=1 µ k |P k (T) v n | 2 + m k=1 µ k (1 -|z k | 2 ) |P k (T) v n | 2 . (16)
The problem which we are facing is that there are too few polynomials P k to span the whole space C ν-1 [X]. The trick consists in adding to the energy part on the right hand side of (16) some nonnegative Hermitian forms in order to gain positive definiteness, while still keeping the corresponding dissipation form nonnegative. This "add and subtract" trick is performed below.

As long as a root z k is at least double (µ k ≥ 2), we introduce the polynomials:

∀ j = 1, . . . , µ k -1 , Q k,j (X) := a (X -z k ) j-1 m ℓ=1 ℓ =k (X -z ℓ ) µ ℓ ,
each of which being of degree ≤ ν -2. (Later we shall use the fact that X Q k,j (X) has degree ≤ ν -1.) We go back to (16) and add/subtract suitable quantities as follows:

2 Re T (P ′ (T) v n ) P (T) v n -ν |P (T) v n | 2 = (T -I)    m k=1 µ k |P k (T) v n | 2 + m k=1 µ k -1 j=1 ε µ k -j (1 -|z k | 2 ) 2(µ k -j) |Q k,j (T) v n | 2    (17) + m k=1 µ k (1 -|z k | 2 ) |P k (T) v n | 2 + m k=1 µ k -1 j=1 ε µ k -j (1 -|z k | 2 ) 2(µ k -j) |Q k,j (T) v n | 2 -|Q k,j (T) v n+1 | 2 ,
where ε > 0 is a parameter to be fixed later on (any choice 0 < ε ≤ 1/4 will do). In (17), it is understood that if µ k = 1 (that is, if the root z k is simple), then we do not add any polynomial Q k,j , the range of indices 1 ≤ j ≤ µ k -1 being empty. Moreover, we recall that if µ k ≥ 2 for some k, then we have |z k | < 1 so the coefficient of the Hermitian form |Q k,j (T) w 0 | 2 on the second line of (17) will be positive.

It remains to show that for some suitably chosen parameter ε > 0, the decomposition (17) yields the result of Lemma 1. Let us first observe that the ν polynomials

Q 1,1 , . . . , Q 1,µ 1 -1 , P 1 , . . . , Q m,1 , . . . , Q m,µm-1 , P m ,
span the space C ν-1 [X] (this is nothing but the classical Hermite interpolation problem). Since the quantity 1 -|z k | 2 is positive as long as µ k is larger than 2, any choice ε > 0 will make the Hermitian form q e defined on C ν by: ∀ (w 0 , . . . , w ν-1 ) ∈ C ν , q e (w 0 , . . . , w ν-1 ) :=

m k=1 µ k |P k (T) w 0 | 2 + m k=1 µ k -1 j=1 ε µ k -j (1 -|z k | 2 ) 2(µ k -j) |Q k,j (T) w 0 | 2 , (18)
positive definite. We thus now define a Hermitian form q d on C ν by:

∀ (w 0 , . . . , w ν-1 ) ∈ C ν , q d (w 0 , . . . , w ν-1 ) := m k=1 µ k (1 -|z k | 2 ) |P k (T) w 0 | 2 + m k=1 µ k -1 j=1 ε µ k -j (1 -|z k | 2 ) 2(µ k -j) |Q k,j (T) w 0 | 2 -|Q k,j (T) w 1 | 2 , (19)
and we are going to show that a convenient choice of ε makes q d nonnegative. (Let us observe here that it is crucial to have the degree of Q k,j less than ν -2 so that the quantity Q k,j (T) w 1 is a linear combination of w 1 , . . . , w ν-1 .) With the above definitions (18) and (19) for q e and q d , the energy balance law (17) reads as claimed in Lemma 1, so the only remaining task is to show that q d is nonnegative for a convenient choice of ε > 0.

for j = 0 obviously gives a nonnegative contribution since it only appears in the first sum). Hence we can keep only the very last term corresponding to j = µ k -1 and we have thus derived the lower bound:

µ k -1 j=1 ε µ k -j (1 -|z k | 2 ) 2(µ k -j) |Q k,j (T) w 0 | 2 -|Q k,j (T) w 1 | 2 ≥ -ε (1 -|z k | 2 ) (1 + |z k | 2 ) |Q k,µ k (T) w 0 | 2 ≥ - 1 2 (1 -|z k | 2 ) |P k (T) w 0 | 2 ,
where we have used |z k | ≤ 1 and ε ≤ 1/4 in the last inequality. Going back to the definition (19) of q d , and summing over the k's, we obtain that the Hermitian form q d is nonnegative for any choice of ε within the interval (0, 1/4]. The proof of Lemma 1 is complete.

The energy-dissipation balance for finite difference schemes

In this Paragraph, we consider the numerical scheme (8). We introduce the following notation:

L := s+1 σ=0 T σ Q σ , M := s+1 σ=0 σ T σ Q σ . (20)
Thanks to Fourier analysis, the following result will be a consequence of Lemma 1.

Proposition 2 (The energy-dissipation balance law for finite difference schemes). Let Assumptions 1 and 2 be satisfied. Then there exist a continuous coercive quadratic form E and a continuous nonnegative quadratic form D on ℓ 2 (Z d ; R) s+1 such that for all sequences (v n ) n∈N with values in ℓ 2 (Z d ; R) and for all n ∈ N, there holds

2 M v n , L v n -∞,+∞ = (s + 1) ||| L v n ||| 2 -∞,+∞ + (T -I) E(v n , . . . , v n+s ) + D(v n , . . . , v n+s ) .
In particular, for any choice of initial data f 0 , . . . , f s ∈ ℓ 2 (Z d ; R), the solution to (8) satisfies

sup n∈N E(u n , . . . , u n+s ) ≤ E(f 0 , . . . , f s ) ,
and (8) is (ℓ 2 -)stable.

Proof. We use the same notation v n for the sequence (v n j ) j∈Z d and the corresponding step function on R d whose value on the cell [j 1 ∆x 1 , (j

1 + 1) ∆x 1 ) × • • • × [j d ∆x d , (j d + 1) ∆x d ) equals v n j for any j ∈ Z d . Then Plancherel's Theorem gives the identity 2 M v n , L v n -∞,+∞ -(s + 1) ||| L v n ||| 2 -∞,+∞ = R d 2 Re T (P ′ κ (T) v n (ξ)) P κ (T) v n (ξ) -(s + 1) P κ (T) v n (ξ) 2 dξ (2 π) d , ( 21
)
where v n denotes the Fourier transform (in L 2 (R d )) of the function v n , and where we have let

P κ (z) := s+1 σ=0
Q σ κ 1 , . . . , κ d z σ , κ j := e i ξ j ∆x j ∈ S 1 , and P ′ κ (z) denotes the derivative of P κ with respect to z. The construction of the quadratic forms E and D is made, as in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF], of the superposition of appropriate energy and dissipation Hermitian forms for each frequency κ ∈ (S 1 ) d , each coordinate κ j being a placeholder for exp(i ξ j ∆x j ). Here, unlike [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF], the polynomial P κ either only has simple roots in D or it has one multiple root in D and all other roots are simple. We cannot therefore construct the energy and dissipation forms in a unified manner. Below we shall use the analysis of Lemma 1 in the neighborhood of finitely many points in (S 1 ) d where P κ has a multiple root and we shall use [Cou15, Lemma 1] in the neighborhood of all points where P κ only has simple roots. (This is the reason why we have recalled the proof of Lemma 1 in the case where all roots are simple.) We shall eventually glue things together thanks to a suitable partition of unity.

Let us first consider the point κ (1) ∈ (S 1 ) d for which P κ (1) has one multiple root (of multiplicity m 1 ) in D and in the neighborhood of which we have a smooth splitting of the eigenmodes z 1 , . . . , z m 1 . The other roots z m 1 +1 , . . . , z s+1 are simple and can thus be determined holomorphically with respect to κ in the neighborhood of κ (1) . Keeping in mind that the dominant coefficient of the polynomial P κ (z) equals Q s+1 (κ) (which is nonzero for κ ∈ (S 1 ) d ), we consider some κ ∈ (S 1 ) d sufficiently close to κ (1) and introduce the Lagrange polynomials:

∀ k = 1, . . . , s + 1 , P k,κ (z) := Q s+1 (κ) s+1 j=1 j =k z -z j (κ) .
We then introduce the following energy and dissipation Hermitian forms on C s+1 (below, κ always denotes an element of (S 1 ) d that is sufficiently close to κ (1) so that all considered quantities are well-defined): ∀ (w 0 , . . . , w s ) ∈ C s+1 , q e,κ (w 0 , . . . , w s ) :=

s+1 k=1 |P k,κ (T) w 0 | 2 (22) + m 1 k=1 m 1 -1 j=1 ε m 1 -j 1 -|z k (κ)| 2 2(m 1 -j) Q s+1 (κ) T -z k (κ) j-1 s+1 ℓ=m 1 +1 T -z ℓ (κ) w 0 2 , q d,κ (w 0 , . . . , w s ) := s+1 k=1 1 -|z k (κ)| 2 |P k,κ (T) w 0 | 2 (23) + m 1 k=1 m 1 -1 j=1 ε m 1 -j 1 -|z k (κ)| 2 2(m 1 -j) × Q s+1 (κ) T -z k (κ) j-1 s+1 ℓ=m 1 +1 T -z ℓ (κ) w 0 2 -Q s+1 (κ) T -z k (κ) j-1 s+1 ℓ=m 1 +1 T -z ℓ (κ) w 1 2 ,
where ε > 0 is a parameter to be fixed later on. Using the decomposition (13) which we have recalled in the proof of Lemma 1, we have the decomposition 2 Re T (P ′ κ (T) w 0 ) P κ (T) w 0 -(s+1) |P κ (T) w 0 | 2 = (T-I) (q e,κ (w 0 , . . . , w s )) + q d,κ (w 0 , . . . , w s ) , (24) for all vectors (w 0 , . . . , w s ) ∈ C s+1 , because we have just added and subtracted some Hermitian forms to the energy-dissipation balance law (13). It remains to prove that q d,κ in ( 23) is nonnegative and that q e,κ in ( 22) is positive definite. Let us start with q e,κ . If κ does not equal κ (1) , we know from Assumption 2 that the roots z 1 (κ), . . . , z s+1 (κ) are pairwise distinct so the Lagrange polynomials P k,κ form a basis of C s [X]. Hence q e,κ in ( 22) is positive definite because we have added a nonnegative form to a positive definite one. We thus now consider the case κ = κ (1) for which the m 1 first roots z 1 , . . . , z m 1 all collapse to z (1) and the m 1 first Lagrange polynomials P 1,κ (1) , . . . , P m 1 ,κ (1) are all equal. At the base point κ = κ (1) , the definition (22) thus reduces to:

q e,κ (1) (w 0 , . . . , w s ) = m 1 |P 1,κ (1) (T) w 0 | 2 + s+1 k=m 1 +1 |P k,κ (1) (T) w 0 | 2 + m 1 m 1 -1 j=1 ε m 1 -j 1 -|z (1) | 2 2(m 1 -j) Q s+1 (κ (1) ) T -z (1) j-1 s+1 ℓ=m 1 +1 T -z ℓ (κ (1) ) w 0 2
, which (up to the harmless positive constant m 1 in the second line) coincides with our definition of the Hermitian form in (18). Since the polynomials:

P 1,κ (1) (X) , P m 1 +1,κ (1) (X) , . . . , P s+1,κ (1) (X) , Q s+1 (κ (1) ) s+1 ℓ=m 1 +1 X -z ℓ (κ (1) ) , Q s+1 (κ (1) ) X -z (1) s+1 ℓ=m 1 +1 X -z ℓ (κ (1) ) , . . . , Q s+1 (κ (1) ) X -z (1) m 1 -2 s+1 ℓ=m 1 +1 X -z ℓ (κ (1) ) ,
form a basis of C s [X] (this is again the classical Hermite interpolation problem), the form q e,κ (1) is positive definite as long as the parameter ε is a fixed positive constant (the choice ε = 1/8 that is made below will do). Moreover, once ε is fixed, the form q e,κ depends in a C ∞ way on κ in the neighborhood of κ (1) .

We now show that the form q d,κ in ( 23) is nonnegative for a well-chosen parameter ε > 0 and κ ∈ (S 1 ) d sufficiently close to κ (1) . The argument is quite similar to what we have done in the proof of Lemma 1 but we now need to take into account that the m 1 first eigenmodes z 1 , . . . , z m 1 split for κ = κ (1) , which will make us choose ε > 0 slightly smaller than in the proof of Lemma 1 in order to absorb an additional error. Before going on, let us recall that the eigenmodes z 1 (κ), . . . , z s+1 (κ) belong to D for κ ∈ (S 1 ) d close to κ (1) with κ = κ (1) . By continuity, this implies that they also belong to D for κ = κ (1) . Hereafter, we shall consider κ ∈ (S 1 ) d close to κ (1) and shall therefore feel free to use the inequality |z ℓ (κ)| ≤ 1 for all ℓ = 1, . . . , s + 1 (the so-called von Neumann condition).

Let us consider some vector (w 0 , . . . , w s ) ∈ C s+1 and let us introduce the notation:

∀ k, j = 1, . . . , m 1 , W k,j := Q s+1 (κ) T -z k (κ) j-1 s+1 ℓ=m 1 +1 T -z ℓ (κ) w 0 , ( 25
)
where the complex numbers W k,j (which, according to (25), are linear combinations of w 0 , . . . , w s ) also depend on κ but there is no need to keep track of this in what follows. We start from the definition (23)

and derive the lower bound:

q d,κ (w 0 , . . . , w s ) ≥ m 1 k=1 1 -|z k (κ)| 2 |P k,κ (T) w 0 | 2 + m 1 k=1 m 1 -1 j=1 ε m 1 -j 1 -|z k (κ)| 2 2(m 1 -j) |W k,j | 2 -|W k,j+1 + z k (κ) W k,j | 2 .
Expanding the square modulus |W k,j+1 + z k (κ) W k,j | 2 and using Young's inequality under the form:

2 Re z k (κ) W k,j W k,j+1 ≤ 1 2 1 -|z k (κ)| 2 |W k,j | 2 + 2 |z k (κ)| 2 1 -|z k (κ)| 2 |W k,j+1 | 2 ≤ 1 2 1 -|z k (κ)| 2 |W k,j | 2 + 1 + |z k (κ)| 2 1 -|z k (κ)| 2 |W k,j+1 | 2 ,
we get:

q d,κ (w 0 , . . . , w s ) ≥ m 1 k=1 1 -|z k (κ)| 2 |P k,κ (T) w 0 | 2 + m 1 k=1 m 1 -1 j=1 ε m 1 -j 1 -|z k (κ)| 2 2(m 1 -j) 1 -|z k (κ)| 2 2 |W k,j | 2 - 2 1 -|z k (κ)| 2 |W k,j+1 | 2 .
After shifting indices, we end up with:

q d,κ (w 0 , . . . , w s ) ≥ m 1 k=1 1 -|z k (κ)| 2 |P k,κ (T) w 0 | 2 + m 1 k=1 m 1 -1 j=1 ε m 1 -j 2 1 -|z k (κ)| 2 2(m 1 -j)+1 |W k,j | 2 - m 1 k=1 m 1 j=2 2 ε ε m 1 -j 1 -|z k (κ)| 2 2(m 1 -j)+1 |W k,j | 2 .
Instead of choosing ε ∈ (0, 1/4] as in the proof of Lemma 1, we make the more restrictive choice ε ∈ (0, 1/8] and thus obtain:

q d,κ (w 0 , . . . , w s ) ≥ m 1 k=1 1 -|z k (κ)| 2 |P k,κ (T) w 0 | 2 - 1 -|z k (κ)| 2 4 |W k,m 1 | 2 + m 1 k=1 m 1 -1 j=1 ε m 1 -j 4 1 -|z k (κ)| 2 2(m 1 -j)+1 |W k,j | 2 . ( 26
)
We go back to the definition of the Lagrange polynomial P k,κ and of the complex numbers W k,j . For k = 1, . . . , m 1 , we have:

P k,κ (z) = Q s+1 (κ) m 1 j=1 j =k z -z j (κ) s+1 ℓ=m 1 +1 z -z ℓ (κ) .
The goal is to absorb in (26) the only negative term by means of all other positive quantities. To do this, we observe that we can expand the polynomial

X -z k (κ) m 1 -1 ,
on the basis of C m 1 -1 [X] formed by the polynomials :

1 , X -z k (κ) , X -z k (κ) m 1 -2 , m 1 j=1 j =k
X -z j (κ) .

The linear system for determining the coefficients is lower triangular and has determinant 1 so we can write for each k = 1, . . . , m 1 :

X -z k (κ) m 1 -1 = m 1 j=1 j =k X -z j (κ) + m 1 -1 j=1 a k,j (κ) X -z k (κ) j-1 , (27)
with holomorphic functions a k,j defined in the neighborhood of κ (1) and that vanish at κ (1) . The decomposition ( 27) gives (just use the definition (25) and the expression of the Lagrange polynomial P k,κ ):

W k,m 1 = P k,κ (T) w 0 + m 1 -1 j=1 a k,j (κ) W k,j ,
and we now apply the Cauchy-Schwarz inequality twice to get:

|W k,m 1 | 2 ≤ 2 |P k,κ (T) w 0 | 2 + 2 (m 1 -1) m 1 -1 j=1 |a k,j (κ)| 2 |W k,j | 2 .
Fixing from now on ε = 1/8 and using the latter inequality in (26), we find that q d,κ is nonnegative for κ sufficiently close to κ (1) (recall that |z k (κ)| < 1 uniformly with respect to κ in the neighborhood of κ (1) since the multiple eigenvalue z (1) lies in D). Moreover, we observe on the defining equation ( 23)that the Hermitian form q d,κ depends in a C ∞ way on κ in the neighborhood of κ (1) .

The above analysis close to κ (1) can be repeated word for word in the neighborhood of any other point κ (2) , . . . , κ (K) where the dispersion relation (4) has a multiple root. Now, if κ ∈ (S 1 ) d is such that the dispersion relation (4) only has simple roots at κ = κ, the analysis is much simpler since we know in that case that the roots z 1 , . . . , z s+1 locally depend holomorphically on κ and the energy and dissipation forms can be simply defined as: ∀ (w 0 , . . . , w s ) ∈ C s+1 , q e,κ (w 0 , . . . , w s ) :=

s+1 k=1 |P k,κ (T) w 0 | 2 , q d,κ (w 0 , . . . , w s ) := s+1 k=1 1 -|z k (κ)| 2 |P k,κ (T) w 0 | 2 ,
with the same notation as above for the Lagrange polynomials P k,κ . At this stage, we have shown that for any base point κ in the compact manifold (S 1 ) d , there exists an open neighborhood V of κ in (S 1 ) d and there exists a C ∞ mapping q e,κ , resp. q d,κ , on V with values in the set of positive definite, resp. nonnegative, Hermitian forms, such that the decomposition (24) holds for all κ ∈ V and all vectors (w 0 , . . . , w s ) ∈ C s+1 . By compactness of (S 1 ) d , we can take a finite covering of (S 1 ) d by such neighborhoods and glue the local definitions of the energy and dissipation forms thanks to a subordinate partition of unity. We have thus constructed a positive definite, resp. nonnegative, Hermitian form q e,κ , resp. q d,κ , on C s+1 which depends in a C ∞ way on κ ∈ (S 1 ) d and such that there holds:

2 M v n , L v n -∞,+∞ -(s + 1) ||| L v n ||| 2 -∞,+∞ = (T -I) R d q e,κ v n (ξ), . . . , v n+s (ξ) dξ (2 π) d + R d q d,κ v n (ξ), . . . , v n+s (ξ) dξ (2 π) d ,
where we recall that κ is a placeholder for (exp(i ξ 1 ∆x 1 ), . . . , exp(i ξ d ∆x d )). The conclusion of Proposition 2 follows as in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] by a standard compactness argument for showing continuity of the quadratic forms E and D, and coercivity for E.

The C ∞ regularity of the Hermitian forms q e,κ , q d,κ with respect to κ is not needed in the proof of Proposition 2 (continuity with respect to κ would be enough) but we have paid attention to that particular issue since it is a crucial step for later extending this construction to variable coefficients problems and applying symbolic calculus rules as in [START_REF] Lax | On stability for difference schemes: A sharp form of G ȧrding's inequality[END_REF]. This is left to a future work.

Semigroup estimates for fully discrete initial boundary value problems

It remains to prove Theorem 1 with the help of Proposition 2. The strategy is exactly the same as in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] since the analysis in that earlier work shows that the cornerstone of the proof of Theorem 1 is the existence of a multiplier for the fully discrete Cauchy problem on Z d . Let us emphasize that the relation (21) is of the exact same form as in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF]. The multiplier M v n is the same. The only difference is in the definition of the energy and dissipation forms E and D, but their precise expression is not useful in what follows. What matters is that D is nonnegative, and E is coercive and therefore yields a control of ℓ 2 norms on Z d . Hence we can apply the same arguments as in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] as long as the proof of Theorem 1 only uses the result of Proposition 2 and not the behavior of the roots of the dispersion relation (4). We thus follow the proof of [Cou15, Theorem 1] and explain where the same arguments can be applied without any modification.

The case with zero initial data

The first step in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] is to prove the validity of (6) for zero initial data (

f 0 = • • • = f s = 0 in (1)).
This part of the proof only uses the relation ( 21) and the fact that the multiplier M has the same stencil as the original difference operator L. Hence we can repeat the arguments in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] word for word and obtain the validity of (6) when the iteration (1) is considered with zero initial data. It then remains to consider (1) with nonzero initial data in ℓ 2 and zero interior/boundary forcing terms.

Construction of dissipative boundary conditions

This was the most technical step of the analysis in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF]. The goal here is to construct an auxiliary set of numerical boundary conditions for which, with arbitrary initial data in ℓ 2 , we can derive an optimal semigroup estimate and a trace estimate for the solution. Our result here is the same as in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] but it now holds in the broader framework of Assumption 2. (Theorem 2 is the place where Assumption 3 is needed.)

Theorem 2. Let Assumptions 1, 2 and 3 be satisfied. Then for all P 1 ∈ N, there exists a constant C P 1 > 0 such that, for all initial data f 0 , . . . , f s ∈ ℓ 2 (Z d ) and for all source term (g n j ) j 1 ≤0,j ′ ∈Z d-1 ,n≥s+1 that satisfies the integrability condition:

∀ Γ > 0 , n≥s+1 e -2 Γ n j 1 ≤0 g n j 1 ,• 2 ℓ 2 (Z d-1 ) < + ∞ ,
there exists a unique sequence

(u n j ) j∈Z d ,n∈N in ℓ 2 (Z d ) N solution to the iteration      L u n j = 0 , j ∈ Z d , j 1 ≥ 1 , n ≥ 0 , M u n j = g n+s+1 j , j ∈ Z d , j 1 ≤ 0 , n ≥ 0 , u n j = f n j , j ∈ Z d , n = 0, . . . , s . (28)
Moreover for all γ > 0 and ∆t ∈ (0, 1], this solution satisfies

sup n≥0 e -2 γ n ∆t ||| u n ||| 2 -∞,+∞ + γ γ ∆t + 1 n≥0 ∆t e -2 γ n ∆t ||| u n ||| 2 -∞,+∞ + n≥0 ∆t e -2 γ n ∆t P 1 j 1 =1-r 1 u n j 1 ,• 2 ℓ 2 (Z d-1 ) ≤ C P 1    s σ=0 ||| f σ ||| 2 -∞,+∞ + n≥s+1 ∆t e -2 γ n ∆t j 1 ≤0 g n j 1 ,• 2 ℓ 2 (Z d-1 )    . (29)
Proof. Unsurprisingly, most of the proof of Theorem 2 is the same as in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] but there is one specific point where the behavior of the roots to the dispersion relation ( 4) is used so we review the main steps of the proof and simply refer to [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] when no modification is needed. First, the existence and uniqueness of a solution to (28) follows from the invertibility of Q s+1 on ℓ 2 (Z d ). Then, using Proposition 1, we can derive the same estimate as in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] for the solution to (28):

sup n≥0 e -2 γ n ∆t ||| u n ||| 2 -∞,+∞ + γ γ ∆t + 1 n≥0 ∆t e -2 γ n ∆t ||| u n ||| 2 -∞,+∞ + n≥0 ∆t e -2 γ (n+s+1) ∆t j 1 ∈Z L u n j 1 ,• 2 ℓ 2 (Z d-1 ) ≤ C    s σ=0 ||| f σ ||| 2 -∞,+∞ + n≥s+1 ∆t e -2 γ n ∆t j 1 ≤0 g n j 1 ,• 2 ℓ 2 (Z d-1 )    , (30)
where the constant C is independent of γ, ∆t and on the source terms in (28). It remains to derive the trace estimate for the solution (u n j ) to (28) (that is showing that the third term in the sum on the left hand side of the inequality (29) is controlled by the right hand side).

The derivation of the trace estimate when γ ∆t is large enough is done as in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] since it only uses the invertibility of the operator Q s+1 on ℓ 2 (Z d ). We can thus assume from now on γ ∆t ∈ (0, ln R 0 ] for some fixed constant R 0 > 1. Then we can deduce from (30) that for any j 1 ∈ Z, the Laplace-Fourier transform u j 1 of the step function

u j 1 : (t, y) ∈ R + × R d-1 -→ u n j if (t, y) ∈ [n ∆t, (n + 1) ∆t) × d k=2 [j k ∆x k , (j k + 1) ∆x k ) ,
is well-defined on the half-space {τ ∈ C , Re τ > 0} × R d-1 . The dual variables to (t, y) are denoted τ = γ+i θ, γ > 0, and η = (η 2 , . . . , η d ) ∈ R d-1 . We also use below the notation η ∆ := (η 2 ∆x 2 , . . . , η d ∆x d ).

The following result, which is proved in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF], is used here as a blackbox since its proof is merely based on the validity of (30) and Plancherel's Theorem.

Lemma 2. With R 0 > 1 fixed as above, there exists a constant C > 0 such that for all γ > 0 and ∆t ∈ (0, 1] satisfying γ ∆t ∈ (0, ln R 0 ], there holds

j 1 ∈Z R×R d-1 p 1 ℓ 1 =-r 1 a ℓ 1 e (γ+i θ) ∆t , η ∆ u j 1 +ℓ 1 (γ + i θ, η) 2 dθ dη + j 1 ≤0 R×R d-1 p 1 ℓ 1 =-r 1 e (γ+i θ) ∆t ∂ z a ℓ 1 e (γ+i θ) ∆t , η ∆ u j 1 +ℓ 1 (γ + i θ, η) 2 dθ dη ≤ C    s σ=0 ||| f σ ||| 2 -∞,+∞ + n≥s+1 ∆t e -2 γ n ∆t j 1 ≤0 g n j 1 ,• 2 ℓ 2 (Z d-1 )    . (31)
Recall that the functions a ℓ 1 , ℓ 1 = -r 1 , . . . , p 1 , are defined in (5).

The conclusion now relies on the following crucial result. (This is the place where the behavior of the roots to the dispersion relation (4) matters, and where we therefore need to be careful.)

Lemma 3 (The trace estimate). Let Assumptions 1, 2 and 3 be satisfied. Let R 0 > 1 be fixed as above and let P 1 ∈ N. Then there exists a constant C P 1 > 0 such that for all z ∈ U with |z| ≤ R 0 , for all η ∈ R d-1 and for all sequence (w j 1 ) j 1 ∈Z ∈ ℓ 2 (Z; C), there holds

P 1 j 1 =-r 1 -p 1 |w j 1 | 2 ≤ C P 1    j 1 ∈Z p 1 ℓ 1 =-r 1 a ℓ 1 (z, η ∆ ) w j 1 +ℓ 1 2 + j 1 ≤0 p 1 ℓ 1 =-r 1 z ∂ z a ℓ 1 (z, η ∆ ) w j 1 +ℓ 1 2    . (32)
As in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF], Lemma 3 yields the conclusion of Theorem 2 by integrating (32) for the sequence ( u j 1 (γ + i θ, η)) j 1 ∈Z with respect to (θ, η) (taking z = e (γ+i θ) ∆t accordingly), using the inequality (31) from Lemma 2 and applying Plancherel's Theorem. We thus focus on the proof of Lemma 3 from now on.

Proof of Lemma 3. We reproduce most of the proof that can already be found in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] in order to highlight where Assumption 2 (in its new form) is used. We argue by contradiction and assume that the conclusion to Lemma 3 does not hold. Therefore, up to normalizing and extracting subsequences, there exist three sequences (indexed by k ∈ N):

• a sequence (w k ) k∈N with values in ℓ 2 (Z; C) such that (w k -r 1 -p 1 , . . . , w k P 1 ) belongs to the unit sphere of C P 1 +r 1 +p 1 +1 for all k, and (w k -r 1 -p 1 , . . . , w k P 1 ) converges towards (w -r 1 -p 1 , . . . , w P 1 ) as k tends to infinity,

• a sequence (z k ) k∈N with values in U ∩ {ζ ∈ C , |ζ| ≤ R 0 }, which converges towards z ∈ U ,

• a sequence (η k ) k∈N with values in [0, 2 π] d-1 , which converges towards η ∈ [0, 2 π] d-1 , and these sequences satisfy:

lim k→+∞ j 1 ∈Z p 1 ℓ 1 =-r 1 a ℓ 1 (z k , η k ) w k j 1 +ℓ 1 2 + j 1 ≤0 p 1 ℓ 1 =-r 1 z k ∂ z a ℓ 1 (z k , η k ) w k j 1 +ℓ 1 2 = 0 . (33)
We are going to show that (33) implies that the vector (w -r 1 -p 1 , . . . , w P 1 ) must be zero, which will yield a contradiction since this vector has norm 1.

• We already know that (w k -r 1 -p 1 , . . . , w k P 1 ) converges towards (w -r 1 -p 1 , . . . , w P 1 ) as k tends to infinity, and arguing by induction as in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF], we can show that (33) and Assumption 3 imply that each component (w k j 1 ) k∈N , j 1 ∈ Z, has a limit as k tends to infinity. This limit is denoted w j 1 for any j 1 ∈ Z. Then (33) implies that the sequence w, which does not necessarily belong to ℓ 2 (Z; C), satisfies the two recurrence relations (observe that the recurrence relation (35) only holds on (-∞, 0) and not on Z):

∀ j 1 ∈ Z , p 1 ℓ 1 =-r 1 a ℓ 1 (z, η) w j 1 +ℓ 1 = 0 , (34)
∀ j 1 ≤ 0 ,

p 1 ℓ 1 =-r 1 z ∂ z a ℓ 1 (z, η) w j 1 +ℓ 1 = 0 . (35)
• We define the source terms:

∀ j 1 ∈ Z , F k j 1 := p 1 ℓ 1 =-r 1 a ℓ 1 (z k , η k ) w k j 1 +ℓ 1 , G k j 1 := p 1 ℓ 1 =-r 1 z k ∂ z a ℓ 1 (z k , η k ) w k j 1 +ℓ 1 ,
which, according to (33), satisfy lim

k→0 j 1 ∈Z | F k j 1 | 2 = 0 , lim k→0 j 1 ≤0 | G k j 1 | 2 = 0 . (36)
We also introduce the vectors (here T denotes transposition)

∀ j 1 ∈ Z , W k j 1 := w k j 1 +p 1 , . . . , w k j 1 +1-r 1 T , W j 1 := w j 1 +p 1 , . . . , w j 1 +1-r 1 T ,
and the matrices:

L(z, η) :=      -a p 1 -1 (z, η)/a p 1 (z, η) . . . . . . -a -r 1 (z, η)/a p 1 (z, η) 1 0 . . . 0 0 . . . . . . . . . 0 0 1 0      ∈ M p 1 +r 1 (C) , (37)
M(z, η) :=      -∂ z a p 1 -1 (z, η)/∂ z a p 1 (z, η) . . . . . . -∂ z a -r 1 (z, η)/∂ z a p 1 (z, η) 1 0 . . . 0 0 . . . . . . . . . 0 0 1 0      ∈ M p 1 +r 1 (C) . (38)
The matrix L is well-defined on U × R d-1 thanks to Assumption 3. The matrix M is also well-defined on U × R d-1 because for any η ∈ R d-1 , Assumption 3 asserts that a p 1 (•, η) is a nonconstant polynomial whose roots lie in D. From the Gauss-Lucas Theorem, the roots of ∂ z a p 1 (•, η) lie in the convex hull of those of a p 1 (•, η), hence in D. Therefore ∂ z a p 1 (•, η) does not vanish on U . In the same way, ∂ z a -r 1 (•, η) does not vanish on U . With our above notation, the vectors W k j 1 , W j 1 , satisfy the one step recurrence relations:

∀ j 1 ∈ Z , W k j 1 +1 = L(z k , η k ) W k j 1 + F k j 1 +1 /a p 1 (z k , η k ), 0, . . . , 0 T , (39)
W j 1 +1 = L(z, η) W j 1 , (40)
∀ j 1 ≤ -1 , W k j 1 +1 = M(z k , η k ) W k j 1 + G k j 1 +1 /(z k ∂ z a p 1 (z k , η k )), 0, . . . , 0 T , (41)
W j 1 +1 = M(z, η) W j 1 . (42)
The recurrence relations (40), (42) are just an equivalent way of writing (34), (35).

• From Assumption 3 and the above application of the Gauss-Lucas Theorem, we already know that both matrices L(z, η) and M(z, η) are invertible for (z, η) ∈ U × R d-1 . Furthermore, a quick analysis shows that κ ∈ C \{0} is an eigenvalue of L(z, η) if and only if z is a solution to the dispersion relation (4). Assumption 2 therefore shows that L(z, η) has no eigenvalue on S 1 for (z, η) ∈ U × R d-1 for otherwise the von Neumann condition would not hold. (This eigenvalue splitting property dates back at least to [START_REF] Kreiss | Stability theory for difference approximations of mixed initial boundary value problems[END_REF].) However, central eigenvalues on S 1 may occur for L when z belongs to S 1 (see [START_REF] Coulombel | Stability of finite difference schemes for hyperbolic initial boundary value problems[END_REF] for a thorough analysis of the leap-frog scheme).

As in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF], the crucial point for proving Lemma 3 is that Assumption 2 in its new form still precludes central eigenvalues of M for all z ∈ U . Namely, let us show that for all z ∈ U and all η ∈ R d-1 , M(z, η) has no eigenvalue on S 1 . This property holds because otherwise, for some (z, η) ∈ U × R d-1 , there would exist a root κ 1 ∈ S 1 to the characteristic polynomial of M(z, η), that is (up to multiplying by a nonzero factor):

p 1 ℓ 1 =-r 1 z ∂ z a ℓ 1 (z, η) κ ℓ 1 1 = 0 .
For convenience, the coordinates of η are denoted (η 2 , . . . , η d ). Using the definition (5) of a ℓ 1 , and defining κ := (κ 1 , e i η 2 , . . . , e i η d ) ∈ (S 1 ) d , we have found a root z ∈ U to the relation

s+1 σ=1 σ Q σ (κ) z σ-1 = 0 . ( 43
)
This is where the new form of Assumption 2 matters. Namely, we know that for all κ ∈ (S 1 ) d , the roots of the polynomial equation (4) lie in D and if there are roots on the boundary S 1 , then they must necessarily be simple. Applying again the Gauss-Lucas Theorem, we know that the roots to (43) lie in the convex hull of those to (4) and therefore belong to D (because the only possibility for (43) to have a root on the boundary S 1 would be that (4) admits a double root on S 1 but this degeneracy is precluded by Assumption 2). The Gauss-Lucas Theorem thus shows that the roots to the relation (43) do not belong to U . Hence M(z, η) has no eigenvalue on S 1 for any (z, η) ∈ U × R d-1 .

• At this stage, we know that for (z, η) ∈ U × R d-1 , the eigenvalues of M(z, η) split into two groups: those in U , which we call the unstable ones, and those in D, which we call the stable ones. For (z, η) ∈

End of the proof

The end of the proof of Theorem 1 follows, as in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF], from a superposition argument, see [BGS07, chapter 4] for a similar argument in the context of continuous problems. The solution to (1) with nonzero initial data is decomposed as the sum of a solution to an auxiliary problem (28) (that auxiliary problem incorporates the initial data) and of a solution to a problem of the form (1) with zero initial data (hence our earlier treatment of that case). The analysis in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] can be applied again word for word so we feel free to refer the reader to that earlier work.

A Numerical schemes with two time levels

As we have seen in the proof of Proposition 2, the construction of energy and dissipation functionals for finite difference operators is dictated, through the Plancherel Theorem, by the analogous construction for recurrence relations. The inconvenience in the proof of Lemma 1 is that the construction of the forms q e and q d depends on whether the roots of the polynomial P are simple. There is however one case that can be dealt with in a unified way and for which the coefficients of the forms q e and q d depend in a very simple and explicit way on the coefficients of P . Namely, we have the following result in the case of degree two polynomials2 (the case of degree one polynomials is actually even simpler).

Lemma 4 (The energy-dissipation balance law for second order recurrence relations). Let P = a X 2 + b X + c ∈ C[X] be a polynomial of degree 2 (a = 0), that satisfies the following two properties:

• The roots of P are located in D.

• If P has a double root, then it is located in D.

Then the Hermitian form q e , resp. q d , defined on C 2 by: ∀ (x 1 , x 2 ) ∈ C 2 , q e (x 1 , x ∀ n ∈ N , 2 Re T (P ′ (T) v n ) P (T) v n = 2 |P (T) v n | 2 + q e (v n+1 , v n+2 ) -q e (v n , v n+1 ) + q d (v n , v n+1 ) .

(44)

The defining equations for q e and q d in Lemma 4 show that, if P is a polynomial whose coefficients are trigonometric polynomials on R d , then the coefficients of q e and q d can also be chosen as trigonometric polynomials on R d (this was not the case, in general, with our construction from Lemma 1).

Proof. The validity of (44) is a mere algebra exercise. One can for instance expand the left hand side, which reads: 2 Re (2

a v n+2 + b v n+1 ) (a v n+2 + b v n+1 + c v n ) ,
and verify that it coincides with the right hand side of (44) (a good starting point for this calculation is first to subtract 2 |P (T) v n | 2 to the latter quantity and factorize P (T) v n within the real part before

  2 ) := 2 |a| 2 |x 2 | 2 + 2 Re a x 2 b x 1 + |a| 2 + |c| 2 |x 1 | 2 , q d (x 1 , x 2 ) := |a| 2 -|c| 2 |x 2 | 2 + 2 Re a x 2 b x 1 -2 Re b x 2 c x 1 + |a| 2 -|c| 2 |x 1 | 2 ,is positive definite, resp. nonnegative. Furthermore, for any sequence (v n ) n∈N with values in C, there holds:

From Assumption 1, we know that Qs+1 is an isomorphism on ℓ

(Z d ), which implies by Fourier analysis that Qs+1(κ) does not vanish for κ ∈ (S 1 ) d . In particular, the dispersion relation (4) is a polynomial equation of degree s + 1 in z for any κ ∈ (S 1 ) d .

Our attempts to obtain an analogue of Lemma 4 with 'explicit' Hermitian forms for degree three polynomials have been unsuccessful so far, not mentioning higher degrees.

was supported by ANR project Nabuco, ANR-17-CE40-0025

We use below the convention Q k,µ k := P k , which is compatible with the above definition of P k and of the Q k,j 's. Observing that there holds:

we have for any k = 1, . . . , m:

We use Young's inequality as follows:

and thus derive the lower bound:

Shifting indices, we get:

Restricting from now on to 0 < ε ≤ 1/4, we have 1/(2 ε) ≥ 2 ≥ 1 + |z k | 2 and all terms corresponding to the indices j = 1, . . . , µ k -2 in the above two sums match to give a nonnegative quantity (the first one U ×R d-1 , we then introduce the spectral projector Π s M (z, η), resp. Π u M (z, η), of M(z, η) on the generalized eigenspace associated with eigenvalues in D, resp. U . These projectors are analytic with respect to (z, η) on U × R d-1 . We can integrate from -∞ to 0 the recurrence relation (41) and get

The projector Π s M depends analytically on (z, η) ∈ U × R d-1 . Furthermore, since the spectrum of M does not meet S 1 for (z, η) ∈ U × R d-1 , there exists a constant C > 0 and a parameter δ ∈ (0, 1) that are independent of k ∈ N and such that

We thus get a uniform estimate with respect to k:

Passing to the limit and using (36), we get Π s M (z, η) W 0 = 0, or in other words W 0 = Π u M (z, η) W 0 . Furthermore, since (W j 1 ) j 1 ≤0 satisfies the recurrence relation (42) with W 0 in the generalized eigenspace of M(z, η) associated with eigenvalues in U , we find that (W j 1 ) j 1 ≤0 decays exponentially at -∞ and thus belongs to ℓ 2 (-∞, 0).

• The sequence (W j 1 ) j 1 ≤0 satisfies both recurrence relations (40) and ( 42), which equivalently means that the complex valued sequence (w j 1 ) j 1 ≤0 satisfies the two recurrence relations (34) and (35) for j 1 ≤ 0. Hence (w j 1 ) j 1 ≤0 satisfies the recurrence relation associated with the greatest common divisor of the polynomials associated with (34) and (35). In other words, the vector W 0 belongs to the generalized eigenspace (of either L or M) associated with the common eigenvalues of M(z, η) and L(z, η). Since we already know that M(z, η) has no eigenvalue on S 1 and that W 0 belongs to the generalized eigenspace of M(z, η) associated with eigenvalues in U (the unstable ones), we can conclude that W 0 also belongs to the generalized eigenspace of L(z, η) associated with those common eigenvalues of M(z, η) and L(z, η) in U .

The final argument is the following. The matrix L(z, η) has N u eigenvalues in U , N s in D and N c on S 1 (all eigenvalues are counted with multiplicity). (Since z may belong to S 1 , N c is not necessarily zero.) With rather obvious notations, we let Π u,s,c L (z, η) denote the corresponding spectral projectors of L for (z, η) sufficiently close to (z, η). In particular, the N u eigenvalues corresponding to Π u L (z, η) lie in U uniformly away from S 1 for (z, η) sufficiently close to (z, η). We can then integrate (39) from +∞ to 0 and derive (for k sufficiently large):

Using the uniform exponential decay of L(z k , η k ) -j 1 -1 Π u L (z k , η k ) (with respect to j 1 ) and the convergence (36), we finally end up with Π u L (z, η) W 0 = 0 . Since W 0 belongs to the generalized eigenspace of L associated with those common eigenvalues of M(z, η) and L(z, η) in U , we can conclude that W 0 equals zero. Applying the recurrence relation (40), the whole sequence (W j 1 ) j 1 ∈Z is zero, which yields the expected contradiction. expanding). The relation (44) can be also derived by noting that the above forms q e and q d in Lemma 4 differ from those given in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] (which we have recalled in the proof of Lemma 1) in terms of the Lagrange polynomials associated with P by the standard telescopic "add and subtract" trick. Namely, if z 1 , z 2 denote the two roots of P , then q e equivalently reads:

where the two first terms in the sum on the right hand side correspond to the Lagrange polynomials P 1 (T)x 1 and P 2 (T)x 1 (see ( 14)), and the very last term in the sum on the right hand side has been added in order to keep q e positive definite in case z 1 and z 2 coincide, in which case they belong to D (this last term was absent in [START_REF] Coulombel | The Leray-Gårding method for finite difference schemes[END_REF] since the roots were assumed to be simple). The link with the defining equation for q e in Lemma 4 is made by using the relations:

It is clear from the above alternative definition (45) that q e is positive definite under the assumptions we have made for the polynomial P .

Let us now turn to the dissipation form q d . In agreement with the alternative expression (45) for q e , the reader can check that the form q d given in Lemma 4 can be alternatively defined by the expression:

where the two first terms in the sum on the right hand side read as in (15), and the very last term on the right hand side has been added in order to keep the balance law (44) valid. Expanding the square moduli in the expression of q d , we find that it can be represented by the Hermitian matrix:

whose trace is clearly nonnegative since z 1 and z 2 belong to D. Furthermore, up to the positive |a| 4 factor, its determinant equals:

Expanding the square modulus and factorizing, the latter quantity is found to be equivalently given by:

which is bounded from below by the nonnegative quantity:

Hence the determinant of q d is nonnegative, so q d is nonnegative. The proof of Lemma 4 is complete.

Lemma 4 has an important consequence for the Cauchy problem (8) with s = 1 (finite difference operators with two time levels). Namely, if we follow the proof of Proposition 2 with the aim of constructing some energy and dissipation functionals for (8), we introduce the multiplier M as in (20) and obtain the relation (21). In the case s = 1, the polynomial P κ reads:

If the Cauchy problem (8) is ℓ 2 -stable, then the polynomial P κ satisfies the conditions of Lemma 4 for any κ ∈ (S 1 ) d . Hence we can apply Lemma 4 and rewrite (21) as:

with (here we apply the Plancherel Theorem 'backwards'):

and, similarly:

The interesting feature of these expressions is that both E and D correspond to the sum, with respect to j ∈ Z d , of local energy and dissipation densities E j (v n , v n+1 ), resp. D j (v n , v n+1 ), which depend on finitely many values of the sequences v n , v n+1 near j. For instance, the local density E j (v n , v n+1 ) can be defined by:

Hence there is now a genuine hope to extend the definition of E and D to more general domains (by means of sums of local quantities which do not rely on the Fourier transform) and/or to take the energydissipation balance law (46) as a starting point for deriving stability estimates for finite volume space discretizations on unstructured meshes. This is left to a future work.