QUADRATIC STABILITY OF FLUX LIMITERS

. We propose a novel approach to study the quadratic stability of 2D flux limiters for non expansive transport equations. The theory is developed for the constant coefficient case on a cartesian grid. The convergence of the fully discrete nonlinear scheme is established in 2D with a rate not less than 𝑂 (Δ 𝑥 12 ) in quadratic norm. It is a way to bypass the Goodman–Leveque obstruction Theorem. A new nonlinear scheme with corner correction is proposed. The scheme is formally second-order accurate away from characteristics points, satisfies the maximum principle and is proved to be convergent in quadratic norm. It is tested on simple numerical problems.


Introduction
A famous Godunov Theorem [20] states that a linear scheme for advection which satisfies the maximum principle is only first-order accurate, except in trivial cases. Subsequently, Harten proposed in [25] a constructive notion for Total Variation Diminishing (TVD) numerical schemes which allows the construction of formally second-order accurate nonlinear numerical schemes for advection in space dimension = 1 (1D). To be precise, we follow the literature [21,25] for which formal second-order and high-order accuracy means accuracy in the sense of local Taylor expansion for exact smooth solutions, away from extrema because extrema yield obstruction even for non linear schemes ( [31], Lem. 2.3). The theory of nonlinear TVD schemes is nowadays well established, see [16,19,35] and references therein.
However, in contrast with the success of the TVD theory in 1D, Goodman and Leveque proved in [21] the obstruction Theorem summarized as follows: A TVD numerical scheme with compact stencil for the discretization of conservation laws in space dimension ≥ 2 is at most first-order accurate, except in trivial cases. The proof [21] is by contradiction: one starts from a scheme which satisfies the TVD property and which has a compact stencil in space dimension = 2 (2D); then one proves that this scheme is 1 contracting for one dimensional profiles; it turns into the fact that the scheme is just first-order accurate by another important result [26]. It is a non constructive proof which does not propose any way to go beyond this limitation (nevertheless we immediately remark that 2 stability of linear schemes for 1D advection can already be achieved at any order [11], contrary to linear stability in 1 norm).
Since then, the development of numerical solvers has pursued its route [2,19,35] in space dimension ≥ 2, but without the compactness offered by TVD inequalities. This regrettable situation where there is a divorce between the development of non linear high order numerical methods and the theory of numerical convergence is a direct corollary of the Goodman-Leveque Theorem. Quoting [21]: While it is not logically necessary for a scheme to be TVD in order to be total variation stable, or to converge, we know of no method for computing weak solutions that is not TVD and yet can be shown to converge. Optimum positive schemes for transport by Roe and Sidilkover [32] are linear and first order so cannot constitute a solution. The Barth-Jespersen technique [1] on unstructured meshes guarantees the maximum principle, but the scheme is ultimately first-order accurate. One finds in Barth-Ohleberger [2] a very nice and comprehensive review on the topic of limiter techniques and on the consequences of the Goodman-Leveque obstruction theorem on the development of numerical methods, where the emphasis moves on the preservation of the maximum principle which is a weaker notion than TVD estimates (actually the Godunov Theorem [20] is about the maximum principle). Some technical issues are solved in [7] where a family of formally maximum-preserving and second-order schemes for conservation laws is studied. Provided the schemes in the family satisfy an additional hypothesis (15)(16) which ensures the control of the correction term in the numerical flux, a convergence estimate (∆ 1 4 ) is obtained in 1 norm. Other progresses have been made by systematically using bound-preserving or invariant-domain-preserving numerical methods. Bound-preserving numerical linear and nonlinear methods are extensively studied and developed in the works of Shu [6,30,33] and therein. Refer to [36] for the principles of bound-preserving Finite Volume methods, but even in the latest works, the convergence of fully discrete in time numerical schemes is not established. Domain invariant techniques and positivity-preserving techniques are also studied and developed in the contributions of Guermond-Popov et al. [22,23], still without proof of convergence for fully discrete schemes (see e.g. [22], Thm. 4.8 for quite a technical statement on the stability of a family of bound preserving method for transport). In summary the Goodman-Leveque obstruction Theorem is still in the background, even in the the latest development.
The purpose of this work is to show that new quadratic stability estimates for nonlinear flux limiter techniques offer a strategy to go beyond the Goodman-Leveque obstruction Theorem, at least for transport equations on a cartesian grid. With quadratic stability, a control of nonlinear fully discrete numerical schemes is achieved in 2D for a large class of multidimensional flux limiters, allowing formally second-order accurate (away from characteristics points) and maximum-principle-satisfying numerical schemes for which one can proves numerical convergence by means of quantitative estimates. With respect to [7], the improvement is that the rate of convergence is (∆ 1 2 ) in quadratic norm, so it is a stronger result than (∆ 1 4 ) in 1 norm. Moreover we do need the hypothesis (15-16) of [7] because the new quadratic stability estimate provides control of these terms.
The problems embraced in our analysis must have some non expansive properties in quadratic norm, because the goal is to reproduce these properties at the discrete level. Such a general non expansive equation in general dimension writes where the given velocity field is a ∈ 1,∞ (R ) and the source is ∈ ∞ (R ). For convection dominated flows in applied sciences, the calculation of accurate numerical solutions to this equation is of fundamental importance [18,19,29,35]. If the velocity field is divergent free, that is ∇ · a = 0, then the equation is non expansive in all (R ) norms and in particular in quadratic norm (this is 2 (R ) norm). If the divergence is bounded ∇ · a ∈ ∞ (R ) the expansion is bounded by classical estimates. An equation similar to (1) can be obtained via the linearization of a truly nonlinear equation tot + ∇ · (b(x) ( tot )) = 0 where the final equation for the perturbation defined for exemple by tot = ref + + ( 2 ) can be written as (x, ) + ∇ · (a(x, ) (x, )) = 0: the velocity field a(x, ) is a function of space and time; non expansive estimates can be obtained depending on the space-time regularity of the velocity field. Provided the regularity of the velocity field is sufficient, it is possible to approach weak solutions with smooth solutions (non smooth velocity fields low regularity such as in [5] are not considered in this work). For the simplicity of the theoretical developments, we will assume that the initial data has compact support with three bounded derivatives To concentrate on the main features of our approach and to minimize the amount of notations in 2D, the equation is simplified further. The model problem is 2D advection (that is transport at constant velocity) discretized on a cartesian grid Additionally the advection velocity field is normalized ≥ 0, ≥ 0 and + = 1.
One dimensional configurations aligned with the discretization grid correspond to the limit cases ( , ) = (1, 0) or ( , ) = (0, 1). The equation (3) is already relevant for discretization of kinetic equations, for which = (x, v, ) ≥ 0 represents the density of particles in function of the space variable, the time variable and the velocity v ∈ R variable. The correspondance is v = ( , ) where ≥ 0 is a rescaling factor to account for the normalization + = 1. Considering the fierce activity nowadays in the field of numerical plasma physics [3,8,14,17,28], it is possible that the original methods for the numerical analysis of the advection equation (3) and perhaps the original scheme constructed at the end of this work for the purposes of numerical illustrations can already find immediate applications in this field. We think also of the applications of flux limiters with quadratic stability for the remap stage of Lagrange+remap schemes [9] applied to the numerical discretization of compressible Euler equations (a different technique is [27]). Radiation transport equations [22] could also benefit from the techniques of this work. Extension to non constant coefficients are evoked in the last Section.
The principle of the new quadratic estimates is exposed in 1D in Section 2 and the application to the convergence of the fully discrete schemes is explained in Section 3. The extension to the 2D case is performed in Section 4 where the quadratic stability and the convergence are proved in Theorem 12 which is the main theoretical contribution of this work. A new scheme is constructed in Section 5 with a notion of corner interaction which is natural in the context of this work. This scheme is formally second-order accurate except at characteristics points and satisfies the maximal principle. It converges in quadratic norm by virtue of the theoretical estimates, at a rate not less than (∆ 1 2 ). Numerical results obtained with the new schemes are shown in Section 6, in particular it is clear that the theoretical estimate of convergence of Section 4 is largely suboptimal. It is also clear that the new non linear scheme is more accurate than the linear Lax-Wendroff for truly weak solutions with BV regularity. Some natural extensions are evoked in last Section 7.
Warning. All material in this work is restricted to cartesian meshes with constant space step ∆ = ∆ . It is possible to consider cartesian meshes with different space steps ∆ ̸ = ∆ . Extension to general meshes is an open problem, in particular because the theory of basic non linear schemes for transport on general meshes is trickier as explained in [24], see also [1,2,7,19,29,30,35,36].

Quadratic estimates in 1D
Quadratic stability of flux limiters is easier to explain 1D than in 2D and the result of the analysis can be compared with the classical 1D theory of flux limiters [18,19,25,35]. A generic scheme for the advection equation with In these notations, stand for which is the numerical solution at time step = ∆ and stands for +1 which is the numerical solution at time step +1 = ( + 1)∆ . An explicit formulation is )︁ where the CFL number is = ∆ ∆ ∈ (0, 1]. As usual the CFL number is in a restricted range to insure the numerical stability. The scalar + 1 2 is the flux limiter between two consecutive cells (or mesh points). It is well known [29,35] that if + 1 2 ≡ 0 for all then one obtains the upwind scheme, if + 1 2 ≡ 1 for all then one obtains the Lax-Wendroff scheme. The theory of TVD flux limitation has been developed in [25,34] and is exposed in [18,19,29,35]. For example the minmod limiter is defined by The minmod function is defined by minmod( , ) = 0 for ≤ 0 and minmod( , ) = sign( ) min(| |, | |) for ≥ 0. Contrary to the classical references [18,19,29,35] where the analysis is 1 -based or TVD-based, our analysis is based on the quadratic discrete Lebesgue space In the next preliminary result, the 2 norm of the numerical solution at the end of the time step is compared with the 2 norm of the numerical solution at the beginning of the time step. We note ∆ = The parameters of the quadratic form are noted = (︁ Proof. The scheme (5), (6) is written explicitly as where the upwind scheme is up = (1 − ) + −1 and the difference is ∆ + 1 2 = +1 − for all . Denoting = | | 2 − | up | 2 , one can write with these expressions, one can calculate ∑︀ The differences in (11), (12) become telescopic by summation, so they vanish and the claim is obtained.
The upwind scheme corresponds to + 1 2 ≡ 0: in this case the quadratic form 0 from (10) is non negative for all ∆ provided 0 < ≤ 1. One recovers without surprise the quadratic stability of the upwind scheme under CFL condition (8). The Lax-Wendroff scheme corresponds to + 1 2 ≡ 1. One obtains the quadratic form 1 where One recovers, also without surprise, the quadratic stability of the Lax-Wendroff scheme under CFL condition (8). The natural question is to determine if there exists a general condition on the limiters = (︁ A trivial answer is interpolation of + 1 2 between 0 and 1. If the interpolation is uniform with respect to , that is + 1 2 = ∈ [0, 1] for all ∈ Z, then the result is evident. The key fact below is that the interpolation can be taken for + 1 2 independently to + 1 2 for all ̸ = . The main technical difficulty concerns the quadratic form This quadratic form is a part of (∆) since one has the formula Proposition 2. Assume the CFL condition (8) and assume 0 ≤ Proof. The condition ∆ ∈ 2 guarantees the convergence of the sum in (10), so (∆) is finite. The proof uses a rearrangement of (10), which is valid for ∆ ∈ 2 . Since (∆) is linear with respect to , we note An expansion shows that So (∆) = (∆). Therefore (∆) = (∆) + (∆), that is The conditions of the claim insure (∆) ≥ 0.
Theorem 3. Assume the conditions of Proposition 2. Then the scheme is stable in quadratic norm: Proof. The identity (14) yields the claim.
Proof. Taking into account the strict CFL condition 0 < < 1, one gets two inequalities from (14) and (15). The first one is Remark 1 (Control of discrete gradients). The bound (16) is similar to what is called weak BV estimates in [16]. This estimate shows a quadratic control of the second discrete gradient ∆ + 1 2 − ∆ − 1 2 . It also shows a quadratic control of the first gradient, multiplied by 1 − + 1 2 .

Convergence in 1D
The convenient space for the numerical analysis is the 2 space equipped with a quadratic norm weighted with the mesh size The associated scalar product is denoted ( , ) Δ = ∆ ∑︀ ∈Z for , ∈ Δ . The strategy of the proof of convergence relies on two ideas. The first idea corresponds to the semi-discrete limit regime which correspond to ( ) terms in (∆): involved manipulations are similar to classical ones in Finite Volume techniques [16,33]. The second idea is fully discrete and is necessary to treat the term which come from the discretization in time: it relies on sharp estimates which are new with respect to the literature.

The semi-discrete scheme
The semi-discrete scheme is continuous in time. It is obtained by letting ∆ → 0 or → 0 in (5) and (6). It writes We remind the reader that, for simplicity, we consider a smooth initial data (2). Then the solution is smooth as well. The interpolation of the smooth function ∈ 3 (R × R + ) (with compact support in space) is denoted as ( ) = ( ∆ , ). One defines the truncation error ( ) One has the decomposition from which it is easy to verify that The (∆ 2 ) is uniform with respect to the index and the time . Since is compact in space, then only a finite number of ( ) can be non zero at a given time . The number of non zero terms is (∆ −1 ) uniformly with respect to the time . Let us define the numerical error ( ) = ( ) − ( ) which vanishes at initial time (0) = 0. The error satisfies Proof. Multiply (20) by ( ) and sum over all . One obtains where the residual can be estimated with (19). One gets One has ⃒ ⃒ ⃒ and only a (∆ −1 ) terms are non zero. One gets for small ∆ The Gronwall Lemma yields the claim.

The fully discrete scheme
One combines the analysis of Section 3.1 and the inequality (16). The interpolation of the exact solution on the space-time grid is = ( ∆ , ∆ ). The truncation error is It is decomposed in two parts = + . The first part is the truncation error of the Lax-Wendroff scheme . The second part comes from the flux limitation The error = − vanishes at initial time (that is 0 = 0 for all ) and satisfies By comparison with (21), it is convenient to rewrite it as where the quadratic form is evaluated with respect to ∆ , = (︁ ∆ , (24) does not pose any problem because by definition with the notation The first contribution is .

Remark 1 and a Cauchy-Schwarz inequality yield
Similarly, Remark 1 and simple inequalities yield 2 ≤ ∆ 1 2 . The dependance of the constant with respect to the data is indicated in the Lemma.
Proof. Consider Lemma 6 and take a small parameter > 0. A Cauchy-Schwarz inequality yields For small enough, the coefficient is less than the coefficient 1 2 in front of the quadratic form in the left hand side of (24). One obtains where we used the same kind of manipulations as in (22), ‖ ‖ Δ = (1) and the CFL condition to bound ∆ with respect to ∆ . All constants depend on ∈ (0, 1). The end of the proof with a discrete Gronwall inequality is standard.
For example one obtains that the minmod scheme (5)- (7) with the minmod flux (9) is quadratically stable and convergent in quadratic norm. This is already a progress with respect to the proof with TVD estimates in [10] where the convergence of the minmod scheme is shown to be not less than (∆ 1 2 ), but in the 1 norm which is weaker than the quadratic norm.

Quadratic estimates in 2D
Now that quadratic stability of flux limiters has been established in 1D, we study a way to use this method for the design of flux limiters in higher dimension. At the end of the construction, it will provide a strategy to bypass the Goodman-Leveque obstruction by using quadratic stability instead of TVD stability. We will use the same approach as in 1D, that is we will modify a 2D Lax-Wendroff linear scheme with limiters and study the resulting scheme. However, in a preliminary stage, we need to define what we denote as the 2D Lax-Wendroff scheme. Due to its simplicity, the scheme is probably not new, even we do not know a reference in the literature.
One considers the same Finite Volume structure as in the Goodman-Leveque work [21] , − , ∆ + where the explicit numerical fluxes at time step are denoted as + 1 2 , and , + 1 2 for all possible and . The mesh size is the same in the horizontal and vertical directions, that is ∆ = ∆ > 0. The coefficients ≥ 0 and ≥ 0 are normalized as in (4). The basic first-order accurate numerical method is the upwind scheme for which the fluxes are Upwind fluxes: up The upwind scheme (28), (29) is stable in quadratic norm (and in all Lebesgue norms) under the CFL condition (8). The upwind scheme is naturally consistant at first order (∆ + ∆ ) = (∆ ) with the advection equation.
A way to derive the Lax-Wendroff scheme is to analyze the numerical dissipation of the upwind scheme with the modified equation technique, and then to subtract the first order numerical dissipation.
Lemma 8. The modified equation of the upwind scheme (28), (29) is or equivalently Proof. With the notations = ( , ) and = ( , +1 ) for all and , one writes the expansions For solution to (28), one has = −( + ) and = ( + ) 2 . Replacing the discrete solution by the exact one in (28) and (29), one gets which is the modified equation (30) after discarding the second order terms. The tensorial nature of the numerical diffusion is better revealed after rearrangement of the second terms in the modified equation. Using + = 1, one has It gives (31).
The modified equation (31) couples an advection equation with an anisotropic diffusion equation. The theory of discrete anisotropic equations made recent progress on the basis of Selling's decomposition of symmetric tensors [4]. In this work, we use a more direct approach to analyze the tensorial nature of the numerical diffusion. The first contribution (1 − ) Δ 2 ( + ) 2 is the diffusion in the direction of the flow, it vanishes for = 1. The second contribution − Δ 2 ( − ) 2 is the diffusion in the direction at angle 3 4 with the horizontal direction. It is independent of the time step and is purely a two-dimensional grid effect. This term has no counterpart in dimension one.
We define the 2D Lax-Wendroff type scheme for which the fluxes are Lax-Wendroff fluxes: It can be rewritten in another equivalent form which is more adapted to our purposes. We note the difference in the direction of the flow and we define modified fluxes Modified fluxes: We also consider the corner difference It allows to defined the scheme So the total incoming flux of (28)-(33) is that is it is equal to the total incoming flux of (34)- (36). The total outgoing flux of (28)-(33) is which is the sum of the total outgoing fluxes of (34)-(36) plus the corner correction in (36). So both schemes are the same. The Lax-Wendroff fluxes (33) are easily justified by the modified equation, so the scheme is second-order accurate (a confirmation is in the numerical Sect. 6). For = 1 and = 0, the scheme is identical to the classical Lax-Wendroff scheme which does not respect the maximum principle.
The 2D spaces for the numerical analysis are extension of the 1D spaces
Lemma 10. One has ∑︀ | , | 2 = ∑︀ | , | 2 − (∆) where the quadratic form is Proof. The structure of the calculation is the same as in the proof of Lemma 1. One decomposes -Calculation of 1 .
One has But one also has , −1 .
One gets by addition -Calculation of 2 .

One has
So the term 3 is Other calculations are immediate.
To analyze the sign of the quadratic form in Lemma 10, we rely on the one-dimensional result of Proposition 2 for the quadratic form (13). Let us define four quantities 1,2,3,4 . The first one 1 concerns the cell-based differences analyzed in the horizontal direction by means of the 1D quadratic form (13) The second one 2 concerns the cell-based differences analyzed in the vertical direction by means of (13) The third one 3 concerns the corner-based differences analyzed in the horizontal direction Finally the fourth one 4 concerns the corner-based differences analyzed in the vertical direction One has the decomposition where the additional term is The term 5 is by definition ( 2 ), while 1,2,3,4 are affine with respect to the Courant number . Moreover 5 is homogeneous of degree 2 with respect to all products ( , ∆ , ) and (︁ )︁ . Therefore one can consider that 5 in (41) is a correction of slightly different nature with respect to 1,2,3,4 . Under the conditions (40) and assuming the CFL condition, one already has 1,2,3,4 ≥ 0 as a corollary of Proposition 2. It remains to study the sign of 5 .
2 is defined at vertices. Then one writes Since + = 1, one can use the identity 2 + 2 = ( + ) 2 + ( − ) 2 to transform the first two lines 5 . One obtains An expansion of the last square yields because the double product is rearranged as ( , ) 2 = ( , ) 2 using the commutation property = . So one has the formula Since it is the sum of three non negative terms, the claim is obtained. Proof. Stability is a consequence of (41) and Lemma 11. To show convergence, it is sufficient to adapt Section 3.2. One considers that the truncation error = + is a sum of the truncation error = (∆ 2 ) of the Lax-Wendroff scheme plus the contribution (43) of the limiters: it denoted as = (︀ , )︀ , like (23) in 1D. The latter is function of all 1 − , and all 1 − + 1 2 , + 1 2 . The identity (24) also holds in 2D and it is sufficient to estimate ( , ) Δ as in (25). We give below the main steps of the proof.
-First step. One writes up, The total residual = + is split in two terms. The first is the residual of the linear Lax-Wendroff scheme so it is (∆ 2 ). The second one is the departure to the Lax-Wendroff residual (compare with (23)) -Second step. One has by rearrangement (note the ∆ in from the sums, in accordance with the definition of the scalar product in 2D) -Third step. From the decomposition (41) and Remark 1, one obtains the bound One has readily that ∑︁ -Fourth step. One bounds (44) with a Cauchy-Schwarz inequality where the right hand side is estimated with (45-46). One obtains -Fifth step. One has more directly ( , ) Δ ≤ ( ) (∆ ) 1 2 ∆ 1 2 . One obtains the 2D equivalent of Lemma 6. It yields the proof after an estimate (27) which is the same as in Theorem 7.

Construction of a 2D minmod limiter
The objective in this section is to show that the above considerations can be used to construct a practical procedure for the definition of the limiters ((39) and (40)). The proposed method relies in a first step on the LBV method [12,13] which is an extension of the TVD criterion on general grid. The LBV framework is able to bring, on the one hand some compactness, and on the other hand some nice algebra for the manipulations of multidimensional limiters. In the sequel, the compactness offered by LVD estimates will not be used, only the algebraic manipulations. This method is well suited to design the limiters ( , ) centered in the cells. The second step is devoted to the definition of the limiters (︁ )︁ for the corner interactions.

Construction of cell-based limiters ,
The LBV semi-norm is defined in [13] for general meshes. It is a discrete 1 norm for the discrete gradient in the direction of the flow. We review the LBV framework then we use it to define the cell-based limiters , s. Another idea will be used to define the corner-based limiters.

The LBV semi-norm on a cartesian mesh
A simple definition of the LBV semi-norm on a cartesian mesh is the following.
Let us assume that the numerical solution at the end of time step verifies for all , a double inequality which can be rewritten as Proof. From (50), one has One obtains the inequality Summation over all indices yields which yields the claim since + = 1.
The upwind scheme is LVD, indeed it admits the explicit formulation , = (1− ) , + ( −1, + , −1 ) which satisfies the criterion of Lemma 14. The two-dimensional Lax-Wendroff fluxes (33) generate a scheme which is not LVD for 0 < < 1: this is evident, since the Lax-Wendroff scheme is not TVD in dimension one which corresponds to = 1 and = 0. The interesting question is to find a constructive method for the design of the classical Finite Volume numerical fluxes on the edges + 1 2 , and , + 1 2 such that the LVD criterion is satisfies. Below we use the method proposed by Lagoutière.
Lemma 15. Assume the usual CFL condition 0 < ≤ 1. Then the upwind fluxes and Lax-Wendroff fluxes (33) satisfy for all , Proof. We analyze the two fluxes separately.
-Upwind flux. One has (︁ which yields the claim by convex combination. Lemma 16. Consider the explicit version of the scheme (28) Assume the numerical fluxes satisfies two family of inequalities which are, firstly the incoming conditions (52) for all , , and secondly the following outgoing conditions for all , (54) Remark 2. Inequalities (52) are called incoming because they concern incoming fluxes − 1 2 , and , − 1 2 in a given cell with indices and . They can be interpreted as a kind of consistency requirement of the incoming flux. Inequalities (54) are called outgoing because they concern outgoing fluxes + 1 2 , and , + 1 2 . Assuming the incoming flux satisfy inequalities (52), then the proof below shows that the conditions (54) for the outgoing flux are necessary and sufficient to satisfy the LVD criterion (49).
Proof of Lemma 16. Consider (52) and (53). Then the first inequality of (54) is equivalent to that is , ≤ , after simplifications. Similarly the second inequality of (54) is equivalent to that is , ≥ , after simplifications.
and by The latter means that the correction terms at corners are not taken into account. They will be reintroduced later.
Lemma 18. The scheme (38), (39) with (55)-(57) satisfies the special form (49) of the maximum principle, is LVD and is stable in quadratic norm. It is formally only first order accurate.
Proof. By definition, the requirements (40) are satisfied so quadratic stability holds. The scheme is only firstorder accurate because of (57), that is the corner contributions miss with respect to the formulation (34)- (36) of the second-order accurate Lax-Wendroff scheme.
To show that it is LVD and satisfies the maximum principle under the form (49), one relies on the verification of (52)-(54). One has This is a convex combination since , ∈ [0, 1/2], so the double inequality (52) holds for the total incoming flux. Moreover one has by construction (55), (56) that Then one has the lower bound from which the first inequality of (54) is deduced. Similar manipulations yield the second inequality of (54). Therefore the scheme is LVD.
Next we desire to show the asymptotic value , ≈ 1 for smooth solutions. This is the first step is showing that the scheme is formally second order, which is one element of the refutation of the Goodman-Leveque obstruction Theorem. The second step will to construct the corner limiter + 1 One also has ∆ − 1 One gets the approximation It is reorganized as Asymptotically for ∆ → 0 and < 1, one has ≥ 1 > 1. Considering (65) one obtains + for ∆ small enough. Other cases are similar so one gets that + ± 1 2 , ∓ 1 2 = 1 for small ∆ . One recovers the Lax-Wendroff scheme at the limit which is second order in space and time.
At the end of this construction, the scheme (38), (39) with the edge limiters (55)-(57) and the corner limiters (65) is convergent in quadratic norm with an order not less than (∆ 1 2 ), is formally second order away from characteristic points and satisfies the maximum principle.

Numerical illustrations
We provide some numerical illustrations of the capabilities of the numerical method that has been developed in the above theoretical Sections. The results show without ambiguity that the (∆ 1 2 ) is pessimistic. We compare the results calculated with the Lax-Wendroff scheme (33)- (36) and with the new nonlinear scheme ((38), (39))+((55)-(57))+(65). The error is reported in 1 , 2 and ∞ norms with a procedure explained in Section 6.1. For the two first test problems, the error for the gradient is also reported in 1 and 2 norms.
The domain of calculation is the academic square with periodic boundary conditions, i.e. the torus T = [0, 1] 2 per . In what follows we report error measurements in function of the mesh size in different norms, and deduce the order of convergence. One observes that the order of convergence of the nonlinear scheme is between 1 and the order of convergence of the Lax-Wendroff scheme for the first three test problems with smooth exact solutions. For the last test for which the initial data is the indicatrix function of a square, the initial data has not the smoothness required by (2). On the contrary it is a BV profile and the solution to the equation must be interpreted in the weak sense. In this case the new nonlinear scheme is more accurate that the Lax-Wendroff scheme, both in 1 and 2 norms. For the nonlinear scheme, the order of convergence for the gradient is always less than it is for the function, which it is not the case for the Lax-Wendroff scheme.

Simple test
The initial data is 0 ( , ) = cos 2 ( + 2 ). The velocity vector is defined by = = 1/2. The final time of observation is = 2, so the final solution is equal to the initial solution, that is ( ) = 0 . The CFL is = 0.25. The asymptotic order of accuracy is approximated with the formula order ≈ log( (∆ ) − log( (∆ /2) log 2 where (∆ ) and log( (∆ /2) are evaluated with the two finest meshes. The gradient is reconstructed with the centered formula. The results are in Table 1.
One observes order 2 with the Lax-Wendroff scheme and on order between 1 and 1.5 for the nonlinear scheme. This was expected because it is a known behavior of the 1D minmod scheme. Indeed the accuracy is formally shown to be equal to 2, but away from characteristic points. Around characteristic points, the accuracy is much less because the limitation procedure is activated to satisfy the maximum principle.

Stationary solution
The velocity vector is still = = 1/2. Here 0 ( , ) = cos 2 ( − ) which is a stationary solution ( ) = 0 for all ≥ 0 since ( + ) 0 ≡ 0. Such non trivial stationary solution cannot exist in 1D. The final time of observation is = 2. The CFL is = 0.5. The results are provided in Table 2. The order of convergence is around 3 for the 2 schemes in all norms, except for the nonlinear scheme for which the order is ≈ 5 2 in maximum norm.

A Gaussian
This problem is more challenging for the Lax-Wendroff scheme which develops natural oscillations because this scheme is not able to respect the maximum principle for steep profiles. The velocity vector is now given by = 1/3 and = 2/3. The CFL is = 0.5. The initial data is 0 ( , ) = exp(−100 2 ) with 2 = ( − 0.5) 2 + ( − 0.5) 2 . A plot of the solutions is displayed in Figure 1 on the coarse mesh at = 3. It shows the oscillation of the Lax-Wendroff scheme. On the contrary the new nonlinear scheme respects perfectly the maximum principle, in accordance with its theoretical properties. The error are displayed in Table 3 with similar conclusions as for the previous tests.

Indicatrix function
The interest of an unconditional respect of the maximum principle is better exemplified with an initial data which is an indicatrix function. Indeed, high-order linear schemes cannot satisfy the maximum principle and in practice, important oscillations are generated. Here we choose 0 ( , ) = 1 if max(| − 0.5|, | − 0.5|) < 0.2 and 0 ( , ) = 0 otherwise. It is the indicatrix function of a square. The other datas are as in the two first tests. The solutions are shown in Figure 2. The errors are reported in Table 4. No convergence is observed in ∞ because the solution is a weak solution with BV regularity only. The nonlinear scheme is more accurate than the Lax-Wendroff scheme in 1 and 2 norms. It seems to convergence at order 2 3 in 1 and 1 3 in 2 . This fact is easy to interpret by interpolation of 2 between 1 and ∞ .

Possible extensions
Corner corrections were developed in this work only for the purposes of having compact notations amenable for the refutation of the Goodman-Leveque obstruction theorem by means of quadratically stable limiters. We do not know if corner corrections are logically necessary for the obtention of quadratic stability of flux limitation techniques. Hereafter we evoke some possibilities to extend the modified equation to more interesting problems.
This equation is relevant in plasma physics [3,14,17,28]. It is easy to check that the equivalent equation of the upwind scheme on a cartesian grid is With this formulation, it is easy to construct a conservative (or divergent) Lax-Wendroff like scheme with variable coefficients by modifying (38). One can for example discretize and at half index positions. More complicated coefficients ( , ) and ( , ) are possible also. The other situation concerns a 3D equation with constant coefficients on a cartesian grid + + + = 0 where , , ≥ 0 and + + = 1. The modified equation of the upwind scheme is It can be rewritten as The development of corner correction techniques in 3D should be possible from this formulation. We also think of using the recent theory [4] in order to get more insights into the tensorial nature of the anisotropic diffusion operator. The methods and proofs were developed on a cartesian grid for the purposes of the simplicity of the mathematical developments. It is reasonable to foresee that transport on unstructured grids can also be addressed with quadratically stable limiters. It will nevertheless require a convenient mathematical apparatus to generalize the quadratic forms to such more challenging configurations [24].
Finally one can reasonably expect to generalize the techniques developed in this work to convex scalar conservations laws such as + ∇ · (︁ a(x) 2 2 )︁ = 0 where a ∈ 1,∞ (R ) and can be divergence free (see [15,16]). Indeed such equations admit a entropy law for smooth solutions under the form 2 + ∇ · (︁ a(x) 3 3 )︁ + 1 2 2 ∇ · a(x) = 0, which is compatible with the quadratic estimates of this work.