\documentclass[JEP,XML,SOM]{cedram}
\datereceived{2015-06-12}
\dateaccepted{2015-11-06 }
\dateepreuves{2015-11-19}

\usepackage{mathrsfs}
\let\mathcal\mathscr
\multlinegap0pt
\let\oldforall\forall
\renewcommand{\forall}{\oldforall\,}
\newcommand{\bS}{\boldsymbol{S}}
\newcommand{\bT}{\boldsymbol{T}}
\newcommand{\rmd}{\mathrm{d}}
\newcommand{\rme}{\mathrm{e}}
\newcommand{\sbullet}{{\scriptscriptstyle\bullet}}
\DeclareMathOperator{\reel}{Re}

\theoremstyle{plain}
\newtheorem{lemma}{Lemma}
\newtheorem{proposition}{Proposition}
\newtheorem{theorem}{Theorem}
\newtheorem{corollary}{Corollary}
\theoremstyle{remark}
\newtheorem{definition}{Definition}
\newtheorem{assumption}{Assumption}
\newtheorem{remark}{Remark}

\newcommand{\N}{{\mathbb N}}
\newcommand{\Z}{{\mathbb Z}}
\newcommand{\R}{{\mathbb R}}
\newcommand{\C}{{\mathbb C}}
\newcommand{\T}{{\mathbb T}}

\newcommand{\A}{{\mathbb A}}
\newcommand{\B}{{\mathbb B}}
\newcommand{\E}{{\mathbb E}}
\newcommand{\F}{{\mathbb F}}
\newcommand{\LL}{{\mathbb L}}
\newcommand{\M}{{\mathbb M}}
\newcommand{\D}{{\mathbb D}}
\newcommand{\Dbar}{\overline{\mathbb D}}
\newcommand{\U}{{\mathcal U}}
\newcommand{\Ubar}{\overline{\mathcal U}}
\newcommand{\cercle}{{\mathbb S}^1}
\let\eps\varepsilon
\let\dps\displaystyle
\newcommand{\Ng}{| \! | \! |}
\newcommand{\Nd}{| \! | \! |}

\begin{document}
\frontmatter
\title{The Leray-G{\aa}rding method for finite~difference schemes}

\author[\initial{J.-F.} \lastname{Coulombel}]{\firstname{Jean-François} \lastname{Coulombel}}
\address{CNRS \& Université de Nantes, Laboratoire de Mathématiques
Jean Leray (UMR CNRS 6629)\\
2 rue de la Houssinière, BP 92208, 44322 Nantes Cedex 3, France}
\email{jean-francois.coulombel@univ-nantes.fr}
\urladdr{http://www.math.sciences.univ-nantes.fr/~coulombel/}

\thanks{The research of the author was supported by the ANR project BoND, ANR-13-BS01-0009-01.}

\begin{abstract}
In the fifties, Leray and G{\aa}rding have developed a multiplier technique for deriving a priori estimates for solutions to scalar hyperbolic equations. The existence of such a multiplier is the starting point of the argument by Rauch \cite{rauch} for the derivation of semigroup estimates for hyperbolic initial boundary value problems. In this article, we explain how this multiplier technique can be adapted to the framework of finite difference approximations of transport equations. The technique applies to numerical schemes with arbitrarily many time levels. The existence and properties of the multiplier enable us to derive {\it optimal} semigroup estimates for fully discrete hyperbolic initial boundary value problems.
\end{abstract}

\subjclass{65M06, 65M12, 35L03, 35L04}

\keywords{Hyperbolic equations, difference approximations, stability, boundary conditions,
semigroup}

\alttitle{La méthode de Leray et G{\aa}rding pour les schémas aux  différences finies}

\begin{altabstract}
Dans les années 1950, Leray et G{\aa}rding ont développé une  technique de multiplicateur pour obtenir des estimations a priori de solutions d'équations  hyperboliques scalaires. L'existence d'un multiplicateur est le point de départ du travail de  Rauch \cite{rauch} pour montrer des estimations de semi-groupe pour les problèmes aux limites  hyperboliques. Dans cet article, nous expliquons comment cette technique de multiplicateur  peut être adaptée au cadre des schémas aux différences finies pour les équations de  transport. Ce travail s'applique à des schémas numériques multi-pas en temps. L'existence  et les propriétés du multiplicateur nous permettent d'obtenir des estimations de  semi-groupe optimales pour des versions totalement discrètes des problèmes aux limites hyperboliques.
\end{altabstract}

\altkeywords{Équations hyperboliques, différences finies, stabilité, conditions  aux limites, semi-groupe}

\maketitle
\tableofcontents
\mainmatter

\newpage
Throughout this article, we use the notation
\begin{align*}
\U &:= \{\zeta \in \C,|\zeta|>1 \},& \Ubar &:= \{\zeta \in \C,|\zeta| \ge 1 \},\\
\D &:= \{\zeta \in \C,|\zeta|<1 \},& \cercle &:= \{\zeta \in \C,|\zeta|=1 \},& \Dbar &:= \D \cup \cercle.
\end{align*}
We let ${\mathcal M}_n ({\mathbb K})$ denote the set of $n \times n$ matrices with entries in ${\mathbb K}
= \R \text{ or } \C$. If $M \in {\mathcal 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 \in \C^n$ and for the Euclidean product
of two vectors $x,y \in \R^n$. The norm of a vector $x \in \C^n$ is $|x| := (x^* x)^{1/2}$. The induced matrix
norm on ${\mathcal M}_n (\C)$ is denoted $\| \cdot \|$.

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.

In what follows, we let $d \ge 1$ denote a fixed integer, which will stand for the dimension of the space domain
we are considering. We also use the space $\ell^2$ of square integrable sequences. Sequences may be valued
in $\C^k$ for some integer $k$. Some sequences will be indexed by $\Z^{d-1}$ while some will be indexed by
$\Z^d$ or a subset of $\Z^d$. We thus introduce some specific notation for the norms. Let $\Delta x_i>0$ for
$i=1,\dots,d$ be $d$ space steps. We shall make use of the $\ell^2(\Z^{d-1})$-norm that we define as follows:
for all $v \in \ell^2 (\Z^{d-1})$,
\begin{equation*}
\| v \|_{\ell^2(\Z^{d-1})}^2 := \biggl(\prod_{k=2}^d \Delta x_k \biggr) \sum_{i=2}^d \sum_{j_i\in \Z} |v_{j_2,\dots,j_d}|^2.
\end{equation*}
The corresponding scalar product is denoted $\langle \cdot\,,\cdot \rangle_{\ell^2(\Z^{d-1})}$. Then for all integers
$m_1 \le\nobreak m_2$, we set
\begin{equation*}
\Ng u \Nd_{m_1,m_2}^2 := \Delta x_1 \sum_{j_1=m_1}^{m_2} \|u_{j_1,\sbullet} \|_{\ell^2(\Z^{d-1})}^2,
\end{equation*}
to denote the $\ell^2$-norm on the set $[m_1,m_2] \times \Z^{d-1}$ ($m_1$ may equal $-\infty$ and $m_2$
may equal $+\infty$). The corresponding scalar product is denoted $\langle \cdot\,,\cdot \rangle_{m_1,m_2}$.
Other notation throughout the text is meant to be self-explanatory.

\section{Introduction}
\label{intro}

\subsection{The context}

The goal of this article is to derive semigroup estimates for finite difference approximations of hyperbolic initial
boundary value problems. Up to now, the only available general stability theory for such numerical schemes is
due to Gustafsson, Kreiss and Sundström \cite{gks}. It relies on a Laplace transform with
respect to the time variable. For various technical reasons, the corresponding stability estimates are restricted
to zero initial data. A long standing problem in this line of research~is, starting from the GKS stability estimates,
which are {\it resolvent} type estimates, to incorporate nonzero initial data and to derive {\it semigroup} estimates,
see, e.g., the discussion by Trefethen in \cite[\S4]{trefethen3} and the conjecture by Kreiss and
Wu in~\cite{kreiss-wu}. This problem is delicate for the following reason: the validity of the GKS \hbox{stability}
estimate is known to be equivalent to a {\it slightly stronger version} of the resolvent estimate\vspace*{-5pt}
\begin{equation}
\label{resolventkreiss}
\sup_{z \in \U} (|z|-1) \bigl\| (z I-T)^{-1} \bigr\|_{{\mathcal L}(\ell^2(\N))} <+\infty,
\end{equation}
where $T$ is some bounded operator on $\ell^2(\N)$ that incorporates both the discretization of the hyperbolic
equation and the numerical boundary conditions. Deriving an optimal semigroup estimate amounts to showing
that $T$ is power bounded:
\begin{equation*}
\sup_{n \ge 1} \| T^n \|_{{\mathcal L}(\ell^2(\N))} <+\infty.
\end{equation*}
In finite dimension, the equivalence between power boundedness of $T$ and the resolvent condition
\eqref{resolventkreiss} is known as the Kreiss matrix Theorem, but the analogous equivalence is known to
fail in general in infinite dimension. Worse, even the strong resolvent condition\vspace*{-5pt}
\begin{equation*}
\sup_{n \ge 1}\; \sup_{z \in \U} (|z|-1)^n \bigl\| (z I-T)^{-n} \bigr\|_{{\mathcal L}(\ell^2(\N))} <+\infty,
\end{equation*}
does not imply in general that $T$ is power bounded, see, e.g., the review \cite{strikwerda-wade} or
\cite[Chap.\,18]{TE} for details and historical comments.

Optimal semigroup estimates have nevertheless been derived for some discretized hyperbolic initial boundary
value problems. The very first results in this direction date back to Kreiss and Osher
\cite{kreiss1,osher1,osher2}, even though these works precede \cite{gks} but the main results are exactly of
the form we discuss. The first general derivation of semigroup estimates starting from GKS stability is due to
Wu \cite{wu}, whose analysis deals with numerical schemes with two time levels and scalar equations
(as in \cite{kreiss1,osher1,osher2}). The results in \cite{wu} were extended by
Gloria and the author in
\cite{jfcag} to {\it systems} in arbitrary space dimension, but the arguments in \cite{jfcag} are still restricted to
numerical schemes with two time levels. The present article gives, as far as we are aware of, the first systematic
derivation of semigroup estimates for fully discrete hyperbolic initial boundary value problems in the case of
numerical schemes with arbitrarily many time levels. It generalizes the arguments of \cite{wu,jfcag} and provides
new insight for the construction of `dissipative' (sometimes called `absorbing') numerical boundary conditions for
discretized evolution equations. Let us observe that the leap-frog scheme, with some very specific homogeneous
boundary conditions, has been dealt with by Thomas \cite{thomas} by using a {\it multiplier} technique. It is
precisely this technique which we aim at developing in a systematic fashion for numerical schemes with arbitrarily
many time levels. In particular, we shall explain why the somehow magical multiplier $u_j^{n+2}+u_j^n$ for the
leap-frog scheme, see, e.g., \cite{RM}, follows from a general theory that is the analogue of the
Leray--G{\aa}rding method \cite{leray,garding} for hyperbolic partial differential equations.

\subsection{The main result}

We first set a few notations. We let \hbox{$\Delta x_1,\dots,\Delta x_d,\Delta t>0$} denote space and time steps
where the ratios, the so-called Courant-Friedrichs-Lewy parameters, $\lambda_i :=\Delta t/\Delta x_i$,
$i=1,\dots,d$, are fixed positive constants. We~keep $\Delta t \in\nobreak (0,1]$ as a small parameter and let the
space steps $\Delta x_1,\dots,\Delta x_d$ vary accordingly. The $\ell^2$-norms with respect to the space
variables have been previously defined and thus depend on $\Delta t$ and the CFL parameters through
the cell volume ($\Delta x_2 \cdots \Delta x_d$ on~$\Z^{d-1}$, and $\Delta x_1 \cdots \Delta 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. For instance, for a sequence $(w^n)_{n \in \N}$, the step
function $w$ reads:
\begin{equation*}
w(t) :=w^n,\quad \forall t \in [n \Delta t,(n+1) \Delta t),\quad \forall n \in \N.
\end{equation*}
In particular, we shall feel free to take Fourier and/or Laplace transforms of sequences.

For all $j\in \Z^d$, we set $j=(j_1,j')$ with $j':=(j_2,\dots,j_d) \in \Z^{d-1}$. We let $p,q,r \in \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 \in \N$
denote some fixed integer. We consider a recurrence relation of the form:
\begin{equation}
\label{numibvp}
\begin{cases}
{\dps \sum_{\sigma=0}^{s+1}} Q_\sigma u_j^{n+\sigma} =\Delta t F_j^{n+s+1},&
j \in \Z^d,\quad j_1 \ge 1,\quad n\ge 0,\\[-3pt]
u_j^{n+s+1} +{\dps \sum_{\sigma=0}^{s+1}} B_{j_1,\sigma} u_{1,j'}^{n+\sigma} =g_j^{n+s+1},&
j \in \Z^d,\quad j_1=1-r_1,\dots,0,\quad n\ge 0,\\[-3pt]
u_j^n = f_j^n,& j \in \Z^d,\quad j_1\ge 1-r_1,\quad n=0,\dots,s,
\end{cases}
\end{equation}
where the operators $Q_\sigma$ and $B_{j_1,\sigma}$ are given by:
\begin{equation}
\label{defop}
Q_\sigma := \sum_{\ell_1=-r_1}^{p_1} \sum_{\ell'=-r'}^{p'} a_{\ell,\sigma} \bS^\ell,\quad
B_{j_1,\sigma} := \sum_{\ell_1=0}^{q_1} \sum_{\ell'=-q'}^{q'} b_{\ell,j_1,\sigma} \bS^\ell.
\end{equation}
Let us comment a little on \eqref{numibvp}, \eqref{defop}. First of all, in \eqref{defop}, the coefficients
$a_{\ell,\sigma},b_{\ell,j_1,\sigma}$ are {\it real numbers} and are independent of the small parameter $\Delta t$
(they may depend on the CFL parameters though), while $\bS$ denotes the shift operator on the space grid:
$(\bS^\ell v)_j :=v_{j+\ell}$ for $j,\ell \in \Z^d$. We have also used the short notation
\begin{equation*}
\sum_{\ell'=-r'}^{p'} :=\sum_{i=2}^d \sum_{\ell_i=-r_i}^{p_i},\qquad
\sum_{\ell'=-q'}^{q'} :=\sum_{i=2}^d \sum_{\ell_i=-q_i}^{q_i}.
\end{equation*}
This means that each operator $Q_\sigma,B_{j_1,\sigma}$ in \eqref{numibvp} acts on sequences indexed by
the `spatial variable' $j \in \Z^d$. In particular, the `time variable' $n$ enters as a parameter when we apply these
operators and, for instance, $Q_\sigma u_j^{n+\sigma}$ is a short notation for the application of the operator
$Q_\sigma$ to the sequence $u^{n+\sigma}$, this resulting sequence being evaluated at the space index $j$.

The numerical scheme \eqref{numibvp} should then be understood as follows: one starts with~$\ell^2$ initial data
$(f_j^0),\dots, (f_j^s)$ defined for $j_1 \ge 1-r_1$, $j' \in \Z^{d-1}$, some given interior source term $(F_j^n)$
defined for $j_1 \ge 1$, $j' \in \Z^{d-1}$, $n \ge s+1$, and some given boundary source term $(g_j^n)$ defined for
$j_1=1-r_1,\dots,0$, $j' \in \Z^{d-1}$, $n \ge s+1$. The (space) cells associated with $j_1 \ge 1$ correspond to the
{\it interior domain} (the discretized counterpart of the half-space $\{ x_1 >0 \}$ in $\R^d$), while those associated
with $j_1=1-r_1,\dots,0$ represent the {\it discrete boundary} (the discretized counterpart of the hyperplane
$\{ x_1 =\nobreak0 \}$ in $\R^d$). Assuming that the solution $(u_j^n)$ has been defined up to some time index $n+s$ for
all $j_1 \ge 1-r_1$, $j' \in \Z^{d-1}$, with $n \ge 0$, then the first and second equations in \eqref{numibvp} should
uniquely determine $u_j^{n+s+1}$ for $j_1 \ge 1-r_1$, $j' \in \Z^{d-1}$. More precisely, the equation
\begin{equation*}
{\dps \sum_{\sigma=0}^{s+1}} Q_\sigma u_j^{n+\sigma} =\Delta t F_j^{n+s+1},
\end{equation*}
is meant to give an update for the interior values of $u^{n+s+1}$ by possibly inverting the operator $Q_{s+1}$, and
the numerical boundary conditions
\begin{equation*}
u_j^{n+s+1} +{\dps \sum_{\sigma=0}^{s+1}} B_{j_1,\sigma} u_{1,j'}^{n+\sigma} =g_j^{n+s+1}
\end{equation*}
should determine the boundary values $u_j^{n+s+1}$, $j_1=1-r_1,\dots,0$, by possibly giving their expression
in terms of finitely many interior values. More precisely, for each $j_1 =1-r_1,\dots,0$, and independently of the
tangential variable $j' \in \Z^{d-1}$, $u_j^{n+s+1}$ is assumed to be given by a linear combination of interior
values, which amounts to considering a linear combination of $u_{1,j'+\ell'}^{n+\sigma}, \dots,
u_{1+q_1,j'+\ell'}^{n+\sigma}$ with a finite stencil for $\ell'$ and $\sigma=0,\dots,s+1$. Some examples will
be given later on. Let us just emphasize here that some more general numerical boundary conditions might
probably be considered, including convolution type formula with respect to $n$, but we restrict in this article
to this form of numerical boundary conditions for simplicity.

We wish to deal here simultaneously with explicit and implicit schemes and therefore make the following solvability
assumption.

\begin{assumption}[Solvability of \eqref{numibvp}]
\label{assumption0}
The operator $Q_{s+1}$ is an isomorphism on $\ell^2 (\Z^d)$. Moreover, for all $(F_j) \in \ell^2 (\N^* \times \Z^{d-1})$
and for all $g_{1-r_1,\sbullet},\dots,g_{0,\sbullet} \in \ell^2 (\Z^{d-1})$, there exists a unique solution $(u_j)_{j_1 \ge 1-r_1}
\in \ell^2$ to the system
\begin{equation*}
\begin{cases}
Q_{s+1} u_j =F_j,& j \in \Z^d,\quad j_1 \ge 1,\\
u_j +B_{j_1,s+1} u_{1,j'} =g_j,& j \in \Z^d,\quad j_1=1-r_1,\dots,0.
\end{cases}
\end{equation*}
\end{assumption}

\noindent In particular, Assumption \ref{assumption0} is trivially satisfied in the case of explicit schemes for which
$Q_{s+1}$ is the identity ($a_{\ell,s+1} =\delta_{\ell_1,0} \cdots \delta_{\ell_d,0}$ in \eqref{defop}, with $\delta$ the
Kronecker symbol), for in that case, we first determine all interior values by the relation $u_j =F_j$ ($F_j$ is given),
and we then use these values to define the boundary values $u_j$, $j_1=1-r_1,\dots,0$. The situation is more or
less the same as for lower or upper triangular systems. When $Q_{s+1}$ is not the identity, then the values $u_j$,
$j_1 \ge 1-r_1$, should be determined all at once.

Using Assumption \ref{assumption0}, the first and second equations in \eqref{numibvp} uniquely determine
$u_j^{n+s+1}$ for $j_1 \ge 1-r_1$, and one then proceeds to the following time index $n+s+2$. Existence and
uniqueness of a solution $(u_j^n)$ to \eqref{numibvp} thus follows from Assumption \ref{assumption0}, so the
last requirement for well-posedness is continuous dependence of the solution on the three possible source terms
$(F_j^n)$, $(g_j^n)$, $(f_j^n)$. This is a {\it stability} problem for which several definitions can be chosen according
to the functional framework. The following one dates back to \cite{gks} in one space dimension and was also
considered by Michelson \cite{michelson} in several space dimensions. It is specifically relevant when the
boundary conditions are non-homogeneous ($(g_j^n) \not \equiv 0$).

\begin{definition}[Strong stability]
\label{defstab1}
The finite difference approximation \eqref{numibvp} is said to be `strongly stable' if there exists a constant
$C$ such that for all $\gamma>0$ and all $\Delta t \in (0,1]$, the solution $(u_j^n)$ to \eqref{numibvp}
with $(f_j^0) =\dots =(f_j^s) =0$ satisfies the estimate:
\begin{multline}
\label{stabilitenumibvp}
\dfrac{\gamma}{\gamma \Delta t+1} \sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t}
\Ng u^n \Nd_{1-r_1,+\infty}^2\\[-5pt]
+\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t}
\sum_{j=1-r_1}^{p_1} \| u_{j_1,\sbullet}^n \|_{\ell^2 (\Z^{d-1})}^2 \\
\le C \biggl\{ \dfrac{\gamma \Delta t+1}{\gamma}\hspace*{-1mm} \sum_{n\ge s+1}\hspace*{-1mm} \Delta t
\rme^{-2 \gamma n \Delta t} \Ng F^n \Nd_{1,+\infty}^2
+\sum_{n\ge s+1}\hspace*{-1mm} \Delta t \rme^{-2 \gamma n \Delta t} \hspace*{-1mm} \sum_{j_1=1-r_1}^0
\| g_{j_1,\sbullet}^n \|_{\ell^2 (\Z^{d-1})}^2 \biggr\}.
\end{multline}
\end{definition}

The main contributions in \cite{gks,michelson} are to show that strong stability can be characterized by a
certain {\it algebraic condition}, which is usually referred to as the Uniform Kreiss-Lopatinskii Condition,
see \cite{jfcnotes} for an overview of such results. We do not pursue such arguments here but rather assume
from the start that \eqref{numibvp} is strongly stable. We can thus control, with zero initial data, $\ell^2$ type
norms of the solution to \eqref{numibvp}. Our goal is to understand which kind of stability estimate holds for
the solution to \eqref{numibvp} when one now considers nonzero initial data $(f_j^0),\dots,(f_j^s)$ in $\ell^2$.
Our main assumption is the following and is related to stability of the recurrence relation:
\begin{equation*}
{\dps \sum_{\sigma=0}^{s+1}} Q_\sigma u_j^{n+\sigma} =0,
\end{equation*}
when applied on all $\Z^d$. Recall that in that case, stability is usually understood in the $\ell^2$ sense,
and by Fourier analysis, this leads to studying solutions to the latter recurrence relation of the form $z^n
\exp (i j \cdot \xi)$ for arbitrary wave vectors $\xi \in \R^d$, hence the so-called dispersion relation
\eqref{dispersion} below, see, e.g., \cite{RM,gko}.

\begin{assumption}[Stability for the discrete Cauchy problem]
\label{assumption1}
For all $\xi \in \R^d$, the dispersion relation
\begin{equation}
\label{dispersion}
\sum_{\sigma=0}^{s+1} \widehat{Q_\sigma} (\rme^{i \xi_1},\dots,\rme^{i \xi_d}) z^\sigma =0,\quad
\widehat{Q_\sigma}(\kappa) := \sum_{\ell=-r}^p \kappa^\ell a_{\ell,\sigma},
\end{equation}
has $s+1$ simple roots in $\Dbar$. (The von Neumann condition is said to hold when the roots are located
in $\Dbar$.) In \eqref{dispersion}, we have used the classical notation $\kappa^\ell :=\kappa_1^{\ell_1} \cdots
\kappa_d^{\ell_d}$ for $\kappa \in (\C \setminus \{ 0 \})^d$ and $\ell \in \Z^d$.
\end{assumption}

\noindent Examples of numerical schemes that satisfy Assumption \ref{assumption1} are given in Section
\ref{examples}. At the opposite, Assumption \ref{assumption1} excludes numerical schemes that are based
on first performing a space discretization and then using Adams-Bashforth or Adams-Moulton methods of
order $3$ or higher (such methods have $0$ as a root of multiplicity~$2$ or more at the frequency $\xi=0$,
see \cite[Chap.\,III.3]{hnw}).

From Assumption \ref{assumption0}, we know that $Q_{s+1}$ is an isomorphism on $\ell^2$, which
implies by Fourier analysis that $\widehat{Q_{s+1}} (\rme^{i \xi_1},\dots,\rme^{i \xi_d})$ does not vanish
for any $\xi \in \R^d$. In particular, the dispersion relation \eqref{dispersion} is a polynomial equation of degree
$s+1$ in $z$ for any $\xi \in \R^d$. We now make the following assumption, which already appeared in
\cite{gks,michelson} and several other works on the same topic.

\begin{assumption}[Noncharacteristic discrete boundary]
\label{assumption2}
For $\ell_1=-r_1,\dots,p_1$, $z \in \C$ and $\eta \in \R^{d-1}$, let us define
\begin{equation}
\label{defA-d}
a_{\ell_1}(z,\eta) :=\sum_{\sigma=0}^{s+1} z^\sigma \sum_{\ell'=-r'}^{p'}
a_{\ell,\sigma} \rme^{i \ell' \cdot \eta}.
\end{equation}
Then $a_{-r_1}$ and $a_{p_1}$ do not vanish on $\Ubar \times \R^{d-1}$, and they have nonzero degree
with respect to $z$ for all $\eta \in \R^{d-1}$.
\end{assumption}

\noindent Assumption \ref{assumption2} is used when one performs a Laplace-Fourier transform on
\eqref{numibvp}. The Laplace transform refers to the time variable $n \in \N$ and the Fourier transform refers
to the tangential space variables $j' \in \Z^{d-1}$. One is then led to a recurrence relation with respect to the
space normal variable $j_1$ which, thanks to Assumption \ref{assumption2} can be either written as an `evolution'
equation for increasing or decreasing $j_1$. This will be used in Section \ref{section3}.

Our main result is comparable with \cite[Th.\,3.3]{wu} and \cite[Th.\,2.4 \& 3.5]{jfcag} and shows that
strong stability (or GKS stability) is a sufficient condition for incorporating $\ell^2$ initial conditions in \eqref{numibvp}
and proving {\it optimal} semigroup estimates. The main price to pay in Assumption \ref{assumption1} is that
the roots of the dispersion relation \eqref{dispersion}, which are nothing but the eigenvalues of the so-called
{\it amplification matrix} for the Cauchy problem, need to be {\it simple}. This property is satisfied by the leap-frog
and modified leap-frog schemes in two space dimensions under an appropriate CFL condition, see Section
\ref{examples} for more examples. Our main result reads as follows.

\begin{theorem}
\label{mainthm}
Let Assumptions \ref{assumption0}, \ref{assumption1} and \ref{assumption2} be satisfied, and assume
that the scheme \eqref{numibvp} is strongly stable in the sense of Definition \ref{defstab1}. Then there exists
a constant $C$ such that for all $\gamma>0$ and all $\Delta t \in (0,1]$, the solution to \eqref{numibvp}
satisfies the estimate:
\begin{multline}
\label{estim1d}
\sup_{n \ge 0} \rme^{-2 \gamma n \Delta t} \Ng u^n \Nd_{1-r_1,+\infty}^2+\dfrac{\gamma}{\gamma \Delta t+1}
\sum_{n\ge 0} \Delta t \rme^{-2 \gamma n \Delta t} \Ng u^n \Nd_{1-r_1,+\infty}^2 \\[-5pt]
\shoveright{+\sum_{n\ge 0} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1=1-r_1}^{p_1}\| u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2}
\\
\le C \biggl\{ \sum_{\sigma=0}^s \Ng f^\sigma \Nd_{1-r_1,+\infty}^2
+\dfrac{\gamma \Delta t+1}{\gamma}
\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \Ng F^n \Nd_{1,+\infty}^2\\[-5pt]
+\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1=1-r_1}^0
\| g_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \biggr\}.
\end{multline}
In particular, the scheme \eqref{numibvp} is `semigroup stable' in the sense that there exists a constant
$C$ such that for all $\Delta t \in (0,1]$, the solution $(u_j^n)$ to \eqref{numibvp} with $(F_j^n) =(g_j^n)
=0$ satisfies the estimate
\begin{equation}
\label{estimsemigroup}
\sup_{n\ge 0} \Ng u^n \Nd_{1-r_1,+\infty}^2 \le C \sum_{\sigma=0}^s \Ng f^\sigma \Nd_{1-r_1,+\infty}^2.
\end{equation}
The scheme \eqref{numibvp} is also $\ell^2$-stable with respect to boundary data, see
\cite[Def.\,4.5]{trefethen3}, in the sense that there exists a constant $C$ such that for all $\Delta t \in
(0,1]$, the solution $(u_j^n)$ to \eqref{numibvp} with $(F_j^n)=(f_j^n)=0$ satisfies the estimate
\begin{equation*}
\sup_{n\ge 0} \Ng u^n \Nd_{1-r_1,+\infty}^2 \le C \sum_{n\ge s+1} \Delta t \sum_{j_1=1-r_1}^0
\| g_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2.
\end{equation*}
\end{theorem}

\noindent The semigroup estimate \eqref{estimsemigroup} as well as $\ell^2$-stability with respect to boundary
data are indeed trivial consequences of the main estimate \eqref{estim1d} by letting the parameter $\gamma$
tend to zero. Our main contribution in this article is to exhibit a suitable multiplier for the multistep recurrence
relation in \eqref{numibvp}. With this multiplier, we can readily show that, for zero initial data, the (discrete)
derivative of an {\it energy} can be controlled, as in the work by Rauch \cite{rauch} on partial differential
equations, by the trace estimate of $(u_j^n)$ and this is where strong stability comes into play. This first argument
gives Theorem \ref{mainthm} for zero initial data\footnote{It would even give the claim of Theorem \ref{mainthm}
for nonzero initial data provided that the non-glancing condition of \cite{jfc} is satisfied, but we do not wish to make
such an assumption here.}. By linearity we can then reduce to the case of zero forcing terms in the interior and on
the boundary. The next arguments in \cite{rauch} use time reversibility, which basically always fails for numerical
schemes\footnote{With the notable exception of the leap-frog scheme that is time reversible! Some schemes
based on the Crank-Nicolson integration rule are also time reversible.}. Hence we must find another argument
for dealing with nonzero initial data. Hopefully, the properties of our multiplier enable us to construct an auxiliary
problem, where we modify the boundary conditions of \eqref{numibvp}, and for which we can prove optimal
semigroup and trace estimates by `hand-made' calculations. In other words, we exhibit an alternative set of
boundary conditions that yields {\it strict dissipativity}. Using these auxiliary numerical boundary conditions,
the proof of Theorem \ref{mainthm} follows from an easy (though lengthy) superposition argument, see, e.g.,
\cite[\S4.5]{benzoni-serre} for partial differential equations or \cite{wu,jfcag} for numerical schemes.

\begin{remark}
It could seem at first that the general form of \eqref{numibvp} incorporates not only finite difference approximations
of hyperbolic equations but more generally finite difference approximations of any evolutionary constant coefficient
partial differential equation. Hence Theorem \ref{mainthm} could apply to more general situations. However, it
should be kept in mind that we assume here that each ratio $\Delta t/\Delta x_j$ is constant, and therefore we
consider each coefficient in the operators $Q_\sigma,B_{j_1,\sigma}$ as independent of the time and space steps.
This point of view has some technical advantages since we may for instance view the $Q_\sigma$'s as bounded
operators with norms that are independent of the time and space steps, and all estimates in Theorem \ref{mainthm}
are in fact independent of the (only left) small parameter $\Delta t$ (just divide for instance \eqref{estimsemigroup}
by $\Delta t^d$ and use the definitions of the norms on either side to simplify all cell volumes). However, our
assumption is also a clear limitation when the original system is stiff and one uses implicit schemes in order to
get free of CFL constraints (just think of an implicit discretization of the heat equation). In that case, there would
be several small parameters involved and the coefficients in $Q_\sigma,B_{j_1,\sigma}$ could not all necessarily
be considered as constants (or even bounded). We postpone the extension of this work to parabolic or dispersive
equations to a future work.
\end{remark}

\Subsection{Examples}
\label{examples}

\subsubsection{Examples in one space dimension}

Our goal is to approximate the outgoing transport equation ($d=1$ here):
\begin{equation}
\label{transport}
\partial_t u +a \partial_x u=0,\quad u|_{t=0} =u_0,
\end{equation}
with $t,x>0$ and $a<0$. The latter transport equation does not require any boundary condition at $x=0$.
However, discretizing \eqref{transport} usually requires prescribing numerical boundary conditions, unless
one considers an upwind type scheme with a space stencil `on the right' (meaning $r_1=0$ in \eqref{numibvp}).

Let us first emphasize that Assumption \ref{assumption2} excludes the case of explicit two level schemes for which
$s=0$ and $Q_1=I$, for in that case $a_{-r_1}$ and/or $a_{p_1}$ do not depend on $z$. However, this case has
already been dealt with in \cite{wu,jfcag}, and we shall see in Section \ref{section3} where the assumption that
$a_{-r_1}$ and $a_{p_1}$ are {\it not constant} is involved, and why the proof is actually simpler in the case
$s=0$ and $Q_1=I$.

We now detail two possible multistep schemes for discretizing \eqref{transport}. Both are obtained by the
so-called method of lines, which amounts to first discretizing the space derivative $\partial_x u$ and then
choosing an integration technique for discretizing the time evolution, see \cite{gko}.

\subsubsection*{The leap-frog scheme.} It is obtained by approximating the space derivative $\partial_x u$
by the centered difference $(u_{j+1}-u_{j-1})/(2 \Delta x)$, and by then applying the so-called Nyström
method of order $2$, see \cite[Chap.\,III.1]{hnw}. The resulting approximation reads
\begin{equation*}
u_j^{n+2} +\lambda a (u_{j+1}^{n+1}-u_{j-1}^{n+1}) -u_j^n =0,
\end{equation*}
which corresponds to $s=p=r=1$ (here $d=1$ so we write $p$ instead of $p_1$ and so on). Recall that
$\lambda>0$ denotes the fixed ratio $\Delta t/\Delta x$. Even though \eqref{transport} does not require any
boundary condition at $x=0$, the leap-frog scheme stencil includes one point to the left, and we therefore
need to prescribe some numerical boundary condition at $j=0$. One possibility is to prescribe the homogeneous
or inhomogeneous Dirichlet boundary condition. With all possible source terms, the corresponding scheme reads
\begin{equation}
\label{leapfrog}
\begin{cases}
u_j^{n+2} +\lambda a (u_{j+1}^{n+1}-u_{j-1}^{n+1}) -u_j^n =\Delta t F_j^{n+2},&
j\ge 1,\quad n\ge 0,\\
u_0^{n+2} = g_0^{n+2},& n\ge 0,\\
(u_j^0,u_j^1) = (f_j^0,f_j^1),& j\ge 0.
\end{cases}
\end{equation}
Assumption \ref{assumption0} is trivially satisfied because \eqref{leapfrog} is explicit. More precisely, \eqref{leapfrog}
can be written under the form \eqref{numibvp} by setting:
\begin{equation*}
Q_2 := I,\quad Q_1 := \lambda a (\bS -\bS^{-1}),\quad Q_0 := -I,
\end{equation*}
and all operators $B_{j,\sigma}$ are zero ($q=0$ also). The leap-frog scheme satisfies Assumption~\ref{assumption1}
provided that $\lambda |a|<1$. In that case, the two roots to the dispersion relation
\begin{equation}
\label{dispersionlf}
z^2 +2 i \lambda a \sin \xi z -1 =0,
\end{equation}
are simple and have modulus $1$ for all $\xi \in \R$. Assumption \ref{assumption2} is satisfied as long as
the velocity $a$ is nonzero, for in that case $a_1(z)=-a_{-1}(z)=\lambda a z$. The scheme \eqref{leapfrog}
is known to be strongly stable, see \cite{goldberg-tadmor}. In particular, Theorem \ref{mainthm} shows that
\eqref{leapfrog} is semigroup stable. More accurate numerical boundary conditions can be considered, and
we refer to \cite{gko,oliger,sloan,trefethen3} for some other possible choices which might be more meaningful
from a consistency and accuracy point of view.

\subsubsection*{A scheme based on the backwards differentiation rule.} We still start from the transport equation
\eqref{transport}, approximate the space derivative $\partial_x u$ by the centered finite difference $(u_{j+1}-u_{j-1})
/(2 \Delta x)$, and then apply the backwards differentiation formula of order $2$, see \cite[Chap.\,III.1]{hnw}. The
resulting scheme reads:
\begin{equation}
\label{bdf2}
\dfrac{3}{2} u_j^{n+2} +\dfrac{\lambda a}{2} (u_{j+1}^{n+2}-u_{j-1}^{n+2}) -2 u_j^{n+1} +\dfrac{1}{2} u_j^n =0,
\end{equation}
which corresponds to $s=1$ and
\begin{equation*}
Q_2 := \dfrac{3}{2} I +\dfrac{\lambda a}{2} (\bS -\bS^{-1}),\quad
Q_1 := -2 I,\quad Q_0 := \dfrac{1}{2} I.
\end{equation*}
The operator $Q_2$ is an isomorphism on $\ell^2(\Z)$ since $Q_2$ is an isomorphism for any small $\lambda a$
(as a perturbation of $3/2 I$), $Q_2$ depends continuously on $\lambda a$, and there holds (uniformly with
respect to $\lambda a$):
\begin{equation*}
\dfrac{3}{2} \Ng u \Nd_{-\infty,+\infty} \le \Ng Q_2 u \Nd_{-\infty,+\infty}.
\end{equation*}
The operator $Q_2$ is therefore an isomorphism on $\ell^2(\Z)$ for any $\lambda a$ (see, e.g.,
\hbox{\cite[Lem.\,4.3]{jfcsinum}}). Let us now study the dispersion relation \eqref{dispersion}, which reads here
\begin{equation}
\label{dispersionbdf}
\Bigl(\dfrac{3}{2} +i \lambda a \sin \xi \Bigr) z^2 -2 z +\dfrac{1}{2} =0.
\end{equation}
It is clear that the latter equation has two simple roots in $z$ for any $\xi \in \R$. Moreover, if $\sin \xi=0$, the roots
are $1$ and $1/3$ which belong to $\Dbar$. In the case $\sin \xi \neq 0$, none of the roots belongs to $\cercle$ and
examining the case $\lambda a \sin \xi=1$, we find that for $\sin \xi \neq 0$, both roots belong to $\D$ (which is
consistent with the shape of the stability region for the backwards differentiation formula of order $2$, see
\cite[Chap.\,V.1]{hw}). Assumption \ref{assumption1} is therefore satisfied. Assumption \ref{assumption2} is satisfied
as long as $a$ is nonzero since there holds $p=r=1$ and $a_1(z)=-a_{-1}(z)=\lambda a z^2/2$.

Theorem \ref{mainthm} therefore yields semigroup boundedness as long as one uses numerical boundary conditions
for which the numerical scheme is well-defined (this is at least the case for $\lambda a$ small enough) and strong
stability holds. We shall go back later on to the form of our multiplier for the scheme \eqref{bdf2} and compare it with
another technique that is available in the literature, see, e.g., \cite{Emmrich1,Emmrich2} and references therein.

\subsubsection{Examples in two space dimensions}

We now wish to approximate the two-dimensional transport equation ($d=2$):
\begin{equation}
\label{transport2d}
\partial_t u +a_1 \partial_{x_1} u +a_2 \partial_{x_2} u =0,\quad u|_{t=0} =u_0,
\end{equation}
in the space domain $\{ x_1>0, x_2 \in \R \}$. When $a_1$ is negative, the latter problem does not necessitate
any boundary condition at $x_1=0$. Following \cite{ag1}, we use one of the following two-dimensional versions of
the leap-frog scheme, either
\begin{equation}
\label{lf2-1}
u_{j,k}^{n+2} +\lambda_1 a_1 (u_{j+1,k}^{n+1}-u_{j-1,k}^{n+1})
+\lambda_2 a_2 (u_{j,k+1}^{n+1}-u_{j,k-1}^{n+1}) -u_{j,k}^n =0,
\end{equation}
or
\begin{multline}\label{lf2-2}
u_{j,k}^{n+2}+\lambda_1 a_1 \Bigl(\dfrac{u_{j+1,k+1}^{n+1}+u_{j+1,k-1}^{n+1}}{2}
-\dfrac{u_{j-1,k+1}^{n+1}+u_{j-1,k-1}^{n+1}}{2} \Bigr) \\
+\lambda_2 a_2 \Bigl(\dfrac{u_{j+1,k+1}^{n+1}+u_{j-1,k+1}^{n+1}}{2}
-\dfrac{u_{j+1,k-1}^{n+1}+u_{j-1,k-1}^{n+1}}{2} \Bigr) -u_{j,k}^n =0.
\end{multline}
Assumption \ref{assumption0} is trivially satisfied because \eqref{lf2-1} and \eqref{lf2-2} are explicit schemes.
The scheme \eqref{lf2-1} satisfies Assumption \ref{assumption1} if and only if $\lambda_1 |a_1| +\lambda_2
|a_2|<1$, while the scheme \eqref{lf2-2} satisfies Assumption \ref{assumption1} if and only if $\max (\lambda_1
|a_1|,\lambda_2 |a_2|)<1$. Let us now study when Assumption \ref{assumption2} is valid. For the scheme
\eqref{lf2-1}, we have
$r_1=p_1=1$, and
\begin{equation*}
a_1(z,\eta) =\lambda_1 a_1 z,\quad a_{-1}(z,\eta) =-a_1(z,\eta),
\end{equation*}
so Assumption \ref{assumption2} is valid as long as $a_1 \neq 0$. For the scheme \eqref{lf2-2}, we have
again $r_1=p_1=1$, and
\begin{equation*}
a_1(z,\eta) =z (\lambda_1 a_1 \cos \eta +i \lambda_2 a_2 \sin \eta),\quad
a_{-1}(z,\eta) =z (-\lambda_1 a_1 \cos \eta +i \lambda_2 a_2 \sin \eta),
\end{equation*}
so Assumption \ref{assumption2} is valid as long as both $a_1$ and $a_2$ are nonzero. We refer to \cite{ag2}
for the verification of strong stability depending on the choice of some numerical boundary conditions for
\eqref{lf2-1} or \eqref{lf2-2}. If strong stability holds, then Theorem \ref{mainthm} yields semigroup boundedness
and $\ell^2$-stability with respect to boundary data.

\subsubsection*{Acknowledgments}
This article was completed while the author was visiting the Institute
of Mathematical Science at Nanjing University. The author warmly thanks Professor Yin Huicheng and the Institute for their hospitality during this visit.

\section{The Leray-G{\aa}rding method for fully discrete Cauchy problems}
\label{section2}

This section is devoted to proving stability estimates for discretized Cauchy problems, which is the first step
before considering the discretized initial boundary value problem \eqref{numibvp}. More precisely, we consider
the simpler case of the whole space $j \in \Z^d$, and the recurrence relation:
\begin{equation}
\label{numcauchy}
\begin{cases}
{\dps \sum_{\sigma=0}^{s+1}} Q_\sigma u_j^{n+\sigma} =0,&
j\in \Z^d,\quad n\ge 0,\\
u_j^n = f_j^n,& j\in \Z^d,\quad n=0,\dots,s,
\end{cases}
\end{equation}
where the operators $Q_\sigma$ are given by \eqref{defop}. We recall that in \eqref{defop}, the $a_{\ell,\sigma}$
are real numbers and are independent of the small parameter $\Delta t$ (they may depend on the CFL parameters
$\lambda_1,\dots,\lambda_d$), while $\bS$ denotes the shift operator on the space grid: $(\bS^\ell v)_j :=
v_{j+\ell}$ for $j,\ell \in \Z^d$. Stability of \eqref{numcauchy} is defined as follows.

\begin{definition}[Stability for the discrete Cauchy problem]
\label{def2}
The numerical scheme defined by \eqref{numcauchy} is ($\ell^2$-) stable if $Q_{s+1}$ is an isomorphism from
$\ell^2 (\Z^d)$ onto itself, and if furthermore there exists a constant $C_0>0$ such that for all $\Delta t \in (0,1]$,
for all initial conditions $(f_j^0)_{j \in \Z^d},\dots,(f_j^s)_{j \in \Z^d}$ in $\ell^2 (\Z^d)$, there holds
\begin{equation}
\label{estimcauchy}
\sup_{n \in \N} \Nd u^n \Nd_{-\infty,+\infty}^2 \le C_0 \sum_{\sigma=0}^s \Ng f^\sigma \Nd_{-\infty,+\infty}^2.
\end{equation}
\end{definition}

\noindent Let us quickly recall, see e.g. \cite{gko}, that stability in the sense of Definition \ref{def2} is in fact
independent of $\Delta t \in (0,1]$ (because \eqref{numcauchy} does not involve $\Delta t$ and \eqref{estimcauchy}
can be simplified on either side by $\prod_i \Delta x_i$), and can be characterized in terms of the uniform power
boundedness of the so-called amplification matrix
\begin{equation}
\label{defA2pas}
{\mathcal A}(\kappa) :=
\begin{pmatrix}
-\dfrac{\widehat{Q_s}(\kappa)}{\widehat{Q_{s+1}}(\kappa)} & & \hspace*{3mm}& \dots &\hspace*{5mm}& -\dfrac{\widehat{Q_0}(\kappa)}{\widehat{Q_{s+1}}(\kappa)} \\[10pt]
\hspace*{5mm}1& 0 && \dots &&\hspace*{-5mm}0 \\
& \ddots&& \ddots& & \\
\hspace*{5mm}0& && \ddots& \hspace*{2mm}\ddots&\hspace*{-5mm}\vdots \\
\hspace*{5mm}0& 0 && &1&\hspace*{-5mm}0 \end{pmatrix}\in {\mathcal M}_{s+1}(\C),
\end{equation}
where the $\widehat{Q_\sigma}(\kappa)$'s are defined in \eqref{dispersion} and where it is understood that ${\mathcal A}$
is defined on the largest open set of $\C^d$ on which $\widehat{Q_{s+1}}$ does not vanish. Let us also recall that if~$Q_{s+1}$ is an isomorphism from $\ell^2(\Z^d)$ onto itself, then $\widehat{Q_{s+1}}$ does not vanish on~$(\cercle)^d$,
and therefore does not vanish on an open neighborhood of~$(\cercle)^d$. With the above definition \eqref{defA2pas} for
${\mathcal A}$, the following well-known result holds, see e.g. \cite{gko}:

\begin{proposition}[Characterization of stability for the fully discrete Cauchy problem]
\label{prop1}
Assume that $Q_{s+1}$ is an isomorphism from $\ell^2(\Z^d)$ onto itself. Then the scheme \eqref{numcauchy} is
stable in the sense of Definition \ref{def2} if and only if there exists a constant $C_1>0$ such that the amplification
matrix ${\mathcal A}$ in \eqref{defA2pas} satisfies
\begin{equation*}
\forall n \in \N,\quad \forall \xi \in \R^d,\quad
\left\| {\mathcal A}(\rme^{i \xi_1},\dots,\rme^{i \xi_d})^n \right\| \le C_1.
\end{equation*}
In particular, the spectral radius of ${\mathcal A}(\rme^{i \xi_1},\dots,\rme^{i \xi_d})$ should not be
larger than $1$ (the so-called von Neumann condition).
\end{proposition}

The eigenvalues of ${\mathcal A}(\rme^{i \xi_1},\dots,\rme^{i \xi_d})$ are the roots to the dispersion
relation \eqref{dispersion}. When these roots are simple for all $\xi \in \R^d$, the von Neumann condition is
both necessary and {\it sufficient} for stability of \eqref{numcauchy}, see, e.g., \cite[Prop.\,3]{jfcnotes}.
Assumption~\ref{assumption1} is therefore a way to assume that \eqref{numcauchy} is stable for the discrete
Cauchy problem. Let us also recall that for each eigenvalue of ${\mathcal A}$, the corresponding eigenspace
has dimension~$1$ since ${\mathcal A}$ is a companion matrix. Therefore, assuming that the roots to the
dispersion relation \eqref{dispersion} are simple is equivalent to assuming that ${\mathcal A}$ is diagonalizable.

Our goal is to derive the semigroup estimate \eqref{estimcauchy} not by applying Fourier
transform to \eqref{numcauchy} and using uniform power boundedness of ${\mathcal A}$, but rather by
multiplying the first equation in \eqref{numcauchy} by a suitable {\it local} multiplier. The analysis relies
first on the simpler case where one only considers the time evolution and no additional space variable.

\subsection{Stable recurrence relations}

In this section, we consider sequences $(v^n)_{n \in \N}$ with values in $\C$. The index $n$ should be thought
of as the discrete time variable, and we therefore introduce the new notation $\bT$ for the shift operator on the
time grid: $(\bT^m v)^n :=v^{n+m}$ for all $m,n\in \N$. We start with the following elementary but crucial Lemma,
which is the analogue of \cite[Lem.\,1.1]{garding}.

\begin{lemma}[The energy-dissipation balance law]
\label{lem1}
Let $P \in \C[X]$ be a polynomial of degree $s+1$ whose roots are simple and located in $\Dbar$. Then
there exists a positive definite Hermitian form $q_e$ on $\C^{s+1}$, and a nonnegative Hermitian form $q_d$
on $\C^{s+1}$, that both depend in a ${\mathcal C}^\infty$ way on $P$, such that for any sequence $(v^n)_{n \in \N}$
with values in $\C$, there holds
\begin{multline*}
\forall n \in \N, \quad 2 \reel \Big(\overline{\bT (P'(\bT) v^n)} P(\bT) v^n \Big)\\
=(s+1) |P(\bT) v^n|^2 +(\bT-I) (q_e(v^n,\dots,v^{n+s})) +q_d(v^n,\dots,v^{n+s}).
\end{multline*}
In particular, for all sequence $(v^n)_{n \in \N}$ that satisfies the recurrence relation
\begin{equation*}
\forall n \in \N,\quad P(\bT) v^n =0,
\end{equation*}
the sequence $(q_e(v^n,\dots,v^{n+s}))_{n\in \N}$ is non-increasing.
\end{lemma}

The fact that there exists a Hermitian norm on $\C^{s+1}$ that is non-increasing along solutions to the
recurrence relation is not new. In fact, it is easily seen to be a consequence of the Kreiss matrix Theorem,
see \cite{strikwerda-wade}. However, the important point here is that we can construct a multiplier that yields
directly the `energy boundedness' (or~decay). The fact that the coefficients of this multiplier are integer
multiples of the coefficients of $P$ will be crucial in the analysis of Section \ref{section3}, see also Proposition
\ref{prop2} below.

\begin{proof}
We borrow some ideas from \cite[Lem.\,1.1]{garding} and introduce the interpolation polynomials:
\begin{equation*}
\forall k=1,\dots,s+1,\quad P_k(X) :=a \prod_{j \neq k} (X-x_j),
\end{equation*}
where $x_1,\dots,x_{s+1}$ denote the roots of $P$, and $a \neq 0$ its dominant coefficient. Since the roots of
$P$ are pairwise distinct, the $P_k$'s form a basis of $\C_s[X]$ and they depend in a ${\mathcal C}^\infty$ way
on the coefficients of $P$. We have
\begin{equation*}
P'=\sum_{k=1}^{s+1} P_k.
\end{equation*}
We then consider a sequence $(v^n)_{n \in \N}$ with values in $\C$ and compute\footnote{Here we use repeatedly
the property $\bT (\overline{w^n})=\overline{w^{n+1}} =\overline{\bT w^n}$, as well as $\overline{\bT w^n}
\bT w^n =|w^{n+1}|^2 =\bT |w^n|^2$.}
\begin{align*}
2 \reel \Big(&\overline{\bT (P'(\bT) v^n)} P(\bT) v^n \Big) - (s+1) |P(\bT) v^n|^2 \\[-5pt]
&= \sum_{k=1}^{s+1} \overline{\bT (P_k(\bT)) v^n} (\bT -x_k) P_k(\bT) v^n
+\bT (P_k(\bT) v^n) (\bT -\overline{x_k}) \overline{P_k(\bT) v^n} \\[-5pt]
&\hspace*{2cm}{}-\sum_{k=1}^{s+1} (\bT -\overline{x_k}) (\overline{P_k(\bT) v^n})
(\bT -x_k) (P_k(\bT) v^n) \\[-5pt]
&= \sum_{k=1}^{s+1} (\bT -|x_k|^2) |P_k(\bT) v^n|^2.
\end{align*}
The conclusion follows by defining:
\begin{align}
q_e(w^0,\dots,w^s) &:=\sum_{k=1}^{s+1} |P_k(\bT) w^0|^2,
\label{defqe}\\[-9pt]
\tag*{\hspace*{1cm}$\forall (w^0,\dots,w^s) \in \C^{s+1}$,\hspace*{1.5cm}}
\\[-9pt]
q_d(w^0,\dots,w^s) &:=\sum_{k=1}^{s+1} (1-|x_k|^2) |P_k(\bT) w^0|^2.\label{defqd}
\end{align}
The form $q_e$ is positive definite because the $P_k$'s form a basis of $\C_s[X]$. The form $q_d$ is
nonnegative because the roots of $P$ are located in $\Dbar$. Both forms depend in a ${\mathcal C}^\infty$
way on the coefficients of $P$ because the roots of $P$ are simple.
\end{proof}

Lemma \ref{lem1} shows that the polynomial $P'$ yields the good multiplier $\bT P'(\bT) v^n$ for
the recurrence relation $P(\bT) v^n =0$. Of course, $P'$ is not the only possible choice, though it will
be our favorite one in what follows. As in \cite[Lem.\,1.1]{garding}, any polynomial of the form\footnote{The
sign condition here on the coefficients $\alpha_k$ is the analogue of the separation condition for the roots in
\cite{leray,garding}.}
\begin{equation*}
Q :=\sum_{k=1}^{s+1} \alpha_k P_k,\quad \alpha_1,\dots,\alpha_{s+1} >0,
\end{equation*}
provides with an energy balance of the form
\begin{multline*}
2 \reel \Big(\overline{\bT (Q(\bT) v^n)} P(\bT) v^n \Big)\\
=(\alpha_1+\cdots+\alpha_{s+1}) |P(\bT) v^n|^2 +(\bT-I) (q_e(v^n,\dots,v^{n+s}))
+q_d(v^n,\dots,v^{n+s}),
\end{multline*}
with suitable Hermitian forms $q_e,q_d$ that have the same properties as stated in Lemma~\ref{lem1}.

\Subsection{The energy-dissipation balance for finite difference schemes}

In this section, we consider the numerical scheme \eqref{numcauchy}. We introduce the following notation:
\begin{equation}
\label{defLM}
L:= \sum_{\sigma=0}^{s+1} \bT^\sigma Q_\sigma,\quad
M:= \sum_{\sigma=0}^{s+1} \sigma \bT^\sigma Q_\sigma,
\end{equation}
so that the discretized Cauchy problem \eqref{numcauchy} reads
\begin{equation*}
\begin{cases}
L u_j^n =0,& j\in \Z^d,\quad n\ge 0,\\
u_j^n = f_j^n,& j\in \Z^d,\quad n=0,\dots,s.
\end{cases}
\end{equation*}
The operator $M$ will be the `multiplier' associated with $L$. Thanks to Fourier analysis, Lemma \ref{lem1} easily
gives the following result:

\begin{proposition}[The energy-dissipation balance law]
\label{prop2}
Let Assumptions \ref{assumption0} and \ref{assumption1} be satisfied. Then there exist a continuous coercive
quadratic form $E_0$ and a continuous nonnegative quadratic form $D_0$ on $\ell^2(\Z^d;\R)^{s+1}$ such that
for all sequences $(v^n)_{n \in \N}$ with values in $\ell^2(\Z^d;\R)$ and for all $n \in \N$, there holds
\begin{multline*}
2 \langle M v^n,L v^n \rangle_{-\infty,+\infty}\\
=(s+1) \Ng L v^n \Nd_{-\infty,+\infty}^2
+(\bT-I) E_0(v^n,\dots,v^{n+s}) +D_0(v^n,\dots,v^{n+s}).
\end{multline*}
In particular, for all initial data $f^0,\dots,f^s \in \ell^2(\Z^d;\R)$, the solution to \eqref{numcauchy} satisfies
\begin{equation*}
\sup_{n \in \N} E_0(v^n,\dots,v^{n+s}) \le E_0(f^0,\dots,f^s),
\end{equation*}
and \eqref{numcauchy} is ($\ell^2$-)stable.
\end{proposition}

\begin{proof}
We use the same notation $v^n$ for the sequence $(v_j^n)_{j \in \Z^d}$ and the corresponding step function on
$\R^d$ whose value on the cell $[j_1 \Delta x_1,(j_1+1) \Delta x_1) \times \cdots \times [j_d \Delta x_d,(j_d+1)
\Delta x_d)$ equals $v_j^n$. Then Plancherel Theorem gives
\begin{multline*}
2 \langle M v^n,L v^n \rangle_{-\infty,+\infty}-(s+1) \Ng L v^n \Nd_{-\infty,+\infty}^2 \\
=\int_{\R^d} 2 \reel \Big(\overline{\bT (P_\zeta'(\bT) \widehat{v^n}(\xi))} P_\zeta(\bT)
\widehat{v^n} (\xi) \Big) -(s+1) \bigl|P_\zeta(\bT) \widehat{v^n}(\xi)\bigr|^2 \dfrac{\rmd\xi}{(2 \pi)^d},
\end{multline*}
where $\widehat{v^n}$ denotes the Fourier transform of $v^n$, and where we have let
\begin{equation*}
P_\zeta (z):=\sum_{\sigma=0}^{s+1} \widehat{Q_\sigma} \big(\rme^{i \zeta_1},\dots,\rme^{i \zeta_d} \big)
z^\sigma,\quad \zeta_j := \xi_j \Delta x_j,
\end{equation*}
and $P'_\zeta(z)$ denotes the derivative of $P_\zeta$ with respect to $z$.

From Assumption \ref{assumption1}, we know that for all $\zeta \in \R^d$, $P_\zeta$ has degree $s+1$ and has
$s+1$ simple roots in $\Dbar$. We can apply Lemma \ref{lem1} and get
\begin{multline*}
2 \langle M v^n,L v^n \rangle_{-\infty,+\infty}-(s+1) \Ng L v^n \Nd_{-\infty,+\infty}^2 \\
=\int_{\R^d} (\bT-I) q_{e,\zeta} \big(\widehat{v^n}(\xi),\dots,\widehat{v^{n+s}}(\xi) \big)
+q_{d,\zeta} \big(\widehat{v^n}(\xi),\dots,\widehat{v^{n+s}}(\xi) \big) \,\dfrac{\rmd\xi}{(2 \pi)^d},
\end{multline*}
where $q_{e,\zeta},q_{d,\zeta}$ depend in a ${\mathcal C}^\infty$ way on $\zeta \in \R^d$ and are $2 \pi$-periodic
in each $\zeta_j$. Furthermore, $q_{e,\zeta}$ is positive definite and $q_{d,\zeta}$ is nonnegative. We then define
\begin{align*}
\tag*{\hspace*{1mm}}
E_0(w^0,\dots,w^s) &:=\int_{\R^d} q_{e,\zeta} \big(\widehat{w^0}(\xi),\dots,\widehat{w^s}(\xi) \big)
\dfrac{\rmd\xi}{(2 \pi)^d},\\[-5pt]
\tag*{$\forall (w^0,\dots,w^s) \in \ell^2(\Z^d;\R)^{s+1}$,\hspace*{2.2cm}}
\\[-5pt]
\tag*{\hspace*{1mm}}
D_0(w^0,\dots,w^s) &:=\int_{\R^d} q_{d,\zeta} \big(\widehat{w^0}(\xi),\dots,\widehat{w^s}(\xi) \big)
\dfrac{\rmd\xi}{(2 \pi)^d}.
\end{align*}
By compactness of $[0,2 \pi]^d$, the hermitian forms $q_{e,\zeta},q_{d,\zeta}$ are uniformly bounded with respect
to $\zeta$ (because they depend continuously on $\zeta$). Therefore $E_0$ and $D_0$ define continuous quadratic
forms on $\ell^2(\Z^d;\R)^{s+1}$, and $D_0$ is nonnegative as the integral of a nonnegative quantity. Furthermore,
$q_{e,\zeta}$ is a positive Hermitian form for each $\zeta$ so by compactness it is uniformly coercive with respect to
$\zeta$. Hence $E_0$ is a coercive quadratic form and the proof of Proposition \ref{prop2} is complete.
\end{proof}

\subsection{Examples}

Let us clarify Proposition \ref{prop2} in the case of the one-dimensional examples of Section \ref{examples}.
For the leap-frog scheme in one space dimension, there holds
\begin{equation*}
L=\bT^2 +\lambda a \bT (\bS-\bS^{-1}) -I,
\end{equation*}
and our multiplier $M u_j^n$ reads
\begin{equation*}
M u_j^n =2 u_j^{n+2} +\lambda a (u_{j+1}^{n+1}-u_{j-1}^{n+1})
=u_j^{n+2}+u_j^n +\underbrace{L u_j^n}_{=0},
\end{equation*}
where the equality $L u_j^n =0$ is used as long as $(u_j^n)$ corresponds to a solution to the leap-frog scheme.
We thus recover the more classical multiplier $u_j^{n+2}+u_j^n$ used in~\cite{RM}, but we emphasize that both
multipliers coincide only on solutions to the leap-frog scheme. It will appear more clearly in Section \ref{section3}
why our choice for $M u_j^n$ has a major advantage when considering initial boundary value problems because
the main energy-dissipation identity of Proposition \ref{prop2} not only holds for solutions to $L u_j^n =0$ but for
any sequence $(u^n)$ with values in $\ell^2 (\Z^d)$. Let us now look at the energy and dissipation
functionals provided by Proposition \ref{prop2}. Since the (two simple) roots $z_1,z_2$ to \eqref{dispersionlf}
have modulus $1$, if we keep the notation of the proof of Lemma~\ref{lem1} and Proposition \ref{prop2}, we get
the expressions
\begin{align*}
q_{e,\zeta} (w^0,w^1) &= \big| w^1-z_1(\zeta) w^0 \big|^2 +\big| w^1-z_2(\zeta) w^0 \big|^2\\
&=2 (|w^0|^2 +|w^1|^2) +4 \reel \big(i \lambda a \sin \zeta \overline{w^1} w^0 \big),\\
q_{d,\zeta} (w^0,w^1) &= 0,
\end{align*}
for all $\zeta \in [0,2 \pi]$. After substituting $\zeta =\xi \Delta x$ and integrating with respect to $\xi$, we~get
\begin{equation*}
E_0(v^n,v^{n+1}) =2 \sum_{j \in \Z} \Delta x (v_j^n)^2 +\Delta x (v_j^{n+1})^2
+2 \sum_{j \in \Z} \Delta x \lambda a \Big(v_{j+1}^n-v_{j-1}^n \Big) v_j^{n+1},
\end{equation*}
and $D_0 \equiv 0$ (no dissipation). Proposition \ref{prop2} shows that the energy functional $E_0$ is preserved
for solutions to the leap-frog scheme, a fact that also comes more directly from the relation
\begin{equation*}
2 \sum_{j \in \Z} \Delta x (u_j^{n+2} +u_j^n) \Big(
u_j^{n+2} +\lambda a (u_{j+1}^{n+1}-u_{j-1}^{n+1}) -u_j^n \Big) =0,
\end{equation*}
after regrouping
\begin{multline*}
2 \sum_{j \in \Z} \Delta x (u_j^{n+2} +u_j^n) (u_j^{n+2} -u_j^n)\\[-5pt]
=
2 \sum_{j \in \Z} \Delta x (u_j^{n+1})^2 +\Delta x (u_j^{n+2})^2
-2 \sum_{j \in \Z} \Delta x (u_j^n)^2 +\Delta x (u_j^{n+1})^2,
\end{multline*}
and using summation by parts for the remaining terms. What is important here is that both quantities $|z_1(\zeta)|^2
+|z_2(\zeta)|^2$ and $z_1(\zeta)+z_2(\zeta)$ are trigonometric polynomials in~$\zeta$, and therefore the expression
of $E_0$ turns out to give a finite sum of `local' quadratic functionals of one of the forms
\begin{equation*}
\sum_{j \in \Z} \Delta x v_{j+\ell_1}^n v_{j+\ell_2}^n,\quad
\sum_{j \in \Z} \Delta x v_{j+\ell_1}^n v_{j+\ell_2}^{n+1},\quad
\sum_{j \in \Z} \Delta x v_{j+\ell_1}^{n+1} v_{j+\ell_2}^{n+1}.
\end{equation*}
This means that we could also define the energy functional $E_0$ for sequences $(v_j)$ defined only for $j \ge 0$
and not on all $\Z$ by `localizing' in space (this fact is used in~\cite{rauch} in the context of partial differential equations
because the energy functional that arises in \cite{leray,garding} turns out to be a local quantity).

Let us now turn to the scheme \eqref{bdf2} that is based on the backwards differentiation formula of order $2$.
In that case, we have
\begin{equation*}
L=\bT^2 \Bigl(\dfrac{3}{2} I +\dfrac{\lambda a}{2} (\bS -\bS^{-1}) \Bigr)
-2 \bT +\dfrac{1}{2} I,\quad
M=\bT^2 \left(3 I +\lambda a (\bS -\bS^{-1}) \right) -2 \bT,
\end{equation*}
so even if $L u_j^n=0$, our multiplier $M u_j^n$ does not coincide with the `more standard'~$u_j^{n+2}$ used
in \cite{Emmrich1,Emmrich2}. Here the multiplier $M u_j^n$ incorporates some terms that take into account the
spatial discretization while the multiplier $u_j^{n+2}$ is designed to take advantage of the $G$-stability of the
BDF-$2$ integration rule, see \cite[Chap.\,V.6]{hw}. The energy-dissipation functionals provided by Proposition
\ref{prop2} are not as elegant in this case as what they were in the case of the leap-frog scheme. Namely, we
keep the notation of the proofs of Lemma \ref{lem1} and Proposition \ref{prop2}. Letting $z_1,z_2$ denote the
roots to the dispersion relation \eqref{dispersionbdf}, and introducing the notation
\begin{equation*}
{\mathcal Q}(\zeta) := \Bigl | \dfrac{3}{2} +i \lambda a \sin \zeta \Bigr|^2 \bigl(|z_1(\zeta)|^2 +|z_2(\zeta)|^2
\bigr),
\end{equation*}
we compute
\begin{align*}
q_{e,\zeta} (w^0,w^1) &= \Bigl | \dfrac{3}{2} +i \lambda a \sin \zeta \Bigr |^2 \big(
\big| w^1-z_1(\zeta) w^0 \big|^2 +\big| w^1-z_2(\zeta) w^0 \big|^2 \big) \\
&=2 \Bigl | \Bigl(\dfrac{3}{2} +i \lambda a \sin \zeta \Bigr) w^1 \Bigr|^2
+{\mathcal Q}(\zeta) |w^0|^2 -4 \reel \Bigl(\overline{\Bigl(\dfrac{3}{2} +i \lambda a \sin \zeta \Bigr)
w^1} w^0 \Bigr),
\end{align*}
\begin{multline*}
q_{d,\zeta} (w^0,w^1) = \Bigl(2 \Bigl| \dfrac{3}{2} +i \lambda a \sin \zeta \Bigr|^2 -{\mathcal Q}(\zeta) \Bigr)
|w^1|^2 +\bigl({\mathcal Q}(\zeta) -\dfrac{1}{2} \bigr) |w^0|^2 \\
-4 \reel \Bigl(\overline{\Bigl(\dfrac{3}{2} +i \lambda a \sin \zeta \Bigr) w^1} w^0 \Bigr)
-2 \Bigl | \dfrac{3}{2} +i \lambda a \sin \zeta \Bigr|^2 \reel \bigl(\big(|z_1|^2 z_2 +|z_2|^2 z_1 \big)
\overline{w^1} w^0 \bigr),
\end{multline*}
for all $\zeta \in [0,2 \pi]$. The energy functional is then defined as
\begin{equation*}
E_0(v^n,v^{n+1})
=\int_\R q_{e,\xi \Delta x} \big(\widehat{v^n}(\xi),\widehat{v^{n+1}}(\xi) \big) \,\dfrac{\rmd\xi}{2 \pi},
\end{equation*}
and using the above expression for the Hermitian form $q_{e,\zeta}$, we get
\begin{multline*}
E_0(v^n,v^{n+1}) =2 \sum_{j \in \Z} \Delta x \Bigl(\dfrac{3}{2} v_j^{n+1} +\dfrac{\lambda a}{2}
(v_{j+1}^{n+1}-v_{j-1}^{n+1}) \Bigl)^2\\
-4 \sum_{j \in \Z} \Delta x \Bigl(\dfrac{3}{2} v_j^{n+1}
+\dfrac{\lambda a}{2} (v_{j+1}^{n+1}-v_{j-1}^{n+1}) \Bigr) v_j^n
+\int_\R {\mathcal Q}(\xi \Delta x) |\widehat{v^n}(\xi)|^2 \,\dfrac{\rmd\xi}{2 \pi},
\end{multline*}
and rather similar expression for $D_0$. The problem at this stage is that there is no obvious reason for ${\mathcal Q}$
to be a trigonometric polynomial in $\zeta$ and therefore the last term in the decomposition of $E_0$ does not obviously
decompose as a linear combination of local energy functionals
\begin{equation*}
\sum_{j \in \Z} \Delta x v_{j+\ell_1}^n v_{j+\ell_2}^n.
\end{equation*}
Though the energy functional $E_0$ will be sufficient for our purpose here, its nonlocal feature may prevent from
extending the current multiplier technique to obtain stability results on non-Cartesian meshes. We leave this question
to further study.

\section{Semigroup estimates for fully discrete initial boundary value problems}
\label{section3}

We now turn to the proof of Theorem \ref{mainthm} for which we shall use the results of Section~\ref{section2}
as a toolbox. By linearity of \eqref{numibvp}, it is sufficient to prove Theorem \ref{mainthm} separately in the
case $(f_j^0)=\cdots=(f_j^s)=0$, and in the case $(F_j^n)=0$, $(g_j^n)=0$. The latter case is the most difficult
and requires the introduction of an auxiliary set of `dissipative' boundary conditions. Solutions to \eqref{numibvp}
are always assumed to be real valued, which means that the data are real valued. For complex valued initial data
and/or forcing terms, one just uses the linearity of \eqref{numibvp}.

\subsection{The case with zero initial data}

We first assume $(f_j^0)=\cdots=(f_j^s)=0$. By~strong stability, we already know that \eqref{stabilitenumibvp}
holds with a constant $C$ that is independent of $\gamma>0$ and $\Delta t \in (0,1]$. Therefore, proving
Theorem \ref{mainthm} amounts to showing the existence of a constant $C$, that is independent of $\gamma>0$
and $\Delta t \in (0,1]$ such that the solution to \eqref{numibvp} with $(f_j^0)=\cdots=(f_j^s)=0$ satisfies
\begin{multline}
\label{estimnumibvp'}
\sup_{n \ge 0} \rme^{-2 \gamma n \Delta t} \Ng u^n \Nd_{1-r_1,+\infty}^2
\le C \biggl\{ \dfrac{\gamma \Delta t+1}{\gamma}
\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \Ng F^n \Nd_{1,+\infty}^2\\[-5pt]
+\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1=1-r_1}^0
\| g_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \biggr\}.
\end{multline}
We thus consider a parameter $\gamma>0$ and a time step $\Delta t \in (0,1]$, and focus on the numerical scheme
\eqref{numibvp} with zero initial data (that is, $(f_j^0)=\cdots=(f_j^s)=0$). For all $n \in \N$, we extend the sequence
$(u_j^n)$ by zero for $j_1 \le -r_1$:
\begin{equation*}
v_j^n :=\begin{cases}
u_j^n &\text{\rm if } j_1 \ge 1-r_1,\quad j' \in \Z^{d-1},\\
0 &\text{\rm otherwise.}
\end{cases}
\end{equation*}
Observe that $L v_j^n$ is not zero for all $j \in \Z^d$ hence the need for the general framework of Proposition
\ref{prop2}. We thus use Proposition \ref{prop2} and compute:
\begin{multline*}
(\bT-I) E_0(v^n,\dots,v^{n+s}) +D_0(v^n,\dots,v^{n+s})\\
=
2 \langle M v^n,L v^n \rangle_{-\infty,+\infty} -(s+1) \Ng L v^n \Nd_{-\infty,+\infty}^2.
\end{multline*}
Due to the form of the operator $L$, see \eqref{defLM}, and the fact that $v_j^n$ vanishes for $j_1 \le -r_1$, there
holds:
\begin{equation*}
L v_j^n =\begin{cases}
\Delta t F_j^{n+s+1} &\text{\rm if } j_1 \ge 1,\\
0 &\text{\rm if } j_1 \le -r_1-p_1,
\end{cases}
\end{equation*}
and we thus get
\begin{multline*}
(\bT-I) E_0(v^n,\dots,v^{n+s}) + D_0(v^n,\dots,v^{n+s}) \\
= \biggl(\prod_{k=1}^d \Delta x_k \biggr) \sum_{j_1 \ge 1} \sum_{j' \in \Z^{d-1}}
2 \Delta t (M v_j^n) F_j^{n+s+1} -(s+1) \Delta t^2 (F_j^{n+s+1})^2 \\
+ \biggl(\prod_{k=1}^d \Delta x_k \biggr) \sum_{j_1=1-r_1-p_1}^0 \sum_{j' \in \Z^{d-1}}
2 (M v_j^n) L v_j^n -(s+1) (L v_j^n)^2.
\end{multline*}
We multiply the latter equality by $\exp(-2 \gamma (n+s+1) \Delta t)$, sum with respect to $n$
from~$0$ to some $N$ and use the fact that $D_0$ is nonnegative. Recalling that the initial data in
\eqref{numibvp} vanish, we get
\begin{multline}
\label{energy1}
\rme^{-2 \gamma (N+s+1) \Delta t} E_0 \big(v^{N+1},\dots,v^{N+s+1} \big)\\[-5pt]
+\underbrace{\big(1-\rme^{-2 \gamma \Delta t} \big) \sum_{n=1}^N \rme^{-2 \gamma (n+s) \Delta t}
E_0 (v^n,\dots,v^{n+s})}_{\ge 0}\le S_{1,N} +S_{2,N},
\end{multline}
\vspace*{-10pt}%
\par\noindent
with
\begin{multline}
\label{defS1Ngamma}
S_{1,N}\\[-5pt]
:=\sum_{n=0}^N \rme^{-2 \gamma (n+s+1) \Delta t} \Big(2 \Delta t
\langle M v^n,F^{n+s+1} \rangle_{1,+\infty} -(s+1) \Delta t^2 \Ng F^{n+s+1} \Nd_{1,+\infty}^2 \Big),
\end{multline}
and
\begin{multline}
\label{defS2Ngamma}
S_{2,N}\\[-5pt]
:=\biggl(\prod_{k=1}^d \Delta x_k \biggr) \sum_{n=0}^N \rme^{-2 \gamma (n+s+1) \Delta t}
\sum_{j_1=1-r_1-p_1}^0 \sum_{j' \in \Z^{d-1}} 2 (M v_j^n) L v_j^n -(s+1) (L v_j^n)^2.
\end{multline}

Let us now estimate the two source terms $S_{1,N},S_{2,N}$ in \eqref{energy1}. We begin with the term $S_{2,N}$
defined in \eqref{defS2Ngamma}. Let us recall that the ratio $\Delta t/\Delta x_1$ is fixed\footnote{This is one first
occurrence where restricting to the `hyperbolic scaling' for the time and space steps is convenient.}. Furthermore,
the form of the operators $L$ and $M$ in \eqref{defLM} gives the estimate (recall that $v_j^n$ vanishes for
$j_1 \le -r_1$):
\begin{equation*}
S_{2,N} \le C \Delta t \biggl(\prod_{k=2}^d \Delta x_k \biggr) \sum_{n=0}^N
\rme^{-2 \gamma (n+s+1) \Delta t} \sum_{j_1=1-r_1}^{p_1} \sum_{j' \in \Z^{d-1}}
(u_j^n)^2 +\cdots +(u_j^{n+s+1})^2,
\end{equation*}
for a constant $C$ that does not depend on $N$, $\gamma$ nor on $\Delta t$. We thus have, uniformly with respect
to $N \in \N$, $\gamma>0$ and $\Delta t \in (0,1]$:
\begin{align*}
S_{2,N} &\le C \sum_{n=s+1}^{N+s+1} \Delta t \rme^{-2 \gamma n \Delta t}
\sum_{j_1=1-r_1}^{p_1} \| u_{j_1,\sbullet}^n \|_{\ell^2 (\Z^{d-1})}^2\\
&\le C \sum_{n \ge s+1} \Delta t \rme^{-2 \gamma n \Delta t}
\sum_{j_1=1-r_1}^{p_1} \| u_{j_1,\sbullet}^n \|_{\ell^2 (\Z^{d-1})}^2 \\
&\le C \biggl\{ \dfrac{\gamma \Delta t+1}{\gamma}
\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \Ng F^n \Nd_{1,+\infty}^2\\[-10pt]
&\hspace*{5cm}
+\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1=1-r_1}^0
\| g_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \biggr\},
\end{align*}
where we have used the trace estimate \eqref{stabilitenumibvp} that follows from the strong stability assumption.

Let us now focus on the term $S_{1,N}$ in \eqref{energy1}, see the defining equation \eqref{defS1Ngamma}.
We use the Cauchy-Schwarz inequality and derive (using now the interior estimate in \eqref{stabilitenumibvp}
that follows from the strong stability assumption and the fact that the coefficients in the multiplier $M$ are
independent of $\Delta t$):
\begin{align*}
&S_{1,N} \le 2 \sum_{n=0}^N \Delta t \rme^{-2 \gamma (n+s+1) \Delta t}
\Ng M v^n \Nd_{1,+\infty} \Ng F^{n+s+1} \Nd_{1,+\infty} \\
&\le C \sum_{n=0}^N\! \Delta t \rme^{-2 \gamma (n+s+1) \Delta t} \Bigl(
\Ng v^{n+1} \Nd_{1-r_1,+\infty} +\cdots +\Ng v^{n+s+1} \Nd_{1-r_1,+\infty} \Bigr) \Ng F^{n+s+1} \Nd_{1,+\infty} \\
&\le C \dfrac{\gamma}{\gamma \Delta t+1}
\sum_{n=s+1}^{N+s+1}\hspace*{-2mm} \Delta t \rme^{-2 \gamma n \Delta t} \Ng u^n \Nd_{1-r_1,+\infty}^2
+C \dfrac{\gamma \Delta t+1}{\gamma} \sum_{n=s+1}^{N+s+1}\hspace*{-2mm} \Delta t
\rme^{-2 \gamma n \Delta t} \Ng F^n \Nd_{1,+\infty}^2 \\
&\le C \biggl\{ \dfrac{\gamma \Delta t+1}{\gamma}
\sum_{n\ge s+1}\hspace*{-2mm} \Delta t \rme^{-2 \gamma n \Delta t} \Ng F^n \Nd_{1,+\infty}^2
+\sum_{n\ge s+1}\hspace*{-2mm} \Delta t \rme^{-2 \gamma n \Delta t}\hspace*{-2mm} \sum_{j_1=1-r_1}^0
\| g_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \biggl\}.
\end{align*}
Ignoring the nonnegative term on the left-hand side of \eqref{energy1} and using the coercivity of $E_0$, we
have proved that there exists a constant $C>0$ that is uniform with respect to $N,\gamma,\Delta t$ such that:
\begin{multline*}
\rme^{-2 \gamma (N+s+1) \Delta t} \Ng v^{N+s+1} \Nd_{-\infty,+\infty}^2 \le C
\biggl\{ \dfrac{\gamma \Delta t+1}{\gamma}
\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \Ng F^n \Nd_{1,+\infty}^2 \\[-5pt]
+\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1=1-r_1}^0
\| g_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \biggr\},
\end{multline*}
which yields \eqref{estimnumibvp'} and therefore the validity of Theorem \ref{mainthm} in the case of zero
initial data.

\subsection{Construction of dissipative boundary conditions}

In this paragraph, we consider an auxiliary problem for which we shall be able to prove simultaneously
an optimal semigroup estimate and a trace estimate for the solution. The argument here is independent
of the original numerical scheme \eqref{numibvp}, but the auxiliary scheme introduced in Theorem
\ref{absorbing} below will be used later on to decompose the solution to \eqref{numibvp} into two
pieces, each of which being estimated by separate tools. We thus forget temporarily about \eqref{numibvp}
and state the following key result.

\begin{theorem}
\label{absorbing}
Let Assumptions \ref{assumption0}, \ref{assumption1} and \ref{assumption2} be satisfied. Then for all $P_1
\in \N$, there exists a constant $C_{P_1}>0$ such that, for all initial data $(f^0_j),\dots,(f^s_j) \in \ell^2 (\Z^d)$
and for all source term $(g_j^n)_{j_1 \le 0,n \ge s+1}$ that satisfies
\begin{equation*}
\forall \Gamma >0,\quad \sum_{n\ge s+1} \rme^{-2 \Gamma n} \sum_{j_1 \le 0}
\| g_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 < +\infty,
\end{equation*}
there exists a unique sequence $(u_j^n)_{j \in \Z^d,n\in \N}$ solution to
\begin{equation}
\label{numabsorbing}
\begin{cases}
L u_j^n =0,& j_1 \ge 1,\quad j' \in \Z^{d-1},\quad n\ge 0,\\
M u_j^n =g_j^{n+s+1},& j_1 \le 0,\quad j' \in \Z^{d-1},\quad n\ge 0,\\
u_j^n = f_j^n,& j \in \Z^d,\quad n=0,\dots,s.
\end{cases}
\end{equation}
Moreover for all $\gamma>0$ and $\Delta t \in (0,1]$, this solution satisfies
\begin{multline}
\label{estimabsorb}
\sup_{n \ge 0} \rme^{-2 \gamma n \Delta t} \Ng u^n \Nd_{-\infty,+\infty}^2
+\dfrac{\gamma}{\gamma \Delta t+1}
\sum_{n\ge 0} \Delta t \rme^{-2 \gamma n \Delta t} \Ng u^n \Nd_{-\infty,+\infty}^2\\[-8pt]
\shoveright{+\sum_{n\ge 0} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1=1-r_1}^{P_1}
\| u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2}\\
\le C_{P_1} \biggl\{ \sum_{\sigma=0}^s \Ng f^\sigma \Nd_{-\infty,+\infty}^2
+\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1 \le 0}
\| g_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \biggr\}.
\end{multline}
\end{theorem}

Theorem \ref{absorbing} justifies why we advocate the choice $M u_j^n =2 u_j^{n+2} +\lambda a
(u_{j+1}^{n+1} -u_{j-1}^{n+1})$ rather than the more standard $u_j^{n+2}+u_j^n$ as a multiplier for the
leap-frog scheme. Despite repeated efforts, we have not been able to prove the estimate \eqref{estimabsorb}
when using the numerical boundary condition $u_j^{n+2}+u_j^n$ on $j_1 \le 0$, in conjunction with the
leap-frog scheme on $j_1 \ge 1$.

\begin{proof}
Let us first quickly observe that the solution to \eqref{numabsorbing} is well-defined since, as long as we have
determined the solution up to a time index $n+s$, $n \ge 0$, then $u^{n+s+1}$ is sought as a solution to an
equation of the form
\begin{equation*}
Q_{s+1} u^{n+s+1} =F,
\end{equation*}
where $F$ belongs to $\ell^2(\Z^d)$ (this is due to the form of $L$ and $M$, see \eqref{defLM}). Hence $u^n$
is uniquely defined and belongs to $\ell^2(\Z^d)$ for all $n \in \N$.

The proof of Theorem \ref{absorbing} starts again with the application of Proposition \ref{prop2}. Using
the non-negativity of the dissipation form $D_0$, we get\footnote{Since $L u_j^n=0$ for $j_1 \ge 1$, one
could also write $\Ng L u^n \Nd_{-\infty,0}^2$ rather than $\Ng L u^n \Nd_{-\infty,+\infty}^2$ on the left
hand-side of the inequality.}
\begin{multline*}
(\bT-I) E_0(u^n,\dots,u^{n+s}) +(s+1) \Ng L u^n \Nd_{-\infty,+\infty}^2\\
\le 2 \langle M u^n,L u^n \rangle_{-\infty,+\infty} =2 \langle g^{n+s+1},L u^n \rangle_{-\infty,0}.
\end{multline*}
By the Young inequality
\begin{equation*}
2 \langle g^{n+s+1},L u^n \rangle_{-\infty,0} \le \dfrac{s+1}{2} \Ng L u^n \Nd_{-\infty,0}^2
+\dfrac{2}{s+1} \Ng g^{n+s+1} \Nd_{-\infty,0}^2,
\end{equation*}
we get
\begin{equation*}
(\bT-I) E_0(u^n,\dots,u^{n+s}) +\dfrac{s+1}{2} \Ng L u^n \Nd_{-\infty,+\infty}^2
\le \dfrac{2}{s+1} \Ng g^{n+s+1} \Nd_{-\infty,0}^2.
\end{equation*}
We multiply the latter inequality by $\exp (-2 \gamma (n+s+1) \Delta t)$, sum from $n=0$ to some arbitrary
$N$ and already derive the estimate (here we use again the fact that $\Delta t/\Delta x_1$ is a fixed positive
constant):
\begin{multline*}
\sup_{n \ge 1} \rme^{-2 \gamma (n+s) \Delta t} E_0(u^n,\dots,u^{n+s})
+\big(1-\rme^{-2 \gamma \Delta t} \big)
\sum_{n \ge 0} \rme^{-2 \gamma (n+s) \Delta t} E_0(u^n,\dots,u^{n+s}) \\
\shoveright{+\sum_{n\ge 0} \Delta t \rme^{-2 \gamma (n+s+1) \Delta t} \sum_{j_1 \in \Z}
\| L u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2}\\
\le C \biggl\{ \rme^{-2 \gamma s \Delta t} E_0(f^0,\dots,f^s)
+\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1 \le 0}
\| g_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \biggl\}.
\end{multline*}
Using the coercivity of $E_0$ and the inequality
\begin{equation*}
1-\rme^{-2 \gamma \Delta t} \ge \dfrac{\gamma \Delta t}{\gamma \Delta t +1},
\end{equation*}
we have therefore derived the estimate
\begin{multline}
\label{estim1}
\sup_{n \ge 0} \rme^{-2 \gamma n \Delta t} \Ng u^n \Nd_{-\infty,+\infty}^2
+\dfrac{\gamma}{\gamma \Delta t +1} \sum_{n\ge 0} \Delta t \rme^{-2 \gamma n \Delta t}
\Ng u^n \Nd_{-\infty,+\infty}^2 \\
\shoveright{+\sum_{n\ge 0} \Delta t \rme^{-2 \gamma (n+s+1) \Delta t} \sum_{j_1 \in \Z}
\| L u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2}\\
\le C \biggl\{ \sum_{\sigma=0}^s \Ng f^\sigma \Nd_{-\infty,+\infty}^2
+\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1 \le 0}
\| g_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \biggl\},
\end{multline}
where the constant $C$ is independent of $\gamma$, $\Delta t$ and on the solution $(u_j^n)$. In order to
prove \eqref{estimabsorb}, the main remaining task is to derive the trace estimate for $(u_j^n)$. This is done
by first dealing with the case where $\gamma \Delta t$ is `large'.

$\bullet$ From the definition of the operator $L$, see \eqref{defLM}, there exists a constant $C>0$ and an
integer $J$ such that
\begin{equation*}
(L u_j^n)^2 \ge \dfrac{1}{2} (Q_{s+1} u_j^{n+s+1})^2 -C \sum_{\sigma=0}^s \sum_{|\ell| \le J}
(u_{j+\ell}^{n+\sigma})^2.
\end{equation*}
Since $Q_{s+1}$ is an isomorphism, there exists a constant $c>0$ such that
\begin{equation*}
\sum_{j \in \Z^d} (L u_j^n)^2 \ge c \sum_{j \in \Z^d} (u_j^{n+s+1})^2 -\dfrac{1}{c} \sum_{\sigma=0}^s
\sum_{j \in \Z^d} (u_j^{n+\sigma})^2.
\end{equation*}
Multiplying by $\exp (-2 \gamma (n+s+1) \Delta t)$ and summing with respect to $n \in \N$, we get
\begin{multline}
\label{estimaux}
\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1 \in \Z}
\| u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2\\
\le C \biggl\{
\sum_{n\ge 0} \Delta t \rme^{-2 \gamma (n+s+1) \Delta t} \sum_{j_1 \in \Z}
\| L u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2\\
+\rme^{-2 \gamma \Delta t} \sum_{n\ge 0} \Delta t \rme^{-2 \gamma n \Delta t}
\sum_{j_1 \in \Z} \| u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \biggr\}.
\end{multline}
The second term on the right-hand side is decomposed as
\begin{multline*}
\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1 \in \Z} \| u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2
+\sum_{\sigma=0}^s
\Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1 \in \Z} \| f_{j_1,\sbullet}^\sigma \|_{\ell^2(\Z^{d-1})}^2 \\
\le \sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1 \in \Z} \| u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2
+\lambda_1 \sum_{\sigma=0}^s \Ng f^\sigma \Nd_{-\infty,+\infty}^2.
\end{multline*}
Choosing $\gamma \Delta t$ large enough, that is $\gamma \Delta t \ge \ln R_0$ for some numerical constant
$R_0>1$ that depends only on the (fixed) coefficients of the operator $L$, we can absorb the term
\begin{equation*}
\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t}
\sum_{j_1 \in \Z} \| u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2,
\end{equation*}
from right to left in \eqref{estimaux}, and we have therefore derived the estimate
\begin{multline*}
\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1 \in \Z}
\| u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2\\
\le C \biggl\{
\sum_{n\ge 0} \Delta t \rme^{-2 \gamma (n+s+1) \Delta t} \sum_{j_1 \in \Z}
\| L u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2
+\rme^{-2 \gamma \Delta t} \sum_{\sigma=0}^s \Ng f^\sigma \Nd_{-\infty,+\infty}^2 \biggr\}.
\end{multline*}
It remains to use \eqref{estim1} to bound the first term on the right-hand side, and we get an even better estimate
than \eqref{estimabsorb} which we were originally aiming at:
\begin{multline*}
\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1 \in \Z}
\| u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2\\
\le C \biggl\{ \sum_{\sigma=0}^s \Ng f^\sigma \Nd_{-\infty,+\infty}^2
+\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1 \le 0}
\| g_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \biggl\}.
\end{multline*}
This gives a control of infinitely many traces and not only finitely many (this restriction to finitely many traces
will appear in the regime where $\gamma \Delta t$ can be small).

$\bullet$ From now on, we have fixed a constant $R_0>1$ such that \eqref{estimabsorb} holds for $\gamma
\Delta t \ge \ln R_0$ and we thus assume $\gamma \Delta t \in (0,\ln R_0]$. (Getting rid of all large values of
$\gamma \Delta t$ will be used to gain `compactness'.) We also know that the estimate \eqref{estim1} holds,
independently of the value of $\gamma \Delta t$, and we now wish to estimate the traces of the solution
$(u_j^n)$ for finitely many values of $j_1$.

We first observe from \eqref{estim1} that for all $\gamma>0$ and $\Delta t \in (0,1]$, there exists a constant
$C_{\gamma,\Delta t}$ such that
\begin{equation*}
\forall n \in \N,\quad \rme^{-2 \gamma n \Delta t} \sum_{j \in \Z^d} (u_j^n)^2 \le C_{\gamma,\Delta t}.
\end{equation*}
In particular, for any $j_1 \in \Z$, the Laplace-Fourier transforms $\widehat{u_{j_1}}$ of the step functions
\begin{equation*}
u_{j_1} : (t,y) \in \R^+ \times \R^{d-1} \longmapsto u_j^n \quad \text{\rm if } (t,y) \in \bigl[n \Delta t,(n+1) \Delta t\bigr)
\times \prod_{k=2}^d \bigl[j_k \Delta x_k,(j_k+1) \Delta x_k \bigr),
\end{equation*}
is well-defined on $\{ \tau \in \C, \reel \tau >0 \} \times \R^{d-1}$. The dual variables are denoted
$\tau =\gamma +i \theta$, $\gamma>0$, and $\eta =(\eta_2,\dots,\eta_d) \in \R^{d-1}$. It will also be convenient
to introduce the notation $\eta_\Delta := (\eta_2 \Delta x_2,\dots,\eta_d \Delta x_d)$. Given $\Gamma>0$, the
sequence $(\widehat{u_{j_1}}(\Gamma+i \theta,\eta))_{j_1 \in \Z}$ belongs to $\ell^2(\Z)$ for almost every
$(\theta,\eta) \in \R \times \R^{d-1}$.

We first show the following estimate, which is the Laplace-Fourier analogue of \eqref{estim1}.

\begin{lemma}
\label{lem3'}
With $R_0>1$ fixed as above, there exists a constant $C>0$ such that for all $\gamma>0$ and $\Delta t \in (0,1]$
satisfying $\gamma \Delta t \in (0,\ln R_0]$, there holds
\begin{multline}
\label{estimlem3'}
\sum_{j_1 \in \Z} \int_{\R \times \R^{d-1}} \biggl| \sum_{\ell_1=-r_1}^{p_1}
a_{\ell_1} \big(\rme^{(\gamma +i \theta) \Delta t},\eta_\Delta \big)
\widehat{u_{j_1+\ell_1}}(\gamma +i \theta,\eta) \biggr|^2 \rmd\theta \rmd\eta \\
+\sum_{j_1 \le 0} \int_{\R \times \R^{d-1}} \biggl| \sum_{\ell_1=-r_1}^{p_1} \rme^{(\gamma +i \theta) \Delta t}
\partial_z a_{\ell_1} \big(\rme^{(\gamma +i \theta) \Delta t},\eta_\Delta \big)
\widehat{u_{j_1+\ell_1}}(\gamma +i \theta,\eta) \biggr|^2 \rmd\theta \rmd\eta \\
\le C \biggl\{ \sum_{\sigma=0}^s \Ng f^\sigma \Nd_{-\infty,+\infty}^2
+\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1 \le 0}
\| g_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \biggr\}.
\end{multline}
\end{lemma}

\begin{proof}[Proof of Lemma \ref{lem3'}]
Given $\tau =\gamma +i \theta$ and $\eta$, we compute (here $j_1 \in \Z$ is fixed):
\begin{multline}\label{lem3'1}
\sum_{\ell_1=-r_1}^{p_1} a_{\ell_1} \big(\rme^{\tau \Delta t},\eta_\Delta \big)
\widehat{u_{j_1+\ell_1}}(\tau,\eta)\\[-5pt]
= \widehat{L u_{j_1,\sbullet}} (\tau,\eta)
+\dfrac{1-\rme^{-\tau \Delta t}}{\tau}
\sum_{\sigma=1}^{s+1} \sum_{\sigma'=0}^{\sigma-1} \rme^{(\sigma-\sigma') \tau \Delta t}
{\mathcal F}_{j_1}^{\sigma,\sigma'}(\eta),
\end{multline}
\begin{multline}\label{lem3'2}
\sum_{\ell_1=-r_1}^{p_1} \rme^{\tau \Delta t} \partial_z a_{\ell_1} \big(\rme^{\tau \Delta t},\eta_\Delta \big)
\widehat{u_{j_1+\ell_1}}(\tau,\eta)\\[-5pt]
= \widehat{M u_{j_1,\sbullet}} (\tau,\eta)
+\dfrac{1-\rme^{-\tau \Delta t}}{\tau} \sum_{\sigma=1}^{s+1} \sum_{\sigma'=0}^{\sigma-1} \sigma
\rme^{(\sigma-\sigma') \tau \Delta t} {\mathcal F}_{j_1}^{\sigma,\sigma'}(\eta).
\end{multline}
where, in \eqref{lem3'1} and \eqref{lem3'2}, we have set
\begin{equation*}
{\mathcal F}_{j_1}^{\sigma,\sigma'}(\eta) =\sum_{\ell_1=-r_1}^{p_1} \biggl(\sum_{\ell'=-r'}^{p'} a_{\ell,\sigma}
\rme^{i \ell' \cdot \eta_\Delta} \biggl) \widehat{f^{\sigma'}_{j_1+\ell_1,\sbullet}} (\eta),
\end{equation*}
which corresponds to the partial Fourier transform with respect to $y=(x_2,\dots,x_d) \in \R^{d-1}$, of the step
function associated with the sequence $(Q_\sigma f_j^{\sigma'})$ (no Laplace transform here).

We need to estimate integrals with respect to $(\theta,\eta)$ of the right-hand side of \eqref{lem3'1} and
\eqref{lem3'2}. The first term on the right of \eqref{lem3'1} and \eqref{lem3'2} are easy. For instance, we
have (applying Plancherel Theorem):
\begin{align*}
\sum_{j_1 \in \Z} \int_{\R \times \R^{d-1}} \big| \widehat{L u_{j_1,\sbullet}} (\tau,\eta) \big|^2&
\rmd\theta \rmd\eta =(2 \pi)^d \sum_{j_1 \in \Z} \sum_{n \ge 0} \int_{n \Delta t}^{(n+1) \Delta t}
\rme^{-2 \gamma s} \| L u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \rmd s \\
&= (2 \pi)^d \dfrac{1 -\rme^{-2 \gamma \Delta t}}{2 \gamma \Delta t} \sum_{n \ge 0}
\Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1 \in \Z} \| L u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2.
\end{align*}
We now recall that $\gamma \Delta t$ is restricted to the interval $(0,\ln R_0]$, and we use \eqref{estim1}
to derive
\begin{multline*}
\sum_{j_1 \in \Z} \int_{\R \times \R^{d-1}} \big| \widehat{L u_{j_1,\sbullet}} (\tau,\eta) \big|^2
\rmd\theta \rmd\eta\\
\le C \biggl\{ \sum_{\sigma=0}^s \Ng f^\sigma \Nd_{-\infty,+\infty}^2
+\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1 \le 0}
\| g_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \biggr\}.
\end{multline*}
Similarly, we have
\begin{multline*}
\sum_{j_1 \le 0} \int_{\R \times \R^{d-1}} \big| \widehat{M u_{j_1,\sbullet}} (\tau,\eta) \big|^2
\rmd\theta \rmd\eta\\
=(2 \pi)^d \dfrac{1 -\rme^{-2 \gamma \Delta t}}{2 \gamma \Delta t}
\sum_{n \ge 0} \Delta t \rme^{-2 \gamma n \Delta t}
\sum_{j_1 \le 0} \| M u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2,
\end{multline*}
which we can again uniformly estimate by the right-hand side of \eqref{estimlem3'}.

Going back to the right-hand side terms in \eqref{lem3'1} and \eqref{lem3'2}, we find that there only remains
for proving \eqref{estimlem3'} to estimate the integral (here there are finitely many values of $\sigma$ and
$\sigma'$):
\begin{multline*}
\sum_{j_1 \in \Z} \int_{\R \times \R^{d-1}} \Bigl| \dfrac{1-\rme^{-\tau \Delta t}}{\tau} \Bigr|^2
\big| {\mathcal F}_{j_1}^{\sigma,\sigma'}(\eta) \big|^2 \rmd\theta \rmd\eta\\
= \biggl(\int_\R \Bigl| \dfrac{1-\rme^{-\tau \Delta t}}{\tau} \Bigr|^2 \rmd\theta \biggl)
\sum_{j_1 \in \Z} \int_{\R^{d-1}} \big| {\mathcal F}_{j_1}^{\sigma,\sigma'}(\eta) \big|^2 \rmd\eta,
\end{multline*}
where we have applied Fubini Theorem. Applying first Plancherel Theorem with respect to the $d-1$ last space
variables, we get
\begin{equation*}
\sum_{j_1 \in \Z} \int_{\R^{d-1}} \big| {\mathcal F}_{j_1}^{\sigma,\sigma'}(\eta) \big|^2 \rmd\eta \le C
\sum_{j_1 \in \Z} \sum_{j' \in \Z^{d-1}} \biggl(\prod_{k=2}^d \Delta x_k \biggr) (f_j^{\sigma'})^2 \le
\dfrac{C}{\Delta t} \sum_{\sigma=0}^s \Ng f^\sigma \Nd_{-\infty,+\infty}^2.
\end{equation*}
The conclusion then follows by computing
\begin{equation*}
\int_\R \Bigl| \dfrac{1 -\rme^{-\tau \Delta t}}{\tau} \Bigr|^2 \rmd\theta =
2 \pi \Delta t \dfrac{1 -\rme^{-2 \gamma \Delta t}}{2 \gamma \Delta t},
\end{equation*}
and by recalling that $\gamma \Delta t$ belongs to $(0,\ln R_0]$. We can eventually bound the integrals on
the left-hand side of \eqref{estimlem3'} by estimating separately the integrals of each term on the right-hand side
of \eqref{lem3'1} and \eqref{lem3'2}.
\end{proof}

\noindent The conclusion in the proof of Theorem \ref{absorbing} relies on the following crucial result. Here both
Assumptions \ref{assumption1} and \ref{assumption2} are heavily used.

\begin{lemma}[The trace estimate]
\label{lem3}
Let Assumptions \ref{assumption0}, \ref{assumption1} and \ref{assumption2} be satisfied. Let $R_0>1$ be fixed
as above and let $P_1 \in \N$. Then there exists a constant $C_{P_1}>0$ such that for all $z \in \U$ with $|z| \le
R_0$, for all $\eta \in \R^{d-1}$ and for all sequence $(w_{j_1})_{j_1 \in \Z} \in \ell^2(\Z;\C)$, there holds
\begin{multline}
\label{estimlem3}
\sum_{j_1=-r_1-p_1}^{P_1} |w_{j_1}|^2\\
\le C_{P_1} \biggl\{ \sum_{j_1 \in \Z} \biggl| \sum_{\ell_1=-r_1}^{p_1}
a_{\ell_1}(z,\eta) w_{j_1+\ell_1} \biggr|^2
+\sum_{j_1 \le 0} \biggl| \sum_{\ell_1=-r_1}^{p_1} z \partial_z a_{\ell_1}(z,\eta) w_{j_1+\ell_1} \biggr|^2
\biggr\}.
\end{multline}
Recall that the functions $a_{\ell_1}$, $\ell_1=-r_1,\dots,p_1$, are defined in \eqref{defA-d}.
\end{lemma}

The proof of Lemma \ref{lem3} is rather long. Before giving it in full details, we indicate how Lemma \ref{lem3}
yields the result of Theorem \ref{absorbing}. We apply Lemma \ref{lem3} to $z=\exp (\tau \Delta t)$,
$\tau=\gamma +i \theta$ with $\gamma \Delta t \in (0, \ln R_0]$, $\eta_\Delta \in \R^{d-1}$ and
the sequence $(\widehat{u_{j_1}} (\tau,\eta))_{j_1 \in \Z}$. We then integrate \eqref{estimlem3} with respect to
$(\theta,\eta)$ and use Lemma \ref{lem3'} to derive
\begin{multline*}
\sum_{j_1=-r_1-p_1}^{P_1} \int_{\R \times \R^{d-1}} \left| \widehat{u_{j_1}} (\gamma+i \theta,\eta) \right|^2
\rmd\theta \rmd\eta\\
\le C \biggl\{ \sum_{\sigma=0}^s \Ng f^\sigma \Nd_{-\infty,+\infty}^2
+\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1 \le 0}
\| g_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \biggl\}.
\end{multline*}
It remains to apply Plancherel Theorem and we get
\begin{multline*}
\dfrac{1-\rme^{-2 \gamma \Delta t}}{2 \gamma \Delta t} \sum_{j_1=-r_1-p_1}^{P_1} \sum_{n \in \N}
\Delta t \rme^{-2 \gamma n \Delta t} \| u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \\
\le C \biggl\{ \sum_{\sigma=0}^s \Ng f^\sigma \Nd_{-\infty,+\infty}^2
+\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1 \le 0}
\| g_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \biggl\}.
\end{multline*}
Recalling that $\gamma \Delta t$ is restricted to the interval $(0, \ln R_0]$, we have thus derived the trace
estimate
\begin{multline*}
\sum_{n \in \N} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1=-r_1-p_1}^{P_1}
\| u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2\\
\le C \biggl\{ \sum_{\sigma=0}^s \Ng f^\sigma \Nd_{-\infty,+\infty}^2
+\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1 \le 0}
\| g_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \biggl\}.
\end{multline*}
Combined with the semigroup and interior estimate \eqref{estim1}, this gives the estimate \eqref{estimabsorb}
of Theorem \ref{absorbing} for $\gamma \Delta t \in (0,\ln R_0]$.
\end{proof}

\begin{proof}[Proof of Lemma \ref{lem3}]
Let us recall that the functions $a_{\ell_1}$ are $2 \pi$-periodic with respect to each coordinate of $\eta$.
We can therefore restrict to $\eta \in [0,2 \pi]^{d-1}$ rather than considering $\eta \in \R^{d-1}$. We argue
by contradiction and assume that the conclusion to Lemma \ref{lem3} does not hold. This means the following,
up to normalizing and extracting subsequences; there exist three sequences (indexed by $k \in \N$):
\begin{itemize}
\item a sequence $(w^k)_{k \in \N}$ with values in $\ell^2(\Z;\C)$ such that $(w_{-r_1-p_1}^k,\dots,w_{P_1}^k)$
belongs to the unit sphere of $\C^{P_1+r_1+p_1+1}$ for all $k$, and $(w_{-r_1-p_1}^k,\dots,w_{P_1}^k)$
converges to $(\underline{w}_{-r_1-p_1},\dots,\underline{w}_{P_1})$ as $k$ tends to infinity,

\item a sequence $(z^k)_{k \in \N}$ with values in $\U \cap \{ \zeta \in \C, |\zeta| \le R_0 \}$, which converges to $\underline{z} \in \Ubar$,

\item a sequence $(\eta^k)_{k \in \N}$ with values in $[0,2 \pi]^{d-1}$, which converges to $\underline{\eta}
\in [0,2 \pi]^{d-1}$,
\end{itemize}
and these sequences satisfy:
\begin{equation}
\label{lem3-1}
\lim_{k \rightarrow +\infty} \sum_{j_1 \in \Z} \biggl| \sum_{\ell_1=-r_1}^{p_1}
a_{\ell_1}(z^k,\eta^k) w^k_{j_1+\ell_1} \biggr|^2
+\sum_{j_1 \le 0} \biggl| \sum_{\ell_1=-r_1}^{p_1} z^k \partial_z a_{\ell_1}(z^k,\eta^k) w^k_{j_1+\ell_1}
\biggr|^2 =0.
\end{equation}
We are going to show that \eqref{lem3-1} implies that $(\underline{w}_{-r_1-p_1},\dots,\underline{w}_{P_1})$
must be zero, which will yield a contradiction since this vector must have norm $1$.

$\bullet$ Let us first show that each component $(w^k_{j_1})_{k \in \N}$, $j_1 \in \Z$, has a limit as $k$ tends
to infinity. This is already clear for $j_1=-r_1-p_1,\dots,P_1$. For $j_1>P_1$, we argue by induction. From
\eqref{lem3-1}, we have
\begin{equation*}
\lim_{k \rightarrow +\infty} \sum_{\ell_1=-r_1}^{p_1} a_{\ell_1}(z^k,\eta^k) w^k_{P_1-p_1+1+\ell_1} =0,
\end{equation*}
and by Assumption \ref{assumption2}, we know that $a_{p_1}(\underline{z},\underline{\eta})$ is nonzero.
Hence $(w^k_{P_1+1})_{k \in \N}$ converges towards
\begin{equation*}
-\dfrac{1}{a_{p_1}(\underline{z},\underline{\eta})} \sum_{\ell_1=-r_1}^{p_1-1} a_{\ell_1}(\underline{z},\underline{\eta})
\underline{w}_{P_1-p_1+1+\ell_1},
\end{equation*}
which we define as $\underline{w}_{P_1+1}$. We can argue by induction in the same way for all indices
$j_1>P_1+1$, but also for indices $j_1<-r_1-p_1$ because the function $a_{-r_1}$ also does not vanish
on $\Ubar \times \R^{d-1}$.

Using \eqref{lem3-1}, we have thus shown that for each $j_1 \in \Z$, $(w^k_{j_1})_{k \in \N}$ tends towards
some limit $\underline{w}_{j_1}$ as $k$ tends to infinity, and the sequence $\underline{w}$, which does not
necessarily belong to $\ell^2(\Z;\C)$, satisfies the induction relations:
\begin{align}
\forall j_1 \in \Z,\quad & \sum_{\ell_1=-r_1}^{p_1} a_{\ell_1}(\underline{z},\underline{\eta})
\underline{w}_{j_1+\ell_1} =0,\label{induction1}\\
\forall j_1 \le 0,\quad & \sum_{\ell_1=-r_1}^{p_1} \underline{z} \partial_z a_{\ell_1}(\underline{z},\underline{\eta})
\underline{w}_{j_1+\ell_1} =0.\label{induction2}
\end{align}

$\bullet$ The induction relation \eqref{induction1} is the one that arises in \cite{gks,michelson} and all the works that
deal with strong stability. The main novelty here is to use simultaneously \eqref{induction1} for controlling the unstable
components of $(\underline{w}_{-r_1-p_1},\dots,\underline{w}_{-1})$ and \eqref{induction2} for controlling the stable
components of $(\underline{w}_{-r_1-p_1},\dots,\underline{w}_{-1})$. The fact that $\underline{w}$ satisfies
simultaneously \eqref{induction1} and \eqref{induction2} for $j_1 \le 0$ automatically annihilates the central
components. This sketch of proof is made precise below.

We define the source terms:
\begin{equation*}
F_{j_1}^k := \sum_{\ell_1=-r_1}^{p_1} a_{\ell_1}(z^k,\eta^k) w^k_{j_1+\ell_1},\quad
G_{j_1}^k := \sum_{\ell_1=-r_1}^{p_1} z^k \partial_z a_{\ell_1}(z^k,\eta^k) w^k_{j_1+\ell_1},
\end{equation*}
which, according to \eqref{lem3-1}, satisfy
\begin{equation}
\label{lem3-2}
\lim_{k \rightarrow 0} \sum_{j_1 \in \Z} |F_{j_1}^k|^2 =0,\quad
\lim_{k \rightarrow 0} \sum_{j_1 \le 0} |G_{j_1}^k|^2 =0.
\end{equation}
We also introduce the vectors (here $T$ denotes transposition)
\begin{equation*}
W_{j_1}^k := \bigl(w^k_{j_1+p_1},\dots,w^k_{j_1+1-r_1} \bigr)^T,\quad
\underline{W}_{j_1} := \bigl(\underline{w}_{j_1+p_1},\dots,\underline{w}_{j_1+1-r_1} \bigr)^T,
\end{equation*}
and the matrices:
\begin{align}
\LL (z,\eta) &:= \begin{pmatrix}
-\dfrac{a_{p_1-1} (z,\eta)}{a_{p_1} (z,\eta)} & & \hspace*{3mm}& \dots &\hspace*{5mm}& -\dfrac{a_{-r_1} (z,\eta)}{a_{p_1} (z,\eta)} \\[10pt]
\hspace*{6mm}1& 0 && \dots &&\hspace*{-6mm}0 \\
& \ddots&& \ddots& & \\
\hspace*{6mm}0& && \ddots& \hspace*{2mm}\ddots&\hspace*{-6mm}\vdots \\
\hspace*{6mm}0& 0 && &1&\hspace*{-6mm}0\end{pmatrix} \in {\mathcal M}_{p_1+r_1}(\C),\label{defL} \\[10pt]
\M (z,\eta) &:= \begin{pmatrix}
-\dfrac{\partial_z a_{p_1-1} (z,\eta)}{\partial_z a_{p_1} (z,\eta)} & & \hspace*{3mm}& \dots &\hspace*{5mm}&
-\dfrac{\partial_z a_{-r_1} (z,\eta)}{\partial_z a_{p_1} (z,\eta)} \\[10pt]
\hspace*{6mm}1& 0 && \dots &&\hspace*{-6mm}0 \\
& \ddots&& \ddots& & \\
\hspace*{6mm}0& && \ddots& \hspace*{2mm}\ddots&\hspace*{-6mm}\vdots \\
\hspace*{6mm}0& 0 && &1&\hspace*{-6mm}0\end{pmatrix} \in {\mathcal M}_{p_1+r_1}(\C).\label{defM}
\end{align}
The matrix $\LL$ is well-defined on $\Ubar \times \R^{d-1}$ according to Assumption \ref{assumption2}. The
matrix~$\M$ is also well-defined on $\Ubar \times \R^{d-1}$ because for any $\eta \in \R^{d-1}$, Assumption
\ref{assumption2} asserts that $a_{p_1}(\cdot,\eta)$ is a non-constant polynomial whose roots lie in $\D$. From
the Gauss-Lucas Theorem, the roots of $\partial_z a_{p_1}(\cdot,\eta)$ lie in the convex hull of those of
$a_{p_1}(\cdot,\eta)$. Therefore $\partial_z a_{p_1}(\cdot,\eta)$ does not vanish on $\Ubar$. In the same
way, $\partial_z a_{-r_1}(\cdot,\eta)$ does not vanish on~$\Ubar$.

With our above notation, the vectors $W_{j_1}^k$, $\underline{W}_{j_1}$, satisfy the one step induction relations:
\begin{align}\label{induction1'}
&\forall j_1 \in \Z,\quad \begin{aligned}
W_{j_1+1}^k &= \LL(z^k,\eta^k) W_{j_1}^k
+ \bigl(F^k_{j_1+1}/a_{p_1} (z^k,\eta^k),0,\dots,0 \bigr)^T,\\
\underline{W}_{j_1+1} &=\LL(\underline{z},\underline{\eta}) \underline{W}_{j_1},
\end{aligned} \\
\label{induction2'}
&\forall j_1 \le -1,\quad \begin{aligned}
W_{j_1+1}^k &= \M(z^k,\eta^k) W_{j_1}^k
+\bigl(G^k_{j_1+1}/(z^k \partial_z a_{p_1} (z^k,\eta^k)),0,\dots,0 \bigr)^T,\\
\underline{W}_{j_1+1} &=\M(\underline{z},\underline{\eta}) \underline{W}_{j_1}.
\end{aligned}
\end{align}

$\bullet$ From Assumption \ref{assumption2} and the above application of the Gauss-Lucas Theorem, we
already know that both matrices $\LL(z,\eta)$ and $\M(z,\eta)$ are invertible for $(z,\eta) \in \Ubar \times
\R^{d-1}$. Furthermore, Assumption \ref{assumption1} shows that $\LL(z,\eta)$ has no eigenvalue on
$\cercle$ for $(z,\eta) \in \U \times \R^{d-1}$. This property dates back at least to \cite{kreiss1}. However,
central eigenvalues on $\cercle$ may occur for $\LL$ when $z$ belongs to $\cercle$. The crucial point for
proving Lemma \ref{lem3} is that Assumption \ref{assumption1} precludes central eigenvalues of $\M$ for
all $z \in \Ubar$. Namely, for all $z \in \Ubar$ and all $\eta \in \R^{d-1}$, $\M(z,\eta)$ has no eigenvalue on
$\cercle$. This property holds because otherwise, for some $(z,\eta) \in \Ubar \times \R^{d-1}$, there would
exist a solution $\kappa_1 \in \cercle$ to the dispersion relation
\begin{equation*}
\sum_{\ell_1=-r_1}^{p_1} z \partial_z a_{\ell_1}(z,\eta) \kappa_1^{\ell_1} =0.
\end{equation*}
For convenience, the coordinates of $\eta$ are denoted $(\eta_2,\dots,\eta_d)$. Using the definition \eqref{defA-d}
of $a_{\ell_1}$, and defining $\kappa :=(\kappa_1,\rme^{i \eta_2},\dots,\rme^{i \eta_d})$, we have found a
root $z \in \Ubar$ to the relation
\begin{equation}
\label{dispersionder}
\sum_{\sigma=1}^{s+1} \sigma \widehat{Q_\sigma} (\kappa) z^{\sigma-1} =0,
\end{equation}
but this is not possible because the $s+1$ roots (in $z$) to the dispersion relation \eqref{dispersion} are simple
and belong to $\Dbar$. The Gauss-Lucas Theorem thus shows that the roots to the relation \eqref{dispersionder}
belong to $\D$ (and therefore not to $\Ubar$).

At this stage, we know that the eigenvalues of $\M(z,\eta)$, $(z,\eta) \in \Ubar \times \R^{d-1}$, 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,\eta) \in \Ubar \times \R^{d-1}$, we then introduce the spectral projector $\Pi_\M^s(z,\eta)$, resp. $\Pi_\M^u(z,\eta)$,
of $\M(z,\eta)$ on the generalized eigenspace associated with eigenvalues in $\D$, resp.\ $\U$. We can then integrate
the first induction relation in \eqref{induction2'} and get
\begin{equation*}
\Pi_\M^s(z^k,\eta^k) W_0^k =\dfrac{1}{z^k \partial_z a_{p_1} (z^k,\eta^k)} \sum_{j_1 \le 0}
\M(z^k,\eta^k)^{|j_1|} \Pi_\M^s(z^k,\eta^k) \bigl(G^k_{j_1},0,\dots,0 \bigr)^T.
\end{equation*}
The projector $\Pi_\M^s$ depends continuously on $(z,\eta) \in \Ubar \times \R^{d-1}$. Furthermore, since the
spectrum of $\M$ does not meet $\cercle$ even for $z \in \cercle$, there exists a constant $C>0$ and a $\delta \in
(0,1)$ that are independent of $k \in \N$ and such that
\begin{equation*}
\forall j_1 \le 0,\quad \bigl\| \M(z^k,\eta^k)^{|j_1|} \Pi_\M^s(z^k,\eta^k) \bigr\| \le C \delta^{|j_1|}.
\end{equation*}
We thus get a uniform estimate with respect to $k$:
\begin{equation*}
\bigl|\Pi_\M^s(z^k,\eta^k) W_0^k \bigr|^2 \le C \sum_{j_1 \le 0} |G^k_{j_1}|^2.
\end{equation*}
Passing to the limit and using \eqref{lem3-2}, we get $\Pi_\M^s (\underline{z},\underline{\eta}) \underline{W}_0 =0$,
or in other words $\underline{W}_0 =\Pi_\M^u (\underline{z},\underline{\eta}) \underline{W}_0$.

$\bullet$ The sequence $(\underline{W}_{j_1})_{j_1 \le 0}$ satisfies both induction relations \eqref{induction1'} and
\eqref{induction2'}. Due to the form of the companion matrices $\LL$ and $\M$, see \eqref{defL}-\eqref{defM},
we can conclude that the vector $\underline{W}_0$ belongs to the generalized eigenspace (of either $\LL$ or
$\M$) associated with the common eigenvalues of $\M(\underline{z},\underline{\eta})$ and $\LL(\underline{z},
\underline{\eta})$. We have already seen that $\M(\underline{z},\underline{\eta})$ has no eigenvalue on
$\cercle$ and $\underline{W}_0 =\Pi_\M^u (\underline{z},\underline{\eta}) \underline{W}_0$, so we can
conclude that $\underline{W}_0$ belongs to the generalized eigenspace of $\LL$ associated with those
common eigenvalues of $\M(\underline{z},\underline{\eta})$ and $\LL(\underline{z},\underline{\eta})$ in
$\U$.

The matrix $\LL(\underline{z},\underline{\eta})$ has $N^u$ eigenvalues in $\U$, $N^s$ in $\D$ and $N^c$ on
$\cercle$. (Since $\underline{z}$ may belong to $\cercle$, $N^c$ is not necessarily zero.) With obvious notations,
we let $\Pi_\LL^{u,s,c}(z,\eta)$ denote the corresponding spectral projectors of $\LL$ for $(z,\eta)$ sufficiently
close to $(\underline{z},\underline{\eta})$. In particular, the eigenvalues corresponding to $\Pi_\LL^u(z,\eta)$
lie in $\U$ uniformly away from~$\cercle$ for $(z,\eta)$ sufficiently close to $(\underline{z},\underline{\eta})$.
We can then integrate the first induction relation in \eqref{induction1'} and derive (for $k$ sufficiently large):
\begin{equation*}
\Pi_\LL^u(z^k,\eta^k) W_0^k =-\dfrac{1}{a_{p_1} (z^k,\eta^k)} \sum_{j_1 \ge 0}
\LL(z^k,\eta^k)^{-j_1-1} \Pi_\LL^u(z^k,\eta^k) \bigl(F^k_{j_1},0,\dots,0 \bigr)^T.
\end{equation*}
Using the uniform exponential decay of $\LL(z^k,\eta^k)^{-j_1-1} \Pi_\LL^u(z^k,\eta^k)$ and \eqref{lem3-2},
we finally end up with
\begin{equation*}
\Pi_\LL^u(\underline{z},\underline{\eta}) \underline{W}_0=0.
\end{equation*}
Since $\underline{W}_0$ belongs to the generalized eigenspace of $\LL$ associated with those common
eigenvalues of $\M(\underline{z},\underline{\eta})$ and $\LL(\underline{z},\underline{\eta})$ in $\U$, we can
conclude that $\underline{W}_0$ equals zero. Applying the induction relation \eqref{induction1'}, the whole
sequence $(\underline{W}_{j_1})_{j_1 \in \Z}$ is zero which yields the expected contradiction.
\end{proof}

The crucial property that we use in the proof of Lemma \ref{lem3} is the fact that up to $z \in \cercle$, the
eigenvalues of $\M(z,\eta)$ lie either in $\D$ or $\U$. For the leap-frog scheme, this property would not be
true if we had imposed the auxiliary numerical boundary condition $u_j^{n+2}+u_j^n$ rather than $2 u_j^{n+2}
+\lambda a (u_{j+1}^{n+1}-u_{j-1}^{n+1})$.

Let us also observe that we have used the fact that $a_{p_1}$ and $a_{-r_1}$ are non-constant in order to
study the induction relation \eqref{induction2}. There might be some schemes for which~$a_{p_1}$ and/or
$a_{-r_1}$ are constant but for which one can still apply similar arguments as in the previous proof, even
though \eqref{induction2} is an induction relation with fewer steps than \eqref{induction1}. In this respect,
Assumption \ref{assumption2} might be relaxed in specific applications.

\begin{remark}
The auxiliary problem \eqref{numabsorbing} is in general not of the same form as \eqref{numibvp} because
in \eqref{numabsorbing} one has to impose infinitely many numerical boundary conditions. This is due to the
fact that the stencil of $M$ incorporates points `on the left' with respect to the first space variable. A remarkable
exception occurs for explicit schemes with $s=0$, for in that case the multiplier $M v_j^n$ reads $v_j^{n+1}$
and \eqref{numabsorbing} is exactly the auxiliary problem considered in \cite{jfcag} (and labeled (2.7) there)
where one imposes Dirichlet boundary conditions on finitely many boundary meshes (just use $g_j^n =0$ for
$j_1 \le -r_1$). In full generality, there still remains an open problem of constructing a set of dissipative
numerical boundary conditions of the same form as \eqref{numibvp} with $s \ge 1$, that is with finitely
many numerical boundary conditions, and for which one can prove by hand both a semigroup and a trace
estimate as in Theorem \ref{absorbing}.
\end{remark}

\subsection{End of the proof}

As explained in the introduction of Section \ref{section3}, the linearity of \eqref{numibvp} reduces the proof of
Theorem \ref{mainthm} to the case $(F_j^n)=0$, $(g_j^n)=0$, since we have already dealt with the case of zero
initial data. We thus focus on \eqref{numibvp} with $(F_j^n)=0$ and $(g_j^n)=0$, and write the corresponding
solution $(u_j^n)$ as $u_j^n =v_j^n +w_j^n$, where the sequence $(v_j^n)$ solves:
\begin{equation}
\label{defvjn}
\begin{cases}
L v_j^n =0,& j_1 \ge 1,\quad j' \in \Z^{d-1},\quad n\ge 0,\\
M v_j^n =0,& j_1 \le 0,\quad j' \in \Z^{d-1},\quad n\ge 0,\\
v_j^n = f_j^n,& j \in \Z^d,\quad n=0,\dots,s,
\end{cases}
\end{equation}
and $(w_j^n)$ solves:
\begin{equation}
\label{defwjn}
\begin{cases}
L w_j^n =0,& j \in \Z^d,\ j_1 \ge 1,\ n\ge 0,\\
w_j^{n+s+1} +{\dps \sum_{\sigma=0}^{s+1}} B_{j_1,\sigma} w_{1,j'}^{n+\sigma} =\tilde{g}_j^{n+s+1},&
j \in \Z^d,\ j_1=1-r_1,\dots,0,\ n\ge 0,\\
w_j^n = 0,& j \in \Z^d,\ n=0,\dots,s.
\end{cases}
\end{equation}
For $(v_j^n +w_j^n)_{j_1 \ge 1-r_1}$ to coincide with the solution $(u_j^n)$ to \eqref{numibvp}, it is sufficient to extend
the initial data $f_j^0,\dots,f_j^s$ by zero for $j_1 \le -r_1$, which provides with the initial data in \eqref{defvjn} on all
$\Z^d$, and to define the boundary source term in \eqref{defwjn} by:
\begin{equation}
\label{defgtilde}
\tilde{g}_j^{n+s+1} := -v_j^{n+s+1} -{\dps \sum_{\sigma=0}^{s+1}} B_{j_1,\sigma} v_{1,j'}^{n+\sigma}.
\end{equation}

We can estimate the solution $(v_j^n)$ to \eqref{defvjn} by applying Theorem \ref{absorbing}. In particular,
the trace estimate:
\begin{equation*}
\sum_{n\ge 0} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1=1-r_1}^{P_1}
\| v_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \le C \sum_{\sigma=0}^s \Ng f^\sigma \Nd_{1-r_1,+\infty}^2,
\end{equation*}
for $P_1=\max(p_1,q_1+1)$ gives (recall the definition \eqref{defgtilde} of $\tilde{g}_j^{n+s+1}$):
\begin{align*}
\sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1=1-r_1}^0
\| \tilde{g}_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2
&\le C \sum_{n\ge 0} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1=1-r_1}^{\max(p_1,q_1+1)}
\| v_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \\
&\le C \sum_{\sigma=0}^s \Ng f^\sigma \Nd_{1-r_1,+\infty}^2.
\end{align*}
We can apply Theorem \ref{mainthm} to the solution $(w_j^n)$ to \eqref{defwjn} because the initial data
in \eqref{defwjn} vanish. We get:
\begin{multline*}
\sup_{n \ge 0} \rme^{-2 \gamma n \Delta t} \Ng w^n \Nd_{1-r_1,+\infty}^2
+\dfrac{\gamma}{\gamma \Delta t+1}
\sum_{n\ge 0} \Delta t \rme^{-2 \gamma n \Delta t} \Ng w^n \Nd_{1-r_1,+\infty}^2 \\[-8pt]
\shoveright{+\sum_{n\ge 0} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1=1-r_1}^{p_1}
\| w_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2}\\
\le C \sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t}
\sum_{j_1=1-r_1}^0 \| \tilde{g}_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 
\le C \sum_{\sigma=0}^s \Ng f^\sigma \Nd_{1-r_1,+\infty}^2.
\end{multline*}
Combining with the similar estimate provided by Theorem \ref{absorbing} for $(v_j^n)$, we end up with
the expected estimate:
\begin{multline*}
\sup_{n \ge 0} \rme^{-2 \gamma n \Delta t} \Ng u^n \Nd_{1-r_1,+\infty}^2
+\dfrac{\gamma}{\gamma \Delta t+1}
\sum_{n\ge 0} \Delta t \rme^{-2 \gamma n \Delta t} \Ng u^n \Nd_{1-r_1,+\infty}^2 \\
+\sum_{n\ge 0} \Delta t \rme^{-2 \gamma n \Delta t} \sum_{j_1=1-r_1}^{p_1}
\| u_{j_1,\sbullet}^n \|_{\ell^2(\Z^{d-1})}^2 \le C \sum_{\sigma=0}^s \Ng f^\sigma \Nd_{1-r_1,+\infty}^2,
\end{multline*}
which completes the proof of Theorem \ref{mainthm}.

\section{Conclusion and perspectives}

Let us first observe that in \cite{wade}, Wade has constructed symmetrizers for deriving stability
estimates for multistep schemes, even in the case of variable coefficients. His conditions for constructing
a symmetrizer are less restrictive than Assumption \ref{assumption1}. However, the symmetrizer in
\cite{wade} is genuinely nonlocal and it is therefore not clear that it may be useful for boundary value
problems. The main novelty here is to construct a {\it local} multiplier whose properties allow for the
design of an auxiliary {\it dissipative} boundary value problem. This is the key to Theorem \ref{mainthm},
despite the nonlocal feature of our energy functional.

The main possible improvement of Theorem \ref{mainthm} would consist of assuming that only the roots to
\eqref{dispersion} that lie on $\cercle$ are simple. Here we have assumed that all the roots, including those
in $\D$ are simple. If we could manage to deal with multiple roots in $\D$, then Theorem \ref{mainthm} would
be applicable to any stable numerical approximation of the transport equation \eqref{transport} (recall that
uniform power boundedness for the amplification matrix ${\mathcal A}$ given in \eqref{defA2pas} requires
only that eigenvalues of modulus $1$ be simple).

The results in this paper achieve the proof of a `weak form' of the conjecture in \cite{kreiss-wu} that strong
stability, in the sense of Definition \ref{defstab1}, implies semigroup stability. However, an even stronger
assumption was made in \cite{kreiss-wu}, namely that the sole fulfillment of the interior estimate
\begin{equation*}
\dfrac{\gamma}{\gamma \Delta t+1} \sum_{n\ge s+1} \Delta t \rme^{-2 \gamma n \Delta t}
\Ng u^n \Nd_{1-r_1,+\infty}^2 \le C \dfrac{\gamma \Delta t+1}{\gamma} \sum_{n\ge s+1} \Delta t
\rme^{-2 \gamma n \Delta t} \Ng F^n \Nd_{1,+\infty}^2,
\end{equation*}
when {\it both initial and boundary data} for \eqref{numibvp} vanish, does imply semigroup stability.
The analogous conjecture for partial differential equations seems to be still open so far, but we do hope
that our multiplier technique may yield some insight for dealing with the strong form of the conjecture in
\cite{kreiss-wu}. We also hope to extend our multiplier technique to prove some stability estimates for some
multistep finite volume schemes on non-Cartesian meshes.

\vspace*{\baselineskip}%
\enlargethispage{-\baselineskip}%
\backmatter
\bibliographystyle{jepplain}
\bibliography{coulombel}
\end{document}