The Neumann numerical boundary condition for transport equations

In this article, we show that prescribing homogeneous Neumann type numerical boundary conditions at an outflow boundary yields a convergent discretization in $\ell^\infty$ for transport equations. We show in particular that the Neumann numerical boundary condition is a stable, local, and absorbing numerical boundary condition for discretized transport equations. Our main result is proved for explicit two time level numerical approximations of transport operators with arbitrarily wide stencils. The proof is based on the energy method and bypasses any normal mode analysis.


Introduction
It is a well-known fact that transport equations do not require prescription of any boundary condition at an outflow boundary, that is, when the transport velocity is outgoing with respect to the boundary of the spatial domain. This can be easily understood by integrating the equation along the characteristics. However, many discretizations of the transport equation involve a stencil that includes cells of the numerical grid that are located in the downstream region. Such discretizations necessitate the prescription of numerical boundary conditions at an outflow boundary [GKO95,Str04], even though the underlying partial differential operator does not require any boundary condition for determining the solution.
The construction and analysis of transparent versus absorbing numerical boundary conditions for wave propagation problems now has a very long history, going back at least to the fundamental contribution by Engquist and Majda [EM77], see also among numerous other works [Hag99, Hal82, Hig86, Hig94, GKO95, AAB + 08] and references therein. Our goal in this article is to rigorously justify the commonly acknowledged fact that enforcing Neumann type boundary conditions at an outflow boundary "does the job", in the simple one-dimensional case with a constant velocity (both on the half-line and on an interval). In several respects, the result is definitely not new. Stability estimates for Neumann numerical boundary conditions were stated -though without proof-for instance by Kreiss [Kre66], and a detailed proof of the latter result was later provided by Goldberg [Gol77]. This result now enters the more general framework of [GT81]. However, the approach in [Kre66,Gol77,GT81] is rather elaborate since it relies first on the verification of the so-called Uniform Kreiss-Lopatinskii Condition (that is, in the present context of numerical schemes, a refined version of the Godunov-Ryabenkii condition [GKS72,GKO95]), and then on the application of deep general results which show that the latter condition is sufficient for the derivation of optimal semigroup estimates. Such general results first arose in [Kre66] and were later proved in further generality by Kreiss, Osher and followers [Kre68,Osh69,Wu95,CG11]. When combined with the trace estimates provided by the fulfillment of the Uniform Kreiss-Lopatinskii Condition and the general convergence result of Gustafsson [Gus75], one gets a complete -though lengthy (!) and somehow unclear as far as the topology is concerned-justification that enforcing Neumann type numerical boundary conditions at an outflow boundary yields a convergent scheme for discretized transport equations.
Our goal here is to bypass all the arguments of those previous works that were based on the normal mode analysis and to obtain a more direct convergence result in the ℓ ∞ topology with arguments that are as elementary as possible. Our approach is based on the energy method and discrete integration by parts rather than on the Laplace transform. We hope that some of our arguments might be useful to deal with more involved problems such as multidimensional hyperbolic systems. The main gain when enforcing the Neumann numerical boundary condition with respect, for instance, to the Dirichlet numerical boundary condition is to obtain convergence results in ℓ ∞ , while in the case of the Dirichlet boundary condition at an outflow boundary, despite unconditional stability properties [GT81], boundary layer phenomena allow at best for convergence results in ℓ p , p < +∞, only [KL68,CHG01,BC17]. Our main result below gives a rate of convergence in ℓ ∞ for such discretizations of the transport equations. We do not claim that our rate is optimal, but we do not assume the discretization of the transport equation to be stable in ℓ ∞ (Z) either, so if we do not reach optimality we are certainly not too far from it. We plan to study the more favorable case of ℓ ∞ stable schemes in the future, with the aim of improving the rate of convergence.
The rest of this article is organized as follows. In Section 2, we introduce some notation and state our main convergence result for Neumann numerical boundary conditions at an outflow boundary. Based on the standard approach of numerical analysis, our convergence result relies on accurate stability estimates and a consistency analysis. Our main stability estimate, based on a new elementary approach, is stated and proved in Section 3. (A crucial discrete integration by parts lemma is given in Appendix A.) The concluding arguments for proving our main result are given in Section 4.

Main result
In this article, we consider the following transport problem on an interval. We are given a fixed constant velocity a > 0, an interval length L > 0 and consider the problem:      ∂ t u + a ∂ x u = 0 , t ≥ 0 , x ∈ (0, L) , u(0, x) = u 0 (x) , x ∈ (0, L) , u(t, 0) = 0 , t ≥ 0 , (2.1) with, at least, u 0 ∈ L 2 ((0, L)). Actually, further regularity and compatibility requirements will be enforced later on, but let us stick to that simple framework for the moment. We could consider an inhomogeneous Dirichlet boundary condition at x = 0 in (2.1), but we rather stick to this simpler case in order to highlight the main novelty of our approach which rather focuses on the downstream boundary x = L of the interval. The solution to (2.1) is given by the method of characteristics, which yields the explicit representation formula: ∀ (t, x) ∈ R + × (0, L) , u(t, x) = u 0 (x − a t) , where it is understood in (2.2) that the initial condition u 0 has been extended by zero to R − (no extension is needed on (L, +∞) since a is positive and therefore x − a t < L for all relevant values of t and x). It may seem a too much trivial problem to approximate the problem (2.1) for which an explicit solution is given, but one should keep in mind that such a representation formula ceases to be available for hyperbolic systems in several space dimensions, and our goal is to develop analytical tools which do not rely on the fact that (2.1) is a one-dimensional scalar problem. We therefore consider from now on an approximation of (2.1) by means of a finite difference scheme. We are given a positive integer J, that is meant to be large, and define accordingly the space step ∆x and the grid points (x j ) by: The time step ∆t is then defined as ∆t := λ ∆x where λ > 0 is a fixed constant that is tuned so that our main Assumption 2.1 below is satisfied. The interval (0, L) is divided in J cells (x j−1 , x j ), j = 1, . . . , J, as depicted in Figure 2.1. In what follows, we use the notation t n := n ∆t, n ∈ N, and u n j will play the role of an approximation for the solution u to (2.1) at time t n on the cell (x j−1 , x j ) (or at the mid-point (x j−1 + x j )/2). We do not wish to discriminate between finite difference or finite volume schemes for (2.1), so rather than deriving this or that type of numerical scheme, we consider a linear iteration for the u n j 's that reads in the interior domain: In (2.3), r, p are fixed nonnegative integers, and the coefficients a ℓ , ℓ = −r, . . . , p may only depend on the ratio λ and the velocity a. Most of the usual linear explicit schemes, such as the upwind, Rusanov, Lax-Friedrichs and Lax-Wendroff schemes, can be put in that form. The interior domain corresponds to the indices j = 1, . . . , J. However, in (2.3), the determination of (u n+1 1 , . . . , u n+1 J ) requires the prior knowledge of (u n 1−r , . . . , u n J+p ), which corresponds to a larger set of cells. In what follows the cells (x j−1 , x j ) with j = 1 − r, . . . , 0 and j = J + 1, . . . , J + p will be referred to as "ghost cells". They are depicted in red in Figure 2.1 (in that example, p = r = 2).
Before describing the numerical boundary conditions we enforce for (2.3), let us state our main and in fact only assumption on the coefficients in (2.3).
Assumption 2.1 (Consistency and stability). The coefficients a −r , . . . , a p in (2.3) satisfy a −r a p = 0 (normalization), and for some integer k ≥ 1, there holds: Provided that the relations (2.4) are satisfied for m = 0 (conservativity) and m = 1 (consistency of order 1) with a > 0, the stability assumption (2.5) implies r ≥ 1, which we assume from now on. Though we view this observation here, as in [Str62], as a necessarycondition for stability, the condition r ≥ 1 is also known to be necessary for convergence by comparing the numerical and continuous dependency domains, see [CFL28]. Let us observe that (2.5) is a necessary and sufficient condition for stability of the iteration process (2.3) on ℓ 2 (Z) in a strong sense, meaning here that the map is a contraction (it has norm 1) as an operator on ℓ 2 (Z). However, (2.5) is not sufficient to yield stability in ℓ ∞ (Z) for (2.3), see [Tho65,Hed66]. Note that through the dependence of the a ℓ with respect to λ = ∆t/∆x, (2.5) is usually intended to be true only under a so-called Courant-Friedrichs-Lewy condition asking for λ to be less than some constant depending on the scheme and the velocity a. (Indeed, the Bernstein inequality for trigonometric polynomials implies λ |a| ≤ max(p, r), see [Str62].) We postpone the extension of our work to ℓ ∞ (Z) stable discretizations to a future work since the methods to be used in that framework will definitely need to be different from the more "Hilbertian" techniques we use here.
Let us observe that Assumption 2.1 does not include any dissipative behavior for (2.3), meaning that we do not assume a bound of the form: for some suitable integer q and positive constant c. In that respect, the framework of Assumption 2.1 is more general than the works [Kre68,Gol77,GKS72,Osh69] and following works that are based on these pioneering results. We thus expect that our approach may be useful to deal with multidimensional problems in which dissipativity is most of the time excluded (or restrictive).
If the interior cells of the grid are labeled, as in (2.3), by j ∈ {1, . . . , J}, the numerical approximation of (2.1) requires, for passing from one time index n to the next, prescribing r numerical boundary conditions on the left of the interval (that is, close to x = 0), and p numerical boundary conditions on the right (that is, close to x = L). In other words, we need to prescribe the value of the approximate solution (u n j ) in the ghost cells located at the boundary of the interior domain. For simplicity, and in order to be consistent with the continuous problem (2.1), we prescribe Dirichlet homogeneous boundary conditions in conjunction with (2.3) on the left of the interval (0, L): u n ℓ = 0 , ℓ = 1 − r, . . . , 0 .
(2.6) On the right of the interval (0, L), there is nothing to be done if p = 0, that is, in the case of an upwind discretization, for in that case, given the vector (u n 1 , . . . , u n J ), the vector (u n+1 1 , . . . , u n+1 J ) is entirely determined by (2.6) and (2.3), so we can iterate the scheme (2.3)-(2.6) starting from some initial data (u 0 1 , . . . , u 0 J ) to any positive time level n. We therefore assume from now on p ≥ 1, which is the interesting case where the numerical discretization of (2.1) necessitates an outflow numerical boundary condition while the continuous problem does not "obviously" provide with one. In this article, we shall prescribe Neumann type numerical boundary conditions (these are called extrapolation numerical boundary conditions in [Gol77]). For ease of writing, we introduce the difference operator in space which acts on vectors (v j ) j=1−r,...,J+p as follows: Higher order difference operators D m − , m ≥ 2, are defined accordingly by iterating D − . Then given a fixed integer k b ∈ N (b stands for "boundary"), we prescribe the following numerical boundary condition in conjunction with (2.3): (2.7) If k b = 0, this corresponds to prescribing homogeneous Dirichlet boundary conditions, while if k b = 1, this corresponds to the standard Neumann numerical boundary condition: The iteration (2.3), (2.6), (2.7) thus proceeds as follows, see Figure 2.2 for an illustration. Given the vector (u n 1 , . . . , u n J ) for some time level n, one first determines the ghost values (u n 1−r , . . . , u n 0 , u n J+1 , . . . , u n J+p ) by (2.6) and (2.7). The new vector (u n+1 1 , . . . , u n+1 J ) is then determined by applying (2.3). It is assumed that J ≥ 1 in order to make the space step ∆x = L/J meaningful and to have at least one cell in the interval (0, L). Let us emphasize that the same procedure (2.7) is applied at each ghost cell close to the outflow boundary x = L. In the terminology of [GT78,GT81], the numerical boundary conditions are of translatory type.
The scheme (2.3), (2.6), (2.7) is initialized with the piecewise constant projection of the initial condition for (2.1), that is, for the interior cells: Our main result is the following convergence estimate for the scheme (2.3), (2.6), (2.7) supplemented with the initial condition (2.8).
The restriction to numerical schemes with two time levels only in (2.3) is necessary since, for instance, the Neumann numerical boundary condition (2.7) is known to yield violent instabilities when used in conjunction with the leap-frog scheme [Tre84]. We postpone the study of outflow numerical boundary conditions for general multistep schemes to a future work.
The integer k b in (2.7) prescribes the approximation order of the numerical outflow boundary condition. In particular, Theorem 2.2 shows that the Neumann numerical boundary condition (k b = 1) gives a local, stable hence convergent way to approximate the exact discrete transparent boundary condition for (2.3) which is nonlocal in time (see [Cou17] for a general derivation of discrete transparent boundary conditions).
There is a "loss" of 1/2 in the somehow expected rate of convergence min(k, k b ) in (2.9). We emphasize once again that Assumption 2.1 does not imply stability of (2.3) in ℓ ∞ (Z) and therefore even on the whole real line, the rate of convergence k in ℓ ∞ cannot be attained in general. With the additional technical complexity of dealing with numerical boundary conditions, we view this loss of 1/2 as a minor disagreement of our method.
Applying estimate (2.9) to the "theoretical" final time T := L/a after which the exact solution to the transport equation (2.1) becomes zero (everything has flowed out of the interval (0, L) through the outflow boundary x = L), we obtain: meaning that the numerical approximation of the solution to (2.1) has become uniformly small on the interval (0, L) at time level N + 1. At later times, stability estimates for the numerical scheme (2.3), (2.6), (2.7), which we shall prove below, assert that the solution (u n j ) to (2.3), (2.6), (2.7) remains "small" (at least on any given finite time interval, J being large).

Preliminary analysis on a half-line
In this section, we show a preliminary result, which is Theorem 3.5 below, that is entirely analogous to our main result, Theorem 2.2 above, except that the space domain is a half-line with an outgoing velocity at the boundary. The result is interesting on its own since a byproduct of our analysis is the verification of the so-called Uniform Kreiss Lopatinskii Condition though we completely bypass the normal mode analysis that is commonly used to verify it [Gol77]. The convergence analysis on the half-line is more classical. In addition to its own interest, which somehow focuses on the Neumann boundary condition without coupling it to Dirichlet boundary condition, Theorem 3.5 below will be a major building block in the proof of Theorem 2.2.

Stability estimates for Neumann numerical outflow boundary conditions
In this paragraph, we prove Theorem 3.1 below that provides us with stability estimates for the Neumann numerical boundary condition on a half-line. Theorem 3.1 is a key tool for proving stability estimates for the scheme (2.3), (2.6), (2.7) on a finite interval, which in turn yields the convergence result of Theorem 2.2. Let us recall that Theorem 3.1 below is already known to hold true thanks to the joint results of [Kre66,Gol77,Kre68] and in a more general setting [Wu95,CG11,Cou13] (see also references therein). It roughly means that the Neumann numerical boundary conditions satisfy the Uniform Kreiss Lopatinskii Condition at an outflow boundary. However we emphasize that Assumption 2.1 does not imply "dissipativity" for (2.3), which already prohibits using the main results in [Kre68,GKS72] and therefore strongly advocates using the energy method, as we shall do, whenever possible. Independently of the subtleties of the assumptions in those works, we give here an elementary proof of the stability result in [Kre66,Gol77] without any Laplace transform nor "GKS" type arguments, and with more easily checkable assumptions. Our integration by parts approach was probably already well-known for three point schemes and the classical Neumann boundary condition D − u = 0, but our main contribution is to show that the energy method can be used in order to deal with numerical schemes with an arbitrarily wide stencil and for any extrapolation order at the boundary.
Before going on, let us fix the space domain that we consider. Since we deal with a constant coefficient linear problem, by translation invariance, there is no loss of generality in considering the half-line (−∞, L). With the same positive constant λ > 0 as in the previous section, we consider J ∈ N * , and some space and time steps ∆x = L/J and ∆t = λ ∆x, assuming that J is large enough to ensure that ∆t ∈ (0, 1]. We also keep the notation t n := n ∆t, n ∈ N and x j := j ∆x, j ∈ Z. The grid and the associated ghost cells are depicted in Figure 3.1. Our result can be stated as follows. Theorem 3.1. Let a > 0, let k ≥ 1 and k b ∈ N, let λ > 0. Then there exists a constant C > 0 such that for all initial data (f j ) j≤J ∈ ℓ 2 and for all boundary source terms (g n J+1 ) n≥0 , . . . , (g n J+p ) n≥0 verifying the growth condition: the solution (u n j ) j≤J+p,n∈N to the scheme: where the coefficients a ℓ satisfy Assumption 2.1 with integer k, satisfies the estimate: for any γ > 0, and for any ∆t and ∆x such that ∆t/∆x = λ and ∆t ∈ (0, 1]. In particular, the numerical boundary conditions in (3.1) satisfy the Uniform Kreiss Lopatinskii Condition.
Before proving Theorem 3.1, let us recall that we always assume the ratio ∆t/∆x to be constant. This will be used several times below and is reminiscent of the scale invariance properties of the underlying continuous problem. Observe also that in (3.2), the larger the integer k b the better the trace estimate on the left hand side behaves. In particular, the numerical boundary conditions in (3.1) involve the values u n J+1−k b , . . . , u n J+p , and we get "for free" in (3.2) not only the control of those terms but also the extra control of u n J+1−r−k b , . . . , u n J−k b (recall r ≥ 1). The fact that (3.2) implies the Uniform Kreiss Lopatinskii Condition is not our main focus here, so instead of recalling many definitions, we rather refer the interested reader to the review [Cou13].
Proof. We shall use Assumption 2.1 in the proof below only for k = 1, that is, we make the "minimal" consistency requirements for the scheme (2.3). Unlike [Kre66,Gol77], the proof of Theorem 3.1 is done by induction with respect to the index k b ∈ N and relies on the energy method. Let us start with the case k b = 0, which corresponds to Dirichlet numerical boundary conditions.
• The case k b = 0. We consider the numerical scheme: A straightforward proof of the stability estimate (3.2) for k b = 0 was achieved in [CG11] (even in some cases of multidimensional systems), see also [GT81] for an earlier general result based on the theory of [GKS72]. We reproduce here the short proof of (3.2) for the scheme (3.3) for the sake of completeness. Assumption 2.1 implies by the Plancherel Theorem that the mapping: is a contraction, in the sense that the operator norm is not larger than 1. Let us now consider the solution (u n j ) j≤J+p,n∈N to (3.3) at some time index n ∈ N. We extend the sequence (u n j ) j≤J+p by 0 for j ≥ J + p + 1 and still denote u n ∈ ℓ 2 (Z) the resulting sequence. Let us then define Observe that, due to the boundary conditions in (3.3), we do not necessarily have v n+1 j = u n+1 j for j = J + 1, . . . , J + p, nor for j = J + p + 1, . . . , J + p + r (extending also (u n+1 j ) j≤J+p by 0 for j ≥ J + p + 1). Using Assumption 2.1, we get: since (3.4) is a contraction. Equivalently, taking the boundary conditions of (3.3) into account, we get: (3.5) We now derive a bound from below for the trace term arising on the left hand side of (3.5). The real numbers v n+1 J+ℓ , ℓ = 1, . . . , p + r, depend linearly on u n J+1−r , . . . , u n J+p . The coefficients in each linear combination are taken among the a ℓ 's. Hence the quantity can be seen as a nonnegative quadratic form in the variables u n J+1−r , . . . , u n J+p . It is also rather easy to see that this quadratic form is positive definite for we have v n+1 J+p+r = a −r u n J+p , and, therefore, if v n+1 J+1 = · · · = v n+1 J+p+r = 0, then we first have u n J+p = 0 and recursively we can also show u n J+p−1 = · · · = u n J+1−r = 0. Hence there exists a fixed constant c 0 > 0, that only depends on the (fixed) coefficients a ℓ in (2.3), such that: Reporting in (3.5) and using that ∆t/∆x = λ is a fixed positive constant, we get for some constant c > 0 We now apply the following discrete Gronwall type lemma (with the positive parameter Γ := γ ∆t), see [CG11] for repeated use of such summation arguments.
Lemma 3.2. Let (G n ) n≥0 be a sequence of nonnegative real numbers such that: Let (U n ) n≥0 , (B n ) n≥0 be two sequences of nonnegative real numbers such that: Then there holds for all Γ > 0: The proof of Lemma 3.2 is straightforward and therefore omitted. We apply Lemma 3.2 to (3.6) and derive the estimate: We emphasize that when dealing with the case k b = 0, we have only used the stability condition (2.5) of Assumption 2.1, and we have never used (2.4) (not even for m = 0). This is consistent with the result of [GT81] which proves that the Dirichlet boundary condition satisfies the Uniform Kreiss Lopatinskii Condition independently of the nature of the boundary (inflow or outflow).
• The induction argument. We now assume that the stability estimate (3.2) is valid up to some index k b ∈ N, and consider the following numerical scheme: (3.8) We use the induction assumption by defining the following sequence: Applying the stability estimate (3.2) for the index k b and omitting one of the two nonnegative terms on the left hand side of (3.2), we have already derived the preliminary estimate: Let us now turn back to the numerical scheme (3.8) to which we are going to apply the so-called energy method. For a given time index n ∈ N, we compute: (3.10) We use the following discrete integration by parts result, whose proof is postponed to Appendix A.
Lemma 3.3. Let a > 0 and let Assumption 2.1 be satisfied for some k ∈ N * . Then there exist a unique quadratic form Q on R p+r and some real coefficients d 1 , . . . , d p+r such that, for any real numbers v j−r , . . . , v j+p , there holds Furthermore, the quadratic form Q satisfies We apply Lemma 3.3 to the right hand side in (3.10). The sum with respect to j ≤ J makes the Q terms a telescopic sum, and we are left with the energy balance: It remains to estimate the "bulk" term on the right hand side of (3.12), which encodes the "ℓ 2 -dissipation" of the discretization (2.3) of the transport equation. More precisely, if we start from a sequence v ∈ ℓ 2 (Z) and consider its image by the contraction (3.4), we can rewrite Equation (3.10) for v and use the decomposition given in Lemma 3.3 to derive the inequality: Let us now consider the sequence (u n j ) j≤J+p ∈ ℓ 2 for the given integer n ∈ N. We extend the latter sequence as an ℓ 2 sequence on the whole set of integers Z by symmetry with respect to the index J + p: (3.14) For ℓ = 1, . . . , p + r and j = J + 1, . . . , J + 2 (p + r) − ℓ − 1, v j+ℓ−r and v j−r on the right hand side of (3.14) are taken among the real numbers u n J+1−r , . . . , u n J+p . We therefore obtain the upper bound: Using the latter bound in (3.14), we get the estimate We now apply the summation argument of Lemma 3.2 to the latter inequality and derive the estimate: We then combine the latter inequality with the preliminary estimate (3.9), which yields: At this stage, we have almost proved that (3.2) holds up to the index k b + 1. Indeed, if we combine the trace estimates provided by (3.9) and (3.15), we get a full control of the trace of (u n ), that is of (u n J+ℓ ) n≥0 at ℓ = −r − k b , . . . , p: Combining with (3.15), we have completed the proof of (3.2) for the index k b + 1, which also completes the proof of Theorem 3.1.

Convergence estimates for Neumann outflow numerical boundary conditions
In the previous paragraph, we have proved the stability estimate (3.2) in order to highlight the fact that our method automatically yields the verification of the Uniform Kreiss Lopatinskii Condition. However, the exponential weights arising in (3.2) and the fact that no "interior" source term is considered in (3.1) make this estimate hardly applicable as such in view of the convergence analysis below. We therefore state a slightly weakened but more practical version of Theorem 3.1 which will help us proving Theorem 3.5 below.
Proposition 3.4. Let a > 0, let k ≥ 1 and k b ∈ N, let λ > 0. Then there exists a constant C > 0 such that for all initial data (f j ) j≤J ∈ ℓ 2 , for any N ∈ N * , for all boundary source terms (g n J+1 ) 0≤n≤N −1 , . . . , (g n J+p ) 0≤n≤N −1 , and all interior source terms (F n j ) j≤J,1≤n≤N with (F n j ) j≤J ∈ ℓ 2 for all n = 1, . . . , N , the solution (u n j ) j≤J,0≤n≤N to the scheme: Proof. By linearity of (3.17), it is sufficient to examine separately the cases F ≡ 0 (no interior source term), and f ≡ 0, g ≡ 0 (interior forcing only). In the case F ≡ 0, the estimate (3.18) is directly obtained by: • first extending the boundary source terms (g n J+ℓ ) ℓ=1,...,p by 0 for n > N − 1, which does not affect the solution to (3.17) for j ≤ J at time steps earlier than N , • then passing to the limit γ → 0 in (3.2) (and forgetting about the nonnegative trace estimate on the left hand side of (3.2)).
It therefore remains to examine the case f ≡ 0, g ≡ 0, which is done by using the Duhamel formula. Namely, Theorem 3.1 shows that the solution (v n j ) j≤J+p,n∈N to the numerical scheme with homogeneous boundary conditions: satisfies the uniform in time bound Writing v n under the form S n f , this means that the operator S , which is the generator of the discrete semigroup (S n ) n∈N , is power bounded on ℓ 2 (−∞, J). At this stage, it remains to observe that the solution to (3.17) in the case f ≡ 0, g ≡ 0, can be written with the Duhamel formula: Since S is power bounded on ℓ 2 , we end up with: This completes the proof of (3.18).
We are now ready to prove our convergence result for the Neumann boundary condition at an outflow boundary.
Theorem 3.5. Let a > 0, let k ∈ N * and k b ∈ N, let λ > 0. Then there exists C > 0 such that for any T > 0 and J ∈ N * , for any L ≥ 1, and for any u 0 ∈ H k+1 ((−∞, L)), the solution (u n j ) j≤J,n∈N to where the coefficients a ℓ satisfy Assumption 2.1 with integer k, satisfies, with k 0 := min(k, k b ): (3.20) where ∆x = L/J and ∆t = λ ∆x (assumed to be less than 1). In particular, the following global ℓ ∞ convergence estimate is satisfied: Proof. To pass from (3.20) to (3.21) is very crude, and simply relies on the inequality: for any real sequence (b j ) j≤J . We therefore focus on the derivation of the ℓ ∞ n (ℓ 2 j ) estimate (3.20) which unsurprisingly relies on the stability estimate (3.18) provided by Proposition 3.4. In the proof below, we assume k b ≤ k, so that the limiting order of convergence arises from the numerical boundary conditions in (3.19) and not from the discretization of the transport equation. The proof in the case k b > k is quite similar and we leave the corresponding modifications to the interested reader. Since the validity of Assumption 2.1 for some integer k ≥ 1 implies the validity of the relations (2.4) for the restricted subset of indices m = 0, . . . , k b , we can even assume without loss of generality k b = k.
We now denote by U 0 the extension of u 0 as a function in H k+1 (R) by the linear continuous reflexion operator of [DL85], and define, for any j ∈ Z, to be the cell average of the exact solution ((t, x) −→ U 0 (x − a t)) to the transport equation over R.
For the rest of this proof, we will denote the integer part of T /∆t by N T . Then for all n = 0, . . . , N T and j ≤ J + p, we define the error ε n j := u n j − w n j . By definition of the sequences (u n j ) and (w n j ), the error (ε n j ) j≤J,n=0,...,N T is a solution to: Lemma 3.6. Let m, p ∈ N. There exists a constant C > 0 such that for any ∆x > 0 and for any v ∈ H m+1 ((−∞, L + p ∆x)), defining v j := 1 ∆x H 1 ((−∞,L+p ∆x)) , ℓ = 1, . . . , p .
Proof of Lemma 3.6. For ℓ = 1, . . . , p, there holds where we have used Taylor's formula and cancellation properties of the binomial coefficients. Applying the Cauchy-Schwarz inequality to each double integral in the latter expression, we get The proof follows by using the imbedding of H 1 in L ∞ in one space dimension.
We now apply Lemma 3.6 with m = k b to evaluate the boundary errors in (3.23). We get: so the remaining task is to evaluate the size of the interior consistency error (e n j ) j≤J . By its definition (3.22), the value of w n j in each interior cell is the average of the extension of the exact solution, that is the cell average of U 0 (· − a t n ). We thus have, for any j ≤ J and n ≥ 1: where we defined U n−1 ∈ H k+1 (R) as U 0 (· − a t n−1 ). We also define the consistency error on the whole domain: so that we have e n j = E n j for all j ≤ J. We now use Fourier analysis and derive the estimate: We now use Assumption 2.1 (recall k b = k), to derive the bound: uniformly with respect to ∆x and ξ ∈ R. Recalling that the ratio ∆t/∆x is kept fixed, we end up with the consistency estimate: where we have used the fact that U n−1 has been constructed as an extension of u 0 (· − a t n−1 ) boundedly with respect to the H k b +1 norm. Going back to (3.26) and using the estimate (3.27), we finally derive sup n=0,...,N T j≤J which completes the proof of (3.20).

Convergence estimates for Dirichlet numerical boundary conditions
This short paragraph is devoted to the complementary boundary value problem on a half-line that can be extracted from (2.3), (2.6), (2.7), and that focuses on the (homogeneous) Dirichlet boundary condition (2.6) at the inflow boundary x = 0. Hence the half-line is now R + . To be consistent with the notations introduced in Section 2, we let x j := j ∆x, j ∈ Z. The ghost cells correspond to the intervals (x −r , x 1−r ), . . . , (x −1 , x 0 ), see Figure 3.2. The time step ∆t is linked to the space step by keeping ∆t/∆x = λ, λ being the same fixed positive constant as in Section 2. We prove the following result, that is the analogue of Theorem 3.5.
Proposition 3.7. Let a > 0, let k ∈ N * , let λ > 0. There exists C > 0 such that for any T > 0 and J ∈ N * , for any L ≥ 1 and for any u 0 ∈ H k+1 ((0, +∞)) satisfying the flatness conditions: the solution (u n j ) j≥1,n∈N to the scheme: where the coefficients a ℓ satisfy Assumption 2.1 with integer k, satisfies where ∆x = L/J and ∆t = λ ∆x. It is understood in (3.29) that the initial condition u 0 has been extended by 0 to R − . In particular, the following global ℓ ∞ convergence estimate is satisfied: The proof follows more or less the exact same arguments as the proof of Theorem 3.5. One first defines the cell average of the exact solution ((t, x) → u 0 (x − a t)) to the transport equation in the half space R + with homogeneous Dirichlet boundary condition at x = 0. The flatness assumption on u 0 ensures that extending u 0 by 0 on R − yields a function in H k+1 (R). The error between the solution to (3.28) and the cell average of the exact solution will satisfy a system of the form (3.28) with zero initial condition and homogeneous Dirichlet boundary condition but with a nonzero interior forcing term (the so-called consistency error). Estimating this error is done by using a similar (and even easier) extension procedure as in the proof of Theorem 3.5 and Fourier analysis. We leave the details to the interested reader. The stability estimate for (3.28) is analogous to that of Proposition 3.4 and follows from Theorem 3.1 (in the case k b = 0) since we have already observed that dealing with this case does not necessitate the consistency conditions (2.4) but only the stability assumption (2.5).

Proof of Theorem 2.2 4.1 Stability estimates on a finite interval
We now turn to the study of the numerical scheme (2.3), (2.6), (2.7), which is an iteration in a finite dimensional space and therefore really corresponds to a numerical scheme that can be implemented in practice. Let us recall that the space step ∆x is given by ∆x = L/J with J a positive integer, the spatial grid is defined by means of the points x j := j ∆x, j ∈ Z, and the time step is given by ∆t = λ ∆x with λ a fixed positive constant. We first prove a stability estimate for (2.3), (2.6), (2.7), which will have various consequences.
Proposition 4.1. Let a > 0, let k ≥ 1 and k b ∈ N, let λ > 0. Then there exists a constant C 0 > 0 such that for all initial data (f 1 , . . . , f J ), the solution (u n j ) 1≤j≤J,n∈N to the scheme: where the coefficients a ℓ satisfy Assumption 2.1 with integer k, satisfies The estimate (4.2) is compatible with the limit L → +∞ (and ∆x, ∆t fixed) where we formally recover either the stability estimate (3.18) or its analogue for Dirichlet boundary conditions on a half-line, at least in the case of homogeneous boundary conditions and no interior source term. More important, the estimate (4.2) also has a consequence in terms of the spectral radius of the perturbed Toeplitz matrix associated with (4.1). In the absence of the outflow boundary at x = L, the numerical scheme (3.28) with homogeneous Dirichlet boundary conditions at the inflow boundary x = 0 corresponds to iterating the Toeplitz operator represented by the semi-infinite matrix:  a 0 · · · a p 0 · · · · · · · · · . . . . . . . . . . . . a −r · · · a 0 · · · a p . . .
Acting on ℓ 2 (N), the spectrum of this operator is known to contain at least the closed curve { a ℓ exp(i ℓ θ), θ ∈ [0, 2 π]}, see [TE05, Chapter 7] and references therein for a precise statement. In particular, Assumption 2.1 and [TE05, Theorem 7.1] show that the spectrum of this operator lies inside the closed unit disk {z ∈ C , |z| ≤ 1} and accumulates at the point 1. Truncating the above semi-infinite matrix to make it of size J corresponds to prescribing homogeneous Dirichlet boundary conditions at the outflow boundary. In that case, the resulting operator is of the form Π T Π with Π an orthogonal projection and T a contraction. In particular, prescribing homogeneous Dirichlet boundary conditions at the outflow boundary leaves the spectrum (consisting now of finitely many eigenvalues) within the closed unit disk. Prescribing Neumann numerical boundary conditions as in (4.1) corresponds to first truncating the semi-infinite matrix to make it of size J, and then making "large" perturbations of the coefficients in the k b last columns of the last p rows of the truncated matrix. For instance, if r = p = 1 and k b = 1, resp. k b = 2, the corresponding matrix reads:  It is not obvious at first sight that the spectrum of the matrix A J ∈ M J (R) associated with the iteration (4.1) will not deviate from the closed unit disk at a distance O(1), or even O(J −α ) for some α ∈ (0, 1), giving rise to instabilities for (4.1). However, the estimate (4.2) has an important consequence for the matrix A J since it shows in particular that its ℓ 2 induced norm verifies ∀ n ∈ N , A n J 2 ≤ C e C n ∆t/L , so taking the 1/n-th power and passing to the limit (recalling ∆t/L = λ/J, λ > 0 constant), the spectral radius of A J satisfies ρ(A J ) ≤ exp(C/J) .
Hence the spectrum of A J does not deviate from the closed unit disk at a distance larger than O(1/J).
(The bound in (4.2) is even more precise since it even predicts the behavior of the ε-pseudospectrum of A J , see [TE05], prohibiting in particular Jordan blocks associated with eigenvalues of largest modulus.) Proof of Proposition 4.1. The derivation of (4.2) follows from the finite speed of propagation of the numerical scheme (2.3). (Observe that our argument below does not extend to implicit discretizations of the transport equation.) More precisely, let us assume for simplicity that J is even. Let Then for n ≤ N 0 , we can write the solution to (4.1) as the superposition of solutions to two initial boundary value problems of the form (3.19) and (3.28). Indeed, if the initial condition (f j = u in j ) j≤J in (3.19) vanishes for j ≤ J/2, then the solution to (3.19) satisfies the homogeneous Dirichlet boundary condition: ∀ n ≤ N 0 , ∀ ℓ = 1 − r, . . . , 0 , u n ℓ = 0 , because the support of u in is shifted of p cells to the left at each time iteration. In particular, the restriction of the solution to the half-line problem (3.19) with some initial condition u in vanishing for j ≤ J/2 satisfies the numerical scheme (4.1) up to n = N 0 . Similarly, if the initial condition (u in j ) j≥1 in (3.28) vanishes for j ≥ J/2 + 1, then the solution to the half line problem (3.28) satisfies: because the support of u in is shifted of r cells to the right at each time iteration, and the definition of N 0 so the solution to (3.28) vanishes at all points that are used in the calculation of the finite differences (D k b − u n ) J+ℓ . Truncating the initial condition for (4.1) to the left and to the right of the medium cell J/2, we can, thanks to the linearity of the scheme, combine the stability estimate (3.18) and the analogous one for Dirichlet boundary conditions to show that the solution to (4.1) satisfies: with a constant C that is independent of L and J. It then remains to iterate the latter estimate on each interval of integers {0, . . . , N 0 }, {N 0 , . . . , 2 N 0 } and so on to derive (4.2).

Convergence
We now prove Theorem 2.2, whose proof combines the convergence estimates given by Theorem 3.5 and Proposition 3.7, and the stability estimate of Proposition 4.1. Let us recall that we consider for the continuous problem transport an initial condition u 0 ∈ H k+1 ((0, L)) that satisfies the flatness conditions: We therefore extend u 0 by zero on R − , considering now u 0 as an element of H k+1 ((−∞, L)) that vanishes on R − . Let us consider a C ∞ function χ on R that satisfies: and then decompose u 0 as follows: .
The function v 0 belongs to H k+1 ((−∞, L)) and vanishes for x ≤ L/3. Conversely, the function w 0 belongs to H k+1 ((−∞, L)) and vanishes for x ≥ 2 L/3. Moreover, using L ≥ 1, there holds uniformly in L: where C depends only on the fixed function χ (and on k). Let N ∈ N denote the largest integer such that We assume that J is large enough so that N ≥ 1. Applying the same finite speed of propagation argument as in the proof of Proposition 4.1, we can write ∀ n = 0, . . . , N , ∀ j = 1, . . . , J , u n j = v n j + w n j , where (v n j ) j≤J,n=0,...,N satisfies an iteration of the form (3.19), namely: and (w n j ) j≥1,n=0,...,N satisfies an iteration of the form (3.28), namely:  The initial conditions in (4.3) and (4.4) are defined as follows. Following (2.8), we define the discretized initial condition associated with v 0 for (4.3) by letting 1 : To define the initial condition for (4.4), let us recall that 1 − χ(·/L) vanishes on the interval (2 L/3, L) so we can extend w 0 by zero to (L, +∞) and view w 0 as an element of H k+1 ((0, ∞)) that satisfies the flatness conditions of Proposition 3.7. We then define the initial condition for (4.4) by letting: We can apply Theorem 3.5 to (4.3) and get the convergence estimate: with k 0 := min(k, k b ) and where the estimate is uniform with respect to L. Applying Proposition 3.7 to (4.4), we get the other convergence estimate: which we combine with (4.5) to derive: Here we use a specific notation C 1 for the constant in (4.6) in order to emphasize its role in the concluding argument below (in the same way, we use a specific notation C 0 for the constant in the stability estimate (4.2)). It remains more or less to iterate in time (4.6). Namely, at the time level N , (4.6) shows that we can write: To iterate the process, we define (u N,n j ) 1≤j≤J,0≤nN by Note that u 0 (· − a t N ) satisfies the flatness conditions of Theorem 2.2, and that: Thus, on the one side, (4.6) applies and we have On the other side, the sequence (δ n j := u N +n j − u N,n j ) 1≤j≤J,0≤nN satisfies: so that, thanks to the stability estimate (4.2) of Proposition 4.1, for any n ≤ N , Finally, thanks to the triangle inequality, Applying iteratively the same argument, we end up with where we have assumed C 0 ≥ 2 without loss of generality. It remains to choose m := E(N ∆t/T ) + 1, which by definition of N is uniformly bounded with respect to J (N scales like c J with c > 0 constant, and ∆t scales like c ′ /J), and we end up with: The convergence estimate (2.9) of Theorem 2.2 follows by a direct lower bound for the norm on the left hand side.

Numerical experiments
In this section, we illustrate with numerical examples the results proved in the paper. We consider the second order in time and space Lax-Wendroff scheme, that writes which means, in the present formalism, that r = p = 1 and a −1 = (a 2 λ 2 + a λ)/2, a 0 = 1 − a 2 λ 2 , a 1 = (a 2 λ 2 − a λ)/2 and k = 2. This scheme is stable (over Z) and second order accurate provided that λ |a| ≤ 1, and we choose a = 1 and λ = 0.7 in the numerical simulations reported below. We compare the results with the first order Neumann outflow condition (k b = 1) and the second order one (k b = 2), for various initial conditions. Namely, we consider L = 1 for the length of the space interval and we choose as initial data, u 0,1 (x) = ((x − 1/2) + ) 3 , x ∈ (0, 1) , u 0,2 (x) = ((x − 1/2) + ) 2.6 , x ∈ (0, 1) , u 0,3 (x) = ((x − 1/2) + ) 2.5 , x ∈ (0, 1) , which satisfy u 0,1 ∈ H 3 , u 0,2 ∈ H 3 and u 0,3 ∈ H 3−ǫ for any ǫ > 0. The final time is T = 0.5 (after T the exact solution vanishes). The convergence result provided in this paper applies for the initial data u 0,1 and u 0,2 , and forecasts that the l ∞ error is bounded by C ∆x 1.5 if k b = 2 and C ∆x 0.5 if k b = 1. Figure 5.1 represents the initial condition u 0,2 on a grid with 40 cells on (0, 1). On figure 5.2 and 5.3, we plot the numerical solutions (at different times) obtained with k b = 0 (homogeneous Dirichlet outflow condition), k b = 1 (homogeneous Neumann outflow condition) and k b = 2 (homogeneous "second order" Neumann outflow condition). As expected, the Dirichlet condition shows a larger boundary layer, and, especially at time 0.2625, the solution with k b = 2 is much nearer to the exact solution than the others.
Let us now analyze more precisely the error of the schemes. In the following, the computed "error" components actually are not the u n j − x j x j−1 u 0 (y − at n ) dy/∆x: for simplicity, they are replaced with the local values u n j − u 0 ((x j−1 + x j )/2 − at n ). If u 0 ∈ H 2 , |u 0 ((x j−1 + x j )/2 − at n ) − x j x j−1 u 0 (y − at n ) dy/∆x| can be bounded by a constant times ∆x 2 , thus the present formula will not disrupt our evaluation of the order of the scheme, which cannot be greater than 2. By observing Table 1, we see that: • with k b = 2 the order of the scheme indeed seems to be 2, • while with k b = 1, it seems to be 1.
In Table 2 with a less smooth (but still in H 3 ) initial datum, • with k b = 2 the order of the scheme seems to belong to (1.7, 1.75), • and with k b = 1, it seems to be 1.
At last, Table 3 suggests that: • with k b = 2 the order of the scheme be in (1.65, 1.7), • and with k b = 1, the order still be 1.
Clearly the convergence result we have obtained here seem to be non-optimal and there is hope, for numerical schemes such as (5.1) that fit into the framework of [BC17], to fill this discrepancy by using numerical boundary layer expansions. However, one clearly observes that the convergence order does depend not only on the smoothness of the solution but also on the order of the homogeneous Neumann condition for the outflow boundary.
In order to illustrate the discussion about the stability of the scheme (4.1) and the properties of the matrix A J ∈ M J (R) (Section 4.1), we report in Table 4 the spectral radius and the l 2 induced norm of this matrix (both are computed approximately thanks to the functions eig and norm of the library numpy.linalg in the Python language), with respect to the boundary condition and the number of cells. We observe that in all cases the spectral radius is smaller than 1, but that with the second order Neumann condition, the norm of the matrix is larger than 1. This means that stability estimates cannot be obtained by showing that the ℓ 2 -norm does not increase, hence the need for some well-designed analytical tool (here we have used an induction argument with respect to k b combined with the finite speed of propagation).     Table 4: Spectral radii and l 2 induced norms of the linear operator associated with the scheme.

A A discrete integration by parts lemma
In this appendix, we prove the following result, of which Lemma 3.3 is an immediate corollary as explained below.
Lemma A.1. Let S ∈ M m (R), m ≥ 2, be a real symmetric matrix satisfying m i,j=1 S ij = 0 .
Then there exists a unique real symmetric matrix S of size m − 1, and some unique real numbers d 1 , . . . , d m−1 , such that: Proof. The proof is completely elementary. The vector space of real symmetric matrices of size m satisfying the condition m i,j=1 S ij = 0 , has dimension m (m + 1)/2 − 1 = (m − 1) (m + 2)/2. By standard linear algebra, it is therefore sufficient to prove that if a real symmetric matrix S of size m − 1 and some real numbers d 1 , . . . , d m−1 satisfy: Let us now explain how Lemma A.1 gives the result claimed in Lemma 3.3. On R p+r+1 , with vectors written under the form (v −r , . . . , v p ), we consider the quadratic form Because of (2.4) for m = 0, the vector (1, . . . , 1) belongs to the isotropic cone of this quadratic form. Hence we can decompose the real symmetric matrix S associated with this quadratic form by using Lemma A.1. We get a decomposition of the form for some suitable real quadratic form Q on R p+r . It only remains to make the invertible change of variables ( in the argument of Q, which modifies this quadratic form into some Q, and we obtain the result of Lemma 3.3 as announced. The value of Q on the r-th vector of the canonical basis of R p+r is obtained by considering the specific sequence: ∀ j ∈ Z , v j := j , and by identifying the dominant term in j.