Bruno Després
email: despres@ann.jussieu.fr

QUADRATIC STABILITY OF FLUX LIMITERS *

Keywords: 1991 Mathematics Subject Classification. 65M08, 65M12 Flux limiters, quadratic stability, quadratic convergence, Goodman-Leveque Theorem

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 O(∆x 1 2 ) 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 [START_REF] Godunov | A Difference Scheme for Numerical Solution of Discontinuous Solution of Hydrodynamic Equations[END_REF] 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 [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF] 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 d = 1 (1D). The theory of nonlinear TVD schemes is nowadays well established, see [START_REF] Eymard | Finite volume methods[END_REF][START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics-A Practical Introduction[END_REF] and references therein. However, in contrast with the success of the TVD theory in 1D, Goodman and Leveque proved in [START_REF] Goodman | On the accuracy of stable schemes for 2D scalar conservation laws[END_REF] the obstruction Theorem summarized as follows: A TVD numerical scheme with compact stencil for the discretization of conservation laws in space dimension d ≥ 2 is at most first-order accurate, except in trivial cases. The proof [START_REF] Goodman | On the accuracy of stable schemes for 2D scalar conservation laws[END_REF] is by contradiction: one starts from a scheme which satisfies the TVD property and which has a compact stencil in space dimension d = 2 (2D); then one proves that this scheme is L 1 contracting for one dimensional profiles; it turns into the fact that the scheme is just first-order accurate by another important result [START_REF] Harten | On finite-difference approximations and entropy conditions for shocks[END_REF]. It is a non constructive proof which does not propose any way to go beyond this limitation (nevertheless we immediately remark that L 2 stability of linear schemes for 1D advection can already be achieved at any order [START_REF] Despres | Uniform asymptotic stability of Strang's explicit compact schemes for linear advection[END_REF], contrary to linear stability in L 1 norm). To be precise, we follow the literature [START_REF] Goodman | On the accuracy of stable schemes for 2D scalar conservation laws[END_REF][START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF] for which formal second-order and high-order accuracy means accuracy in the sense of local Taylor expansion for exact smooth solutions.

Since then, the development of numerical solvers has pursued its route [START_REF] Barth | Finite Volume methods: foundation and analysis[END_REF][START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics-A Practical Introduction[END_REF] in space dimension d ≥ 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 [START_REF] Goodman | On the accuracy of stable schemes for 2D scalar conservation laws[END_REF]: 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 [START_REF] Roe | Optimum positive linear schemes for advection in two and three dimensions[END_REF] are linear and first order so cannot constitute a solution. The Barth-Jespersen technique [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF] on unstructured meshes guarantees the maximum principle, but the scheme is ultimately first-order accurate. One finds in Barth-Ohleberger [START_REF] Barth | Finite Volume methods: foundation and analysis[END_REF] 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 [START_REF] Godunov | A Difference Scheme for Numerical Solution of Discontinuous Solution of Hydrodynamic Equations[END_REF] is about the maximum principle). 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 [START_REF] Cao | Superconvergence of discontinuous Galerkin method for linear hyperbolic equations in one space dimension[END_REF][START_REF] Lu | An oscillation-free discontinuous Galerkin method for scalar hyperbolic conservation laws[END_REF][START_REF] Shu | Bound-preserving high order finite volume schemes for conservation laws and convection-diffusion equations[END_REF] and therein. Refer to [START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high-order schemes for conservation laws: survey and new developments[END_REF] 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 [START_REF] Guermond | Second-order invariant domain preserving approximation of the compressible Navier-Stokes equations[END_REF][START_REF] Guermond | Positive and Asymptotic Preserving Approximation of the Radiation Transport Equation[END_REF], still without proof of convergence for fully discrete schemes (see for example [START_REF] Guermond | Positive and Asymptotic Preserving Approximation of the Radiation Transport Equation[END_REF][Th. 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.

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 ∂ t u(x, t) + ∇ • (a(x)u(x, t)) = s(x)

(1) where the given velocity field is a ∈ W 1,∞ (R d ) and the source is s ∈ L ∞ (R d ). For convection dominated flows in applied sciences, the calculation of accurate numerical solutions to this equation is of fundamental importance [START_REF] Godlewski | Hyperbolic systems of conservation laws[END_REF][START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF][START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics-A Practical Introduction[END_REF]. If the velocity field is divergent free, that is ∇•a = 0, then the equation is non expansive in all L p (R d ) norms and in particular in quadratic norm (this is

L 2 (R d ) norm). If the divergence is bounded ∇ • a ∈ L ∞ (R d )
the expansion is bounded by classical estimates. An equation similar to (1) can be obtained via the linearization of a truly nonlinear equation ∂ t u tot + ∇ • (b(x)f (u tot )) = 0 where the final equation for the perturbation defined for exemple by

u tot = u ref + εu + O(ε 2 ) can be written as ∂ t u(x, t) + ∇ • (a(x, t)u(x, t)) = 0: the velocity field a(x, t) = b(x)f (u ref (x, t)
) 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 [START_REF] Bouchut | Uniqueness and weak stability for multi-dimensional transport equations with one-sided Lipschitz coefficient[END_REF] 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

u(0) = u 0 ∈ W 3 0 (R d ). ( 2
)
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

∂ t u + p∂ x u + q∂ y u = 0, (x, y, t) ∈ R × R × R + . (3)
Additionally the advection velocity field is normalized p ≥ 0, q ≥ 0 and p + q = 1.

One dimensional configurations aligned with the discretization grid correspond to the limit cases (p, q) = (1, 0) or (p, q) = (0, 1). The equation ( 3) is already relevant for discretization of kinetic equations, for which u = f (x, v, t) ≥ 0 represents the density of particles in function of the space variable, the time variable and the velocity v ∈ R d variable. The correspondance is v = λ(p, q) where λ ≥ 0 is a rescaling factor to account for the normalization p + q = 1. Considering the fierce activity nowadays in the field of numerical plasma physics [START_REF] Bernier | Splitting methods for rotations: application to Vlasov equations[END_REF][START_REF] Crouseilles | Conservative semi-Lagrangian schemes for Vlasov equations[END_REF][START_REF] Einkemmer | Convergence analysis of a discontinuous Galerkin/Strang splitting approximation for the Vlasov-Poisson equations[END_REF][START_REF] Filbet | Comparison of Eulerian Vlasov solvers[END_REF][START_REF] Latu | Field-aligned interpolation for semi-Lagrangian gyrokinetic simulations[END_REF], 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 [START_REF] Després | Numerical Methods for Eulerian and Lagrangian Conservation Laws[END_REF] applied to the numerical discretization of compressible Euler equations (a different technique is [START_REF] Kurganov | Numerical dissipation switch for two-dimensional central-upwind schemes[END_REF]). Radiation transport equations [START_REF] Guermond | Positive and Asymptotic Preserving Approximation of the Radiation Transport Equation[END_REF] 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 11 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 O(∆x 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 performant than the linear Lax-Wendroff for truly weak solutions with BV regularity. Some natural extensions are evoked in last Section 7.

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 [START_REF] Godlewski | Hyperbolic systems of conservation laws[END_REF][START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF][START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics-A Practical Introduction[END_REF]. A generic scheme for the advection equation

∂ t u + ∂ x u = 0 writes U j -U j ∆t + U j+ 1 2 -U j-1 2 ∆x = 0, j ∈ Z, (5)
with

U j+ 1 2 = U j + 1 -ν 2 ϕ j+ 1 2 (U j+1 -U j ) for all j ∈ Z. (6)
In these notations, U j stand for U n j which is the numerical solution at time step t n = n∆t and U j stands for U n+1 j which is the numerical solution at time step t n+1 = (n + 1)∆t. An explicit formulation is

U j = (1 -ν)U j + νU j-1 - ν(1 -ν) 2 ϕ j+ 1 2 (U j+1 -U j ) -ϕ j-1 2 (U j -U j-1 ) . (7)
The scalar ϕ j+ 1 2 is the flux limiter between two consecutive cells (or mesh points). It is well known [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics-A Practical Introduction[END_REF] that if ϕ j+ 1 2 ≡ 0 for all j then one obtains the upwind scheme, if ϕ j+ 1 2 ≡ 1 for all j then one obtains the Lax-Wendroff scheme. The theory of TVD flux limitation has been developed in [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF][START_REF] Sweby | High resolution schemes using flux-limiters for hyperbolic conservation laws[END_REF] and is exposed in [START_REF] Godlewski | Hyperbolic systems of conservation laws[END_REF][START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF][START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics-A Practical Introduction[END_REF]. For example the minmod limiter is defined by

ϕ j+ 1 2 = minmod 1, r j+ 1 2 where r j+ 1 2 = U j -U j-1 U j+1 -U j . (8)
The minmod function is defined by minmod(a, b) = 0 for ab ≤ 0 and minmod(a, b) = sign(a) min(|a|, |b|) for ab ≥ 0.

Contrary to the classical references [START_REF] Godlewski | Hyperbolic systems of conservation laws[END_REF][START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF][START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics-A Practical Introduction[END_REF] where the analysis is l 1 -based or TVD-based, our analysis is based on the quadratic discrete Lebesgue space

l 2 =    U = (U j ) j∈Z such that j∈Z |U j | 2 < ∞    .
In the next preliminary result, the l 2 norm of the numerical solution at the end of the time step is compared with the l 2 norm of the numerical solution at the beginning of the time step. We note ∆ = ∆ j+ 1 2 j∈Z with

∆ j+ 1 2 = U j+1 -U j . Lemma 1. Let U ∈ l 2 . Then one has j∈Z |U j | 2 = j∈Z |U j | 2 -Q ϕ (∆) where Q ϕ (∆) is the quadratic form Q ϕ (∆) = ν(1 -ν) j∈Z |∆ j-1 2 | 2 -ν(1 -ν) j∈Z ϕ j+ 1 2 ∆ j+ 1 2 (1 -ν)∆ j+ 1 2 + ν∆ j-1 2 -ν 2 (1-ν) 2 4 j∈Z ϕ j+ 1 2 ∆ j+ 1 2 -ϕ j-1 2 ∆ j-1 2 2 . ( 9
)
The parameters of the quadratic form are noted ϕ = ϕ j+ 1 2 j∈Z .

Proof. The scheme (5-6) is written explicitly as

U j = U up j - ν(1 -ν) 2 ϕ j+ 1 2 ∆ j+ 1 2 -ϕ j-1 2 ∆ j-1 2
where the upwind scheme is U up j = (1 -ν)U j + νU j-1 and the difference is ∆ j+ 1 2 = U j+1 -U j for all j. Denoting R j = |U j | 2 -|U up j | 2 , one can write

|U j | 2 -|U j | 2 = R j + |U up j | 2 -|U j | 2 = R j -ν(1 -ν)|U j -U j-1 | 2 -ν|U j | 2 + ν|U j-1 | 2 = R j -ν(1 -ν)|∆ j-1 2 | 2 + -ν|U j | 2 + ν|U j-1 | 2 a difference (10)
where

R j = |U j | 2 -|U up j | 2 = 2(U j -U up j )U up j + |U j -U up j | 2 = -ν(1 -ν) ϕ j+ 1 2 ∆ j+ 1 2 -ϕ j-1 2 ∆ j-1 2 U up j + ν 2 (1-ν) 2 4 ϕ j+ 1 2 ∆ j+ 1 2 -ϕ j-1 2 ∆ j-1 2 2 = -ν(1 -ν) ϕ j+ 1 2 ∆ j+ 1 2 U up j+1 -ϕ j-1 2 ∆ j-1 2 U up j +ν(1 -ν)ϕ j+ 1 2 ∆ j+ 1 2 (U up j+1 -U up j ) + ν 2 (1-ν) 2 4 ϕ j+ 1 2 ∆ j+ 1 2 -ϕ j-1 2 ∆ j-1 2 2 = -ν(1 -ν) ϕ j+ 1 2 ∆ j+ 1 2 U up j+1 -ϕ j-1 2 ∆ j-1 2 U up j another difference +ν(1 -ν)ϕ j+ 1 2 ∆ j+ 1 2 (1 -ν)∆ j+ 1 2 + ν∆ j-1 2 + ν 2 (1-ν) 2 4 ϕ j+ 1 2 ∆ j+ 1 2 -ϕ j-1 2 ∆ j-1 2 2 . ( 11
)
With these expressions, one can calculate j∈Z |U j | 2 -|U j | 2 . The differences in [START_REF] Despres | Uniform asymptotic stability of Strang's explicit compact schemes for linear advection[END_REF][START_REF] Després | Generalized Harten formalism and longitudinal variation diminishing schemes for linear advection on arbitrary grids[END_REF] become telescopic by summation, so they vanish and the claim is obtained.

The upwind scheme corresponds to ϕ j+ 1 2 ≡ 0: in this case the quadratic form Q up = Q 0 from ( 9) is non positive for all ∆ provided 0 < ν ≤ 1. One recovers without surprise the quadratic stability of the upwind scheme under CFL condition [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF]. The Lax-Wendroff scheme corresponds to ϕ j+ 1 2 ≡ 1. One obtains he quadratic form

Q lw = Q 1 where Q 1 (∆) = ν(1 -ν) j∈Z |∆ j-1 2 | 2 -∆ j+ 1 2 ((1 -ν)∆ j+ 1 2 + ν∆ j+ 1 2 ) -ν 2 (1-ν) 2 4 j∈Z ∆ j+ 1 2 -∆ j-1 2 2 = ν 2 (1-ν) 2 j∈Z ∆ j+ 1 2 -∆ j-1 2 2 -ν 2 (1-ν) 2 4 j∈Z ∆ j+ 1 2 -∆ j-1 2 2 = ν 2 (1-ν 2 ) 4 j∈Z ∆ j+ 1 2 2 .
One recovers, also without surprise, the quadratic stability of the Lax-Wendroff scheme under CFL condition [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF]. The natural question is to determine if there exists a general condition on the limiters

ϕ = ϕ j+ 1 2 j∈Z such that Q ϕ (∆) ≤ 0 for all ∆ ∈ l 2
. A trivial answer is interpolation of ϕ j+ 1 2 between 0 and 1. If the interpolation is uniform with respect to j, that is ϕ j+ 1 2 = ψ ∈ [0, 1] for all j ∈ Z, then the result is evident. The key fact below is that the interpolation can be taken for ϕ j+ 1 2 independently to ϕ k+ 1 2 for all j = k. The main technical difficulty concerns the quadratic form

H ϕ (∆) = j∈Z |∆ j-1 2 | 2 - j∈Z ϕ j+ 1 2 ∆ j+ 1 2 (1 -ν)∆ j+ 1 2 + ν∆ j-1 2 - ν 4 j∈Z ϕ j+ 1 2 ∆ j+ 1 2 -ϕ j-1 2 ∆ j-1 2 2 . (12)
This quadratic form is a part of Q ϕ (∆) since one has the formula

Q ϕ (∆) = ν(1 -ν)   H ϕ (∆) + ν 2 4 j∈Z ϕ j+ 1 2 ∆ j+ 1 2 -ϕ j-1 2 ∆ j-1 2 2   . (13)
Proposition 2. Assume the CFL condition [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF] and assume 0 ≤ ϕ j+ 1 2 ≤ 1 for all j ∈ Z. Then H ϕ (∆) ≥ 0 for all ∆ ∈ l 2 .

Proof. The condition ∆ ∈ l 2 guarantees the convergence of the sum in [START_REF] Després | Lax Theorem and Finite Volume schemes[END_REF], so H ϕ (∆) is finite. The proof uses a rearrangement of [START_REF] Després | Lax Theorem and Finite Volume schemes[END_REF], which is valid for ∆ ∈ l 2 . Since H ϕ (ν) is linear with respect to ν, we note

H ϕ (ν) = α+νβ. One has α = j∈Z (1 -ϕ j+ 1 2 )|∆ j+ 1 2 | 2 ≥ 0 and β = j∈Z ϕ j+ 1 2 ∆ j+ 1 2 ∆ j+ 1 2 -∆ j-1 2 -1 4 j∈Z ϕ j+ 1 2 ∆ j+ 1 2 -ϕ j-1 2 ∆ j-1 2 2 = j∈Z ϕ j+ 1 2 -1 2 ϕ 2 j+ 1 2 ∆ 2 j+ 1 2 + j∈Z -ϕ j+ 1 2 + 1 2 ϕ j+ 1 2 ϕ j-1 2 ∆ j+ 1 2 ∆ j-1 2 . Consider γ = 1 4 j∈Z ϕ j+ 1 2 ∆ j+ 1 2 -(2 -ϕ j-1 2 )∆ j-1 2 2 - j∈Z (1 -ϕ j+ 1 2 )∆ j+ 1 2 2 .
An expansion shows that

γ = j∈Z 1 4 ϕ 2 j+ 1 2 + 1 4 (2 -ϕ j+ 1 2 ) 2 -(1 -ϕ j+ 1 2 ) 2 ∆ 2 j+ 1 2 -1 2 j∈Z ϕ j+ 1 2 (2 -ϕ j-1 2 )∆ j+ 1 2 ∆ j-1 2 = j∈Z ϕ j+ 1 2 -1 2 ϕ 2 j+ 1 2 ∆ 2 j+ 1 2 + j∈Z -ϕ j+ 1 2 + 1 2 ϕ j+ 1 2 ϕ j-1 2 ∆ j+ 1 2 ∆ j-1 2 . So β = γ. Therefore H ϕ (∆) = α + νγ = j∈Z (1 -ϕ j+ 1 2 ) 1 -ν(1 -ϕ j+ 1 2 ) ∆ j+ 1 2 2 + ν 4 j∈Z ϕ j+ 1 2 ∆ j+ 1 2 -(2 -ϕ j-1 2 )∆ j-1 2 2 . ( 14
)
The conditions of the claim insure H ϕ (∆) ≥ 0. Theorem 3. Assume the conditions of Proposition 2. Then the scheme is stable in quadratic norm:

j∈Z |U j | 2 ≤ j∈Z |U j | 2 .
Proof. The identity (13) yields the claim.

Remark 1 (Control of discrete gradients). It is worthwhile for further developments to notice that [START_REF] Einkemmer | Convergence analysis of a discontinuous Galerkin/Strang splitting approximation for the Vlasov-Poisson equations[END_REF][START_REF] Eymard | Finite volume methods[END_REF] combined with a Cauchy-Schwarz inequality yield the inequality

j∈Z (1 -ϕ j+ 1 2 )|∆ j+ 1 2 | 2 + |∆ j+ 1 2 -∆ j-1 2 | 2 ≤ C(ν)Q ϕ (∆) (15)
for some constant C(ν) > 0 which depends of ν ∈ (0, 1). This estimate shows a quadratic control of the second discrete gradient

∆ j+ 1 2 -∆ j-1 2 .
It also shows a quadratic control of the first gradient, multiplied by 1 -ϕ j+ 1 2 .

Convergence in 1D

The convenient space for the numerical analysis is the l 2 space equipped with a quadratic norm weighted with the mesh size

V ∆x =    U = (U j ) ∈ l 2 with norm U 2 ∆x = ∆x j∈Z |U j | 2 < ∞    . The associated scalar product is denoted (U, V ) ∆x = ∆x j∈Z U j V j for U, V ∈ V ∆x .
The strategy of the proof of convergence relies on two ideas. The first idea corresponds to the semi-discrete limit regime which correspond to O(ν) terms in Q ϕ (∆): involved manipulations are similar to classical ones in Finite Volume techniques [START_REF] Eymard | Finite volume methods[END_REF][START_REF] Shu | Bound-preserving high order finite volume schemes for conservation laws and convection-diffusion equations[END_REF]. 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 ∆t → 0 or ν → 0 in [START_REF] Bouchut | Uniqueness and weak stability for multi-dimensional transport equations with one-sided Lipschitz coefficient[END_REF][START_REF] Cao | Superconvergence of discontinuous Galerkin method for linear hyperbolic equations in one space dimension[END_REF]. It writes

U j (t) + U j (t) -U j-1 (t) ∆x + ∆x 2 ϕ j+ 1 2 (t)(U j+1 (t) -U j (t)) -ϕ j-1 2 (t)(U j (t) -U j-1 (t)) ∆x 2 = 0.
We remind the reader that, for simplicity, we consider a smooth initial data [START_REF] Barth | Finite Volume methods: foundation and analysis[END_REF]. Then the solution u is smooth as well. The interpolation of the smooth function

u ∈ W 3 (R × R + ) (with compact support in space) is denoted as V j (t) = u(j∆x, t).
One defines the truncation error r j (t)

r j (t) = V j (t) + V j (t) -V j-1 (t) ∆x + ∆x 2 ϕ j+ 1 2 (t)(V j+1 (t) -V j (t)) -ϕ j-1 2 (t)(V j (t) -V j-1 (t)) ∆x 2 .
It is easy to verify that

r j (t) = O(∆x 2 ) + ∆x 2 (1 -ϕ j+ 1 2 (t))(V j+1 (t) -V j (t)) -(1 -ϕ j-1 2 (t))(V j (t) -V j-1 (t)) ∆x 2 . ( 16
)
The O(∆x 2 ) is uniform with respect to the index j and the time t. Since u is compact in space, then only a finite number of r j (t) can be non zero at a given time t. The number of non zero terms is O(∆x -1 ) uniformly with respect to the time t.

Let us define the numerical error E j (t) = V j (t) -U j (t) which vanishes at initial time E j (0) = 0. The error satisfies

E j (t) + E j (t) -E j-1 (t) ∆x + ∆x ϕ j+ 1 2 (t)(E j+1 (t) -E j (t)) -ϕ j-1 2 (t)(E j (t) -E j-1 (t)) 2∆x 2 = r j (t). ( 17
)
Theorem 4. Let T > 0. Assume 0 ≤ ϕ j+ 1 2 (t) ≤ 1 for all j ∈ Z and 0 ≤ t ≤ T . Then there exists a constant

C > 0 such that E(t) ∆x ≤ C∆x 1 2 .
Proof. Multiply ( 17) by E j (t) and sum over all j. One obtains

d dt 1 2 E(t) 2 ∆x + ∆x 2 2 j∈Z (1 -ϕ j-1 2 (t)) E j (t) -E j-1 (t) ∆x 2 = (E(t), r(t)) ∆x
where the residual can be estimated with [START_REF] Godlewski | Hyperbolic systems of conservation laws[END_REF]. One gets

d dt 1 2 E(t) 2 ∆x + ∆x 2 2 j∈Z (1 -ϕ j-1 2 (t)) Ej (t)-Ej-1(t) ∆x 2 ≤ C∆x 2 E(t) ∆x -∆x 2 2 j∈Z (1 -ϕ j-1 2 (t)) Ej (t)-Ej-1(t) ∆x Vj (t)-Vj-1(t) ∆x ≤ C∆x 2 E(t) ∆x + ∆x 2 4 j∈Z (1 -ϕ j-1 2 (t)) Ej (t)-Ej-1(t) ∆x 2 + ∆x 2 4 j∈Z (1 -ϕ j-1 2 (t)) Vj (t)-Vj-1(t) ∆x 2 . (18)
One has

Vj (t)-Vj-1(t) ∆x
≤ C and only a O(∆x -1 ) terms are non zero. One gets

d dt 1 2 E(t) 2 ∆x ≤ C∆x 2 E(t) ∆x + C∆x.
Then a Gronwall Lemma yields the claim.

The fully discrete scheme

One combines the analysis of Section 3.1 and the inequality [START_REF] Filbet | Comparison of Eulerian Vlasov solvers[END_REF]. The interpolation of the exact solution u on the space-time grid is V n j = u(j∆x, j). The truncation error is r n j

r n j = V n+1 j -V n j ∆t + V n j -V n j-1 ∆x + (1 -ν) 2 ϕ n j+ 1 2 (V n j+1 -V n j ) -ϕ n j-1 2 (V n j -V n j-1 ) ∆x .
It is decomposed in two parts r n = s n + t n . The first part is the truncation error of the Lax-Wendroff scheme

s n j = V n+1 j -V n j ∆t + V n j -V n j-1 ∆x + (1 -ν) 2 (V n j+1 -V n j ) -(V n j -V n j-1 ) ∆x so s n j = O(∆x 2
). The second part comes from the flux limitation

t n j = - (1 -ν) 2 (1 -ϕ n j+ 1 2 )(V n j+1 -V n j ) -(1 -ϕ n j-1 2 )(V n j -V n j-1 ) ∆x . (19)
The error E n j = V n j -U n j vanishes at initial time (that is E 0 j = 0 for all j) and satisfies

E n+1 j = E n j + ∆tr n j where E n j = (1 -ν)E j + νE j-1 -ν(1-ν) 2 ϕ j+ 1 2 (E j+1 -E j ) -ϕ j-1 2 (E j -E j-1
) . One deduces that

E n+1 2 ∆x = E n 2 ∆x + 2∆t E n , r n ∆x + ∆t 2 r n 2 ∆x .
By comparison with [START_REF] Godunov | A Difference Scheme for Numerical Solution of Discontinuous Solution of Hydrodynamic Equations[END_REF], it is convenient to rewrite it as

E n+1 2 ∆x -E n 2 ∆x 2∆t = E n , r n ∆x + ∆t 2 r n 2 ∆x Using Lemma 1 to evaluate E n 2 ∆x -E n 2 ∆x , one obtains E n+1 2 ∆x -E n 2 ∆x 2∆t + 1 2ν Q ϕ n (∆ E,n ) = E n , r n ∆x + ∆t 2 r n 2 ∆x ( 20
)
where the quadratic form is evaluated with respect to ∆ E,n = ∆ E,n

j+ 1 2 j∈Z with ∆ E,n j+ 1 2 = E n j+1 -E n j .
The last term in [START_REF] Guermond | Second-order invariant domain preserving approximation of the compressible Navier-Stokes equations[END_REF] does not pose any problem because by definition

r n ∆x = O(1) so ∆t r n ∆x = O(∆t). The other term is E n , r n ∆x = E n , s n ∆x + E n , t n ∆x ( 21
)
where

E n , s n ∆x ≤ C∆x 2 E n ∆x and E n , t n ∆x = - (1 -ν) 2 j∈Z (1 -ν)E n j + νE n j-1 - ν(1 -ν) 2 ϕ n j+ 1 2 ∆ E,n j+ 1 2 -ϕ n j-1 2 ∆ E,n j-1 2 (22)
× (1 -ϕ n j+ 1 2 )D n j+ 1 2 -(1 -ϕ n j-1 2 )D n j-1 2 with the notation D n j+ 1 2 = V n j+1 -V n j . Lemma 5. One has E n , t n ∆x ≤ C Q ϕ n ∆ E,n 1 2 ∆x 1 2 .
Proof. From [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF] one splits

E n , t n ∆x = G n 1 + G n 2 with G n 1 = - 1 -ν 2 j∈Z (1 -ν)E n j + νE n j-1 (1 -ϕ n j+ 1 2 )D n j+ 1 2 -(1 -ϕ n j-1 2 )D n j-1 2 and G n 2 = ν(1 -ν) 2 4 j∈Z ϕ n j+ 1 2 ∆ E,n j+ 1 2 -ϕ n j-1 2 ∆ E,n j-1 2 (1 -ϕ n j+ 1 2 )D n j+ 1 2 -(1 -ϕ n j-1 2 )D n j-1 2 .
The first contribution is

G n 1 = 1 -ν 2 j∈Z (1 -ν)∆ E,n j+ 1 2 + ν∆ E,n j-1 2 (1 -ϕ n j+ 1 2 )D n j+ 1 2 = - 1 -ν 2 j∈Z (1 -ϕ n j+ 1 2 )∆ E,n j+ 1 2 D n j+ 1 2 -ν 1 -ν 2 j∈Z ∆ E,n j+ 1 2 -∆ E,n j-1 2 (1 -ϕ n j+ 1 2 )D n j+ 1 2 .
Remark 1 and a Cauchy-Schwarz inequality yield

G n 1 ≤ C Q ϕ n ∆ E,n 1 2   j∈Z |D n j+ 1 2 | 2   1 2 ≤ C Q ϕ n ∆ E,n 1 2 ∆x 1 2 .
Similarly, Remark 1 and simple inequalities yield for the second contribution

G n 2 ≤ CQ ϕ n (∆ E,n ) 1 2 ∆x 1 2
which ends the proof. Theorem 6. Let T > 0. Assume the CFL condition [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF] and assume 0 ≤ ϕ n j+ 1 2 ≤ 1 for all j ∈ Z and all n ∈ N such that n∆t ≤ T . Then there exists a constant C(ν) > 0 such that

max n∆t≤T E n ∆x ≤ C(ν)∆x 1 2 .
Proof. Consider Lemma 5 and take a small parameter ε > 0. A Cauchy-Schwarz inequality yields

E n , t n ∆x ≤ CεQ ϕ n ∆ E,n + c ε ∆x.
For ε small enough, the coefficient Cε is less than the coefficient 1 2ν in front of the quadratic form in the left hand side of [START_REF] Guermond | Second-order invariant domain preserving approximation of the compressible Navier-Stokes equations[END_REF]. One obtains

E n+1 2 ∆x -E n 2 ∆x 2∆t ≤ C∆x 2 E n ∆x + C∆x + ∆t 2 r n 2 ∆x ( 23
)
where all constants depend on ν ∈ (0, 1). The end of the proof with a Gronwall inequality is standard.

For example one obtains that the minmod scheme (5-7) with the minmod flux ( 8) is quadratically stable and convergent in quadratic norm. This is already a progress with respect to the proof with TVD estimates in [START_REF] Després | Lax Theorem and Finite Volume schemes[END_REF] where the convergence of the minmod scheme is shown to be not less than O(∆x 1 2 ), but in the L 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 [START_REF] Goodman | On the accuracy of stable schemes for 2D scalar conservation laws[END_REF] U

i,j -U i,j ∆t + p U i+ 1 2 ,j -U i-1 2 ,j ∆x + q U i,j+ 1 2 -U i,j-1 2 ∆x = 0, for all i, j (24)
where the explicit numerical fluxes at time step n are denoted as U i+ 1 2 ,j and U i,j+ 1 2 for all possible i and j. The mesh size is the same in the horizontal and vertical directions, that is ∆x = ∆y > 0.

The basic first-order accurate numerical method is the upwind scheme for which the fluxes are Upwind fluxes:

U up i+ 1 2 ,j = U up i,j+ 1 2 = U i,j . (25)
The upwind scheme (24-25) is stable in quadratic norm (and in all Lebesgue norms) under the CFL condition

ν = ∆t ∆x ≤ 1. ( 26
)
The upwind scheme is naturally consistant at first order O(∆x + ∆t) = O(∆x) 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 7. The modified equation of the upwind scheme (24-25) is

∂ t u + p∂ x u + q∂ y u - p∆x 2 ∂ xx u - q∆x 2 ∂ yy u + ∆t 2 (p∂ x + q∂ y ) 2 u = 0, (27)
or equivalently

∂ t u + p∂ x u + q∂ y u -(1 -ν) ∆x 2 (p∂ x + q∂ y ) 2 u -pq ∆x 2 (∂ x -∂ y ) 2 u = 0. ( 28
)
Proof. With the notations V n j = u(x j , t n ) for all j and n, one writes the expansions

       p V n i,j -V n i-1,j ∆x = p∂ x u n i,j -p∆x 2 ∂ xx u n i,j + O(∆x 2 ), q V n i,j -V n i,j-1 ∆x = q∂ y u n i,j -q∆x 2 ∂ yy u n i,j + O(∆x 2 ), V n+1 i,j -V n i,j ∆t = ∂ t u n i,j + ∆t 2 ∂ tt u n i,j + O(∆t 2 ). ( 29
)
For u solution to [START_REF] Kurganov | Numerical dissipation switch for two-dimensional central-upwind schemes[END_REF], one has ∂ t u = -(p∂ x + q∂ y )u and ∂ tt u = (p∂ x + q∂ y ) 2 u. Plugging ( 29) in [START_REF] Kurganov | Numerical dissipation switch for two-dimensional central-upwind schemes[END_REF][START_REF] Latu | Field-aligned interpolation for semi-Lagrangian gyrokinetic simulations[END_REF] , one gets

∂ t u + p∂ x u + q∂ y u - p∆x 2 ∂ xx u - q∆x 2 ∂ yy u + ∆t 2 (p∂ x + q∂ y ) 2 u = O(∆x 2 + ∆t 2 )
which is the modified equation ( 27) 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. One has -p∂ xx u -q∂ yy u + ν(p∂

x + q∂ y ) 2 u = -(p∂ x + q∂ y ) 2 u + (p 2 -p)∂ xx u + 2pq∂ xy u + (q 2 -q)∂ yy u + ν(p∂ x + q∂ y ) 2 u = -(1 -ν)(p∂ x + q∂ y ) 2 u -pq (∂ xx -2∂ xy + ∂ yy ) u = -(1 -ν)(p∂ x + q∂ y ) 2 u -pq (∂ x -∂ y ) 2 u.
It gives [START_REF] Roe | Optimum positive linear schemes for advection in two and three dimensions[END_REF].

The modified equation ( 28) 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 [START_REF] Bonnans | Second order monotone finite differences discretization of linear anisotropic differential operators[END_REF]. In this work, we use a more direct approach to analyze the tensorial nature of the numerical diffusion. The first contribution (1 -ν) ∆x 2 (p∂ x + q∂ y ) 2 u is the diffusion in the direction of the flow, it vanishes for ν = 1. The second contribution -pq ∆x 2 (∂ x -∂ y ) 2 u 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:

U lw i+ 1 2 ,j = Ui,j +Ui+1,j 2 -ν Ui+1,j -pUi,j -qUi+1,j-1 2 , U lw i,j+ 1 2 = Ui,j +Ui,j+1 2 -ν Ui,j+1-pUi-1,j+1-qUi,j 2
, ( 30
)
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 ∆ i,j = U i,j -pU i-1,j -qU i,j-1 and we define modified fluxes Modified fluxes:

U mod i+ 1 2 ,j = U i,j + 1-ν 2 ∆ i+1,j , U mod i,j+ 1 2 = U i,j + 1-ν 2 ∆U i,j+1 , (31)
We also consider the corner difference

∆ i+ 1 2 ,j+ 1 2 = U i,j+1 -U i+1,j . (32)
It allows to defined the scheme

U i,j -U i,j ∆t + p U mod i+ 1 2 ,j -U mod i-1 2 ,j ∆x + q U mod i,j+ 1 2 -U mod i,j-1 2 ∆x + pq 2 ∆ i-1 2 ,j+ 1 2 -∆ i+ 1 2 ,j-1 2 ∆x = 0 ( 33
)
Lemma 8. The scheme ( 24)-( 30) is equal to the scheme (33)- [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics-A Practical Introduction[END_REF] with the corner correction [START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high-order schemes for conservation laws: survey and new developments[END_REF]. It is second order accurate in space and time and does not respect the maximum principle.

Proof. The total incoming flux of ( 24)-( 30) is

pU lw i-1 2 ,j + qU lw i,j-1 2 = pUi-1,j +qUi,j-1 2 + Ui,j 2 -ν 2 ∆U i,j = pU i-1,j + qU i,j-1 + 1-ν 2 ∆ i,j = pU mod i-1 2 ,j + qU mod i,j-1 2
, so it is equal to the total incoming flux of [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics-A Practical Introduction[END_REF][START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high-order schemes for conservation laws: survey and new developments[END_REF](33).

The total outgoing flux of ( 24)-( 30) is

pU lw i+ 1 2 ,j + qU lw i,j+ 1 2 = Ui,j 2 + pUi+1,j +qUi,j+1 2
-ν 2 p∆ i+1,j -ν 2 q∆U i,j+1 = U i,j + 1-ν 2 p∆U i+1,j + 1-ν 2 q∆ i,j+1 -q 2 (U i,j+1 -pU i-1,j+1 -qU i,j ) - Ui,j 2 + pUi+1,j +qUi,j+1 2
-p 2 (U i+1,j -pU i,j -qU i+1,j-1 ) = pU mod i-1 2 ,j + qU mod i,j-1 2 -pqU i,j + pq 2 U i+1,j-1 + pq 2 U i-1,j+1 = pU mod i-1 2 ,j + qU mod i,j-1 2 + pq 2 ∆ i-1 2 ,j+ 1 2 -∆ i+ 1 2 ,j-1 2
which is the sum of the total outgoing fluxes of (31-33) plus the corner correction in (33). So both schemes are the same. The Lax-Wendroff fluxes [START_REF] Sweby | High resolution schemes using flux-limiters for hyperbolic conservation laws[END_REF] are easily justified by the modified equation, so the scheme is second-order accurate (a confirmation is in the numerical Section 6). For p = 1 and q = 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

l 2 =    U = (U i,j ) i,j∈Z such that ij |U i,j | 2 < ∞    and V ∆x =    U ∈ l 2 with the norm U 2 ∆x = ∆x ij |U i,j | 2    .

2D limiters

The Lax-Wendroff scheme is written (33) as a Finite Volume scheme with additional corner interactions. The introduction of flux limiters ϕ i+ 1 2 ,j , ϕ i,j+ 1 2 and ϕ i+ 1 2 ,j+ 1 2 is done my mimicking the 1D procedure [START_REF] Crouseilles | Conservative semi-Lagrangian schemes for Vlasov equations[END_REF][START_REF] Després | Numerical Methods for Eulerian and Lagrangian Conservation Laws[END_REF]. It yields the scheme

U i,j -U i,j ∆t + p U lim i+ 1 2 ,j -U lim i-1 2 ,j ∆x + q U lim i,j+ 1 2 -U lim i,j-1 2 ∆x + pq 2
C lim i-1 2 ,j+ 1 2 -C lim i+ 1 2 ,j-1 2 ∆x = 0 (34) where      U lim i+ 1 2 ,j = U i,j + 1-ν 2 ϕ i+1,j ∆ i+1,j , U lim i,j+ 1 2 = U i,j + 1-ν 2 ϕ i,j+1 ∆ i,j+1 , C lim i+ 1 2 ,j+ 1 2 = ϕ i+ 1 2 ,j+ 1 2 ∆ i+ 1 2 ,j+ 1 2 .
(

) 35
The cell-based limiter ϕ i,j is the equivalent of ϕ j+ 1 2 is 1D (while ∆ i,j is the equivalent of ∆ j-1 2 ). The cornerbased limiter limiter ϕ i-1 2 ,j+ 1 2 is a new term related to the tensorial nature of numerical diffusion in 2D. Inspired by the quadratic properties in 1D, we study the properties of the scheme when the flux limiters satisfy the bounds

0 ≤ ϕ i,j , ϕ i+ 1 2 ,j+ 1 2 ≤ 1 for all i, j. (36)
We will assume without restriction that U ∈ V ∆x .

Lemma 9. One has ij |U i,j | 2 = ij |U i,j | 2 -Q ϕ (∆)
where the quadratic form is

Q ϕ (∆) = ν(1 -ν) ij |∆ i,j | 2 + νpq ij ∆ i-1 2 ,j+ 1 2 2 -ν(1 -ν) ij ϕ i,j ∆ i,j ((1 -ν)∆ i,j + νp∆ i-1,j + νq∆ i,j-1 ) -νpq ij ϕ i-1 2 ,j+ 1 2 ∆ i-1 2 ,j+ 1 2 (1 -ν)∆ i-1 2 ,j+ 1 2 + νp∆ i-3 2 ,j+ 1 2 + νq∆ i-1 2 ,j-1 2 - ij ν(1 -ν) 2 (pϕ i+1,j ∆ i+1,j + qϕ i-1,j ∆ i-1, j -ϕ i,j ∆ i,j ) + νpq 2 ϕ i-1 2 ,j+ 1 2 ∆ i-1 2 ,j+ 1 2 -ϕ i+ 1 2 ,j-1 2 ∆ i+ 1 2 ,j-1 2 2
Proof. The structure of the calculation is the same as in the proof of Lemma 1. One decompose

Q ϕ (∆) = ij |U i,j | 2 - ij |U up i,j | 2 =Q 1 + ij |U up i,j | 2 - ij |U i,j | 2 =Q 2 . • Calculation of Q 1 . Since U up i,j = (1 -ν)U i,j + νpU i-1,j + νqU i,j-1 , one has ij |U up i,j | 2 = ij (1 -ν) 2 + ν 2 p 2 + ν 2 q 2 |U i,j | 2 +2 ij (1 -ν)νpU i,j U i-1,j + (1 -ν)νqU i,j U i,j-1 + ν 2 pqU i-1,j U i,j-1 .
But one also has

ν(1 -ν) ij |U i,j -pU i-1,j -qU i,j-1 | 2 = ν(1 -ν) ij 1 + p 2 + q 2 |U i,j | 2 -2ν(1 -ν)pU i,j U i-1,j -2ν(1 -ν)qU i,j U i,j-1 + 2ν(1 -ν)pqU i-1,j U i,j-1 .
One gets by addition

|U up i,j | 2 + ν(1 -ν) ij |U i,j -pU i-1,j -qU i,j-1 | = ij (1 -ν) 2 + ν 2 p 2 + ν 2 q 2 + ν(1 -ν) + ν(1 -ν)p 2 + ν(1 -ν)q 2 |U i,j | 2 + 2νpqU i-1,j U i,j-1 = ij (1 -ν) + νp 2 + νq 2 |U i,j | 2 + 2νpqU i-1,j U i,j-1 = ij |U i,j | 2 -νpq ij |U i-1,j -U i,j-1 | 2 . So Q 1 = ν(1 -ν) ij |U i,j -pU i-1,j -qU i,j-1 | 2 + νpq ij |U i-1,j -U i,j-1 | 2 = ν(1 -ν) ij |∆ i,j | 2 + νpq ij ∆ i-1 2 ,j+ 1 2 2 ≥ 0. • Calculation of Q 2 .
One has

Q 2 = -2 ij (U i,j -U up i,j )U up i,j =Q 3 - ij |U i,j -U up i,j | 2
where

U i,j = U up i,j - ν(1 -ν) 2 (pϕ i+1,j ∆ i+1,j + qϕ i-1,j ∆ i-1, j -ϕ i,j ∆ i,j )- νpq 2 ϕ i-1 2 ,j+ 1 2 ∆ i-1 2 ,j+ 1 2 -ϕ i+ 1 2 ,j-1 2 ∆ i+ 1 2 ,j-1 2 .
So the term Q 3 is

Q 3 = - ν(1 -ν) 2 ij (pϕ i+1,j ∆ i+1,j + qϕ i-1,j ∆ i-1, j -ϕ i,j ∆ i,j ) ((1 -ν)U i,j + νpU i-1,j + νqU i,j-1 ) - νpq 2 ij ϕ i-1 2 ,j+ 1 2 ∆ i-1 2 ,j+ 1 2 -ϕ i+ 1 2 ,j-1 2 ∆ i+ 1 2 ,j-1 2
((1 -ν)U i,j + νpU i-1,j + νqU i,j-1 ) .

Rearrangement yields

Q 3 = ν(1 -ν) 2 ij ϕ i,j ∆ i,j ((1 -ν)∆ i,j + νp∆ i-1,j + νq∆ i,j-1 ) + νpq 2 ij ϕ i-1 2 ,j+ 1 2 ∆ i-1 2 ,j+ 1 2 (1 -ν)∆ i-1 2 ,j+ 1 2 + νp∆ i-3 2 ,j+ 1 2 + νq∆ i-1 2 ,j-1 2 .
Other calculations are immediate.

To analyze the sign of the quadratic form in Lemma 9, we rely on the one-dimensional result of Proposition 2 for the quadratic form [START_REF] Després | A longitudinal variation diminishing estimate for linear advection on arbitrary grids[END_REF]. Let us define four quantities H 1,2,3,4 . The first one H 1 concerns the cell-based differences analyzed in the horizontal direction by means of the 1D quadratic form ( 12)

H 1 = ij |∆ i,j | 2 - ij ϕ i,j ∆ i,j ((1 -ν)∆ i,j + ν∆ i-1,j ) - ν 4 ij |ϕ i,j ∆ i,j -ϕ i-1,j ∆ i-1,j | 2 .
The second one H 2 concerns the cell-based differences analyzed in the vertical direction by means of ( 12)

H 2 = ij |∆ i,j | 2 - ij ϕ i,j ∆ i,j ((1 -ν)∆ i,j + ν∆ i,j-1 ) - ν 4 ij |ϕ i,j ∆ i,j -ϕ i,j-1 ∆ i,j-1 | 2 .
The third one H 3 concerns the corner-based differences analyzed in the horizontal direction

H 3 = ij |∆ i+ 1 2 ,j+ 1 2 | 2 - ij ϕ i+ 1 2 ,j+ 1 2 ∆ i+ 1 2 ,j+ 1 2 (1 -ν)∆ i+ 1 2 ,j+ 1 2 + ν∆ i-1 2 ,j+ 1 2 - ν 4 ij ϕ i+ 1 2 ,j+ 1 2 ∆ i+ 1 2 ,j+ 1 2 -ϕ i-1 2 ,j+ 1 2 ∆ i-1 2 ,j+ 1 2 2
.

Finally the fourth one H 4 concerns the corner-based differences analyzed in the vertical direction

H 4 = ij |∆ i+ 1 2 ,j+ 1 2 | 2 - ij ϕ i+ 1 2 ,j+ 1 2 ∆ i+ 1 2 ,j+ 1 2 (1 -ν)∆ i+ 1 2 ,j+ 1 2 + ν∆ i+ 1 2 ,j-1 2 - ν 4 ij ϕ i+ 1 2 ,j+ 1 2 ∆ i+ 1 2 ,j+ 1 2 -ϕ i+ 1 2 ,j-1 2 ∆ i+ 1 2 ,j-1 2 2
.

One has the decomposition

Q ϕ (∆) = ν(1 -ν) pH 1 + qH 2 + νpq pH 3 + qH 4 + H 5 (37)
where the additional term is

H 5 = ν 2 (1-ν)p 4 ij |ϕ i,j ∆ i,j -ϕ i-1,j ∆ i-1,j | 2 + ν 2 (1-ν)q 4 ij |ϕ i,j ∆ i,j -ϕ i,j-1 ∆ i,j-1 | 2 + ν 2 p 2 q 4 ij ϕ i+ 1 2 ,j+ 1 2 ∆ i+ 1 2 ,j+ 1 2 -ϕ i-1 2 ,j+ 1 2 ∆ i-1 2 ,j+ 1 2 2 + ν 2 pq 2 4 ij ϕ i+ 1 2 ,j+ 1 2 ∆ i+ 1 2 ,j+ 1 2 -ϕ i+ 1 2 ,j-1 2 ∆ i+ 1 2 ,j-1 2 2 - ij ν(1-ν) 2 (pϕ i+1,j ∆ i+1,j + qϕ i-1,j ∆ i-1,j -ϕ i,j ∆ i,j ) + νpq 2 ϕ i-1 2 ,j+ 1 2 ∆ i-1 2 ,j+ 1 2 -ϕ i+ 1 2 ,j-1 2 ∆ i+ 1 2 ,j-1 2 2 (38)
The term H 5 is by definition O(ν 2 ), while H 1,2,3,4 are affine with respect to the Courant number ν. Moreover H 5 is homogeneous of degree 2 with respect to all products (ϕ i,j ∆ i,j ) ij and ϕ i+ 1 2 ,j-

1 2 ∆ i+ 1 2 ,j-1 2 ij
. Therefore one can consider that H 5 in (37) is a correction of slightly different nature with respect to H 1,2,3,4 . Under the conditions (36) and assuming the CFL condition, one already has H 1,2,3,4 ≥ 0 as a corollary of Proposition 2.

It remains to study the sign of H 5 .

Lemma 10. Under the previous conditions, one has H 5 ≥ 0.

Proof. Let us begin with notations which are useful to synthesize the structure of H 5 . For a = (a i,j ) ∈ l 2 defined at centers, one defines the operator L such that (La) i,j = a i,j -pa i-1,j -qa i,j-1 and the operator M such that (M a) i,j = a i-1,j -a i,j-1 . For b = (b

i+ 1 2 ,j+ 1
2 ) ∈ l 2 defined at vertices, one defines with similar notations (Lb)

i+ 1 2 ,j+ 1 2 = b i+ 1 2 ,j+ 1 2 -pb i-1 2 ,j+ 1 2 -qb i+ 1 2 ,j-1 2 and (M b) i+ 1 2 ,j+ 1 2 = b i-1 2 ,j+ 1 2 -b i+ 1 2 ,j-1 2 . Let T the translation operator such that a = T b means a i,j = b i+ 1 2 ,j+ 1 2 .
The adjoints operators L t and M t are defined with respect to the natural quadratic scalar product in l 2 . One has for example (L t a) i,j = a i,j -pa i+1,j -qa i,j+1 .

One uses a = (a i,j ) where a i,j = ϕ i,j ∆ i,j is defined at centers and b 1 2 is defined at vertices. Then one writes

= (b i+ 1 2 ,j+ 1 2 ) where b i+ 1 2 ,j+ 1 2 = ϕ i+ 1 2 ,j+ 1 2 ∆ i+ 1 2 ,j+
H 5 = ν 2 (1-ν) 4 ij p|a i,j -a i+1,j | 2 + q|a i,j -a i,j+1 | 2 + ν 2 pq 4 ij p|b i,j -b i-1,j | 2 + q|b i,j -b i,j-1 | 2 - ν 2 4 (1 -ν)L t a -pqT M b 2 l 2 .
One uses the identity pα 2 + qβ 2 = (pα + qβ) 2 + pq(α -β) 2 to transform the first two lines H 5 . One obtains

H 5 = ν 2 (1-ν) 4 L t a 2 l 2 + pq M t a 2 l 2 + ν 2 pq 4 Lb 2 l 2 + pq M b 2 l 2 -ν 2 4 (1 -ν)L t a -pqT M b 2 l 2 .
An expansion of the last square yields

(1 -ν)L t a -pqT M b 2 l 2 = (1 -ν) 2 L t a 2 l 2 -2(1 -ν)pq(L t a, T M b) l 2 + p 2 q 2 M b 2 l 2 = (1 -ν) 2 L t a 2 l 2 -2(1 -ν)pq(M t a, T Lb) l 2 + p 2 q 2 M b 2 l 2
because the double product is rearranged as (L t a, T M b) l 2 = (M t a, T Lb) l 2 using the commutation property LT M = M T L. So one has the formula

H 5 = ν 3 (1 -ν) 4 L t a 2 l 2 + ν 2 (1 -ν)pq 4 M t a 2 l 2 + ν 2 pq 4 Lb 2 l 2 + ν 2 (1 -ν)pq 2 (M t a, T Lb) l 2 = ν 3 (1 -ν) 4 L t a 2 l 2 + ν 3 (1 -ν)pq 4 M t a 2 l 2 + ν 2 pq 4 (1 -ν)M t a + T Lb 2 l 2 .
Since it is the sum of three non negative terms, the claim is obtained.

Theorem 11. Under the CFL condition [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF], the scheme (35-36) is quadratically stable and is convergent with an order at least O(∆x 1 2 ) in quadratic norm.

Proof. Stability is a consequence of (37) and Lemma 10. To show convergence, it is sufficient to adapt Section 3.2. One considers that the truncation error r n = s n -t n is a sum of the truncation error s n = O(∆x 2 ) of the Lax-Wendroff scheme plus the contribution (39) of the limiters: it denoted as t n = t n i,j , like [START_REF] Goodman | On the accuracy of stable schemes for 2D scalar conservation laws[END_REF] in 1D. The latter is function of all 1 -ϕ i,j and all 1 -ϕ i+ 1 2 ,j+ 1 2 . The identity (20) also holds in 2D and it is sufficient to estimate (E n , t n ) ∆x as in [START_REF] Guermond | Positive and Asymptotic Preserving Approximation of the Radiation Transport Equation[END_REF]. We give below the main steps of the proof.

• First step. One writes E up,n i,j = (1 -ν)E n i,j + νpE n i-1,j + νqE n i,j-1 , E n i,j = E up,n i,j + F i,j and
F n i,j = - ν(1 -ν) 2 pϕ i+1,j ∆ E,n i+1,j + qϕ i-1,j ∆ E,n i-1,j -ϕ i,j ∆ E,n i,j - νpq 2 ϕ i-1 2 ,j+ 1 2 ∆ E,n i-1 2 ,j+ 1 2 -ϕ i+ 1 2 ,j-1 2 ∆ E,n i+ 1 2 ,j-1 2 .
The total residual r n = s n + t n is split in two terms. The first is the residual of the linear Lax-Wendroff scheme so it is O(∆x 2 ). The second one is the departure to the Lax-Wendroff residual (compare with (19))

t n i,j = - ν(1 -ν) 2 p(1 -ϕ i+1,j )D i+1,j + q(1 -ϕ i-1,j )D i-1,j -(1 -ϕ i,j )D i,j ∆x (39)
- νpq 2 (1 -ϕ i-1 2 ,j+ 1 2 )D i-1 2 ,j+ 1 2 -(1 -ϕ i+ 1 2 ,j-1 2 )D i+ 1 2 ,j-1 2 ∆x .
One has

E n , t n ∆x = (E up,n , t n ) ∆x + (F n , t n ) ∆x .
• Second step. One has by rearrangement

(E up,n , t n ) ∆x = ν(1 -ν) 2 ij (1 -ν)∆ E,n i,j + νp∆ E,n i-1,j + νq∆ E,n i,j-1 (1 -ϕ i,j )D i,j + νpq 2 ij (1 -ν)∆ E,n i+ 1 2 ,j-1 2 + νp∆ E,n i-1 2 ,j-1 2 + νq∆ E,n i+ 1 2 ,j-3 2 (1 -ϕ i+ 1 2 ,j-1 2 )D i+ 1 2 ,j-1 2
• Third step. From the decomposition (37) and Remark 1, one obtains the bound

ij (1 -ϕ n i,j )|∆ E,n i,j | 2 + (1 -ϕ n i+ 1 2 ,j+ 1
2

)|∆ E,n i+ 1 2 ,j+ 1 2 | 2 +p|∆ E,n i,j -∆ E,n i-1,j | 2 + q|∆ E,n i,j -∆ E,n i,j-1 | 2 +p|∆ E,n i+ 1 2 ,j+ 1 2 -∆ E,n i-1 2 ,j+ 1 2 | 2 + q|∆ E,n i+ 1 2 ,j+ 1 2 -∆ E,n i+ 1 2 ,j-1 2 | 2 ≤ C(ν)Q ϕ n (∆ E,n ) • Fourth step. One has the fact that |D i,j | + D i+ 1 2 ,j-1 2
≤ C∆x with only O(∆x -1 ) being non zero. So one obtains by a Cauchy-Schwarz inequality that

(E up,n , t n ) ∆x ≤ C(ν)Q ϕ n (∆ E,n ) 1 2 ∆x 1 2 . • Fifth step. One has more directly (F n , t n ) ∆x ≤ C(ν)Q ϕ n (∆ n ) 1 2 ∆x 1 2
. One obtains the 2D equivalent of Lemma 5. It yields the proof after an estimate [START_REF] Harten | On finite-difference approximations and entropy conditions for shocks[END_REF] which is the same as in Theorem 6.

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 (35-36). The proposed method relies in a first step on the LBV method [START_REF] Després | Generalized Harten formalism and longitudinal variation diminishing schemes for linear advection on arbitrary grids[END_REF][START_REF] Després | A longitudinal variation diminishing estimate for linear advection on arbitrary grids[END_REF] which is an extension of the TVD criterion on general grid. This method is well suited to design the limiters (ϕ i,j ) centered in the cells. The second step is devoted to the definition of the limiters ϕ i+ 1 2 ,j+ 1 2 for the corner interactions.

Construction of cell-based limiters ϕ i,j

The LBV semi-norm is defined in [START_REF] Després | A longitudinal variation diminishing estimate for linear advection on arbitrary grids[END_REF] for general meshes. It is a discrete L 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 ϕ i,j 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. Definition 12. On a cartesian mesh, the LBV semi-norm is

|U | LBV = ∆x 2 i,j U i,j -pU i-1,j -qU i,j-1 ∆x = ∆x i,j |U i,j -pU i-1,j -qU i,j-1 | . ( 40
)
The interest of the LBV semi-norm is that it is endowed with an appropriate multidimensional Harten calculus [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF]. Let us consider the local minimal value and the local maximal value m i,j = min (U i,j , pU i-1,j + qU i,j-1 ) and M i,j = max (U i,j , pU i-1,j + qU i,j-1 ) .

Let us assume that the numerical solution at the end of time step verifies for all i, j a double inequality

m i,j ≤ U i,j ≤ M i,j , (41)
which can be rewritten as

U i,j = (1 -d i,j )U i,j + d i,j (pU i-1,j + qU i,j-1 ) (42)
where 0

≤ d i,j ≤ 1. ( 43
)
Lemma 13. Assume (41) for all i, j.

Then |U | LBV ≤ |U | LBV .
Proof. From (42), one has

U i,j -pU i-1,j + qU i,j-1 = U i,j -(pU i-1,j + qU i,j-1 ) - d i,j (U i,j -(pU i-1,j + qU i,j-1 )) + d i-1,j p (U i-1,j -(pU i-2,j + qU i-1,j-1 )) + d i,j-1 q (U i,j-1 -(pU i-1,j-1 + qU i,j-2 )) = (1 -d i,j ) (U i,j -(pU i-1,j + qU i,j-1 )) + pd i-1,j (U i-1,j -(pU i-2,j + qU i-1,j-1 )) + qd i,j-1 (U i,j-1 -(pU i-1,j-1 + qU i,j-2 )) .

One obtains the inequality

U i,j -pU i-1,j + qU i,j-1 ≤ (1 -d i,j ) |U i,j -(pU i-1,j + qU i,j-1 )| + pd i-1,j |U i-1,j -(pU i-2,j + qU i-1,j-1 )| + qd i,j-1 |U i,j-1 -(pU i-1,j-1 + qU i,j-2 )| .

Summation over all indices yields

i,j U i,j -pU i-1,j + qU i,j-1 ≤ i,j ((1 -d i,j ) + pd i,j + qd i,j ) |U i,j -(pU i-1,j + qU i,j-1 )| which yields the claim since p + q = 1.
The upwind scheme is LVD, indeed it admits the explicit formulation U i,j = (1-ν)U i,j +ν (pU i-1,j + qU i,j-1 ) which satisfies the criterion of Lemma 13. The two-dimensional Lax-Wendroff fluxes (30) 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 p = 1 and q = 0. The interesting question is to find a constructive method for the design of the classical Finite Volume numerical fluxes on the edges U n i+ 1 2 ,j and U n i,j+ 1 2 such that the LVD criterion is satisfies. Below we use the method proposed by Lagoutière. Lemma 14. Assume the usual CFL condition 0 < ν ≤ 1. Then the upwind fluxes and Lax-Wendroff fluxes [START_REF] Sweby | High resolution schemes using flux-limiters for hyperbolic conservation laws[END_REF] satisfy for all i, j m i,j ≤ pU i-1 2 ,j + qU i,j-1 2 ≤ M i,j .

Proof. Evident.

Lemma 15. Consider the explicit version of the scheme [START_REF] Kurganov | Numerical dissipation switch for two-dimensional central-upwind schemes[END_REF] U i,j = U i,j -ν pU i+ 1 2 ,j + qU i,j+ The LVD requirement is restricted to the definition ϕ i,j and is a poor interest for the design of ϕ i+ 1 2 ,j+ 1 2 . This is why a different method is used to construct explicitly the corner-based limiters in a way such that ϕ i+ 1 2 ,j+ 1 2 ≈ 1 for smooth solutions. Since it is important for many applications, we concentrate hereafter on the fulfillment of the maximum principle for the numerical solution at the end of the time step (as in many references cited in the introduction). That is we ask that

m abs i,j = min (U i,j , U i-1,j , U i,j-1 ) ≤ U i,j ≤ M abs i,j = max (U i,j , U i-1,j , U i,j-1 ) (51)
where U i,j at the end of the time step if obtained from the numerical solution defined at (47-48) which is denoted as U i,j . That is U i,j takes of the standard fluxes through the edges while U i,j incorporates the corner corrections. Note that m abs i,j ≤ m i,j ≤ U i,j ≤ M i,j ≤ M abs i,j . The formula for U i,j is

U i,j -U i,j ∆t + p U lim i+ 1 2 ,j -U lim i-1 2 ,j ∆x + q U lim i,j+ 1 2 -U lim i,j-1 2 ∆x = 0. ( 52
)
It is followed by corner corrections

U i,j -U i,j ∆t + pq 2
C lim i-1 2 ,j+ 1 2 -C lim i+ 1 2 ,j-1 2 ∆x = 0. (53)
One gets from (52-53)

C lim i-1 2 ,j+ 1 2 -C lim i+ 1 2 ,j-1 2 = 2 pqν U i,j -U i,j .
To fulfill (51) we require the sufficient conditions

   1 pqν U i,j -M abs i,j ≤ C lim i-1 2 ,j+ 1 2 ≤ 1 pqν U i,j -m abs i,j , 1 pqν U i,j -M abs i,j ≤ -C lim i+ 1 2 ,j-1 2 ≤ 1 pqν U i,j -m abs i,j . (54)
It leads to a definition of the corner-based limiters

ϕ i-1 2 ,j+ 1 2 = min ϕ + i-1 2 ,j+ 1 2 , ϕ - i-1 2 ,j+ 1 2 (55)
where

       ϕ + i-1 2 ,j+ 1 2 = minmod 1, 1 pqν ( Ui,j -m abs i,j ) ∆ i-1 2 ,j+ 1 2 + minmod 1, 1 pqν ( Ui,j -M abs i,j ) ∆ i-1 2 ,j+ 1 2 , ϕ - i+ 1 2 ,j-1 2 = minmod 1, 1 pqν ( Ui,j -m abs i,j ) -∆ i+ 1 2 ,j-1 2 + minmod 1, 1 pqν ( Ui,j -M abs i,j ) -∆ i+ 1 2 ,j-1 2 . ( 56
)
Note that U i,j -m abs i,j ≥ 0 and U i,j -M abs i,j ≤ 0. Therefore one term vanishes systematically in both lines of (56).

Lemma 19. The new scheme with corner limitation (56) satisfies the maximum principle. Assume the solution is locally not an extrema

|∂ x u(x i , y j )| + |∂ y u(x i , y j )| > 0.
Then it is formally second-order.

Proof. For the simplicity of the analysis, we consider that the previous stage for the construction of U i,j is performed with optimal limitation, that is ϕ i,j ≡ 1. Then U i,j corresponds to an integration of the the modified equation [START_REF] Roe | Optimum positive linear schemes for advection in two and three dimensions[END_REF], where the diffusion in the direction of the flow is integrated at second-order. That is the transverse diffusion remains. So one can write

U i,j = u(x i , y j ) -∆xν(p∂ x + q∂ y )u(x i , y j ) + pq ∆x 2 ν 2 (∂ x -∂ y ) 2 u(x i , y j ) + O(∆x 3 ).
Let us analyze the first ratio R = 1 pqν ( Ui,j -m abs i,j )

∆ i-1 2 ,j+ 1 2
. For the simplicity of the analysis we assume that ∂ y u(x i , y j ) > ∂ x u(x i , y j ) which is not a restriction. One can write m abs i,j = u(x i , y j ) + ∆x min (0, -∂ x u(x i , y j ), -∂ y u(x i , y j )) + O(∆x 2 ) that is m abs i,j = u(x i , y j ) -∆x max (0, ∂ y u(x i , y j )) + O(∆x 2 ). So U i,j -m abs i,j = ∆x (max (0, ∂ y u(x i , y j )) -ν(p∂ x + q∂ y )u(x i , y j )) + O(∆x 2 ). One also has

∆ i-1 2 ,j+ 1 2 = -∆x (∂ x -∂ y ) u(x i , y j ) + O(∆x 2 ) One gets the approximation R = max (0, ∂ y u(x i , y j )) -ν(p∂ x + q∂ y )u(x i , y j ) pqν (∂ y u(x i , y j ) -∂ x u(x i , y j )) + O(∆x).
It can be reorganized as R = 1 q + max (0, ∂ y u(x i , y j )) -ν∂ y u(x i , y j ) pqν (∂ y u(x i , y j ) -∂ x u(x i , y j )) + O(∆x)

Asymptotically for ∆x → 0 and q < 1, one has R ≥ 1 q > 1. Considering (56) one obtains ϕ

+ i-1 2 ,j+ 1 2 = 1
for ∆x small enough. Other cases are similar so one gets that ϕ + i± 1 2 ,j∓ 1 2 = 1 for small ∆x. 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 (34-35) with the edge limiters (47-49) and the corner limiters (56) is convergent in quadratic norm with an order not less than O(∆x 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 O(∆x 1 2 ) is pessimistic. We compare the results calculated with the Lax-Wendroff scheme (30-33) and with the new nonlinear scheme (34-35)+(47-49)+(56). The error is reported in L 1 , L 2 and L ∞ norms with a procedure explained in Section 6.1.

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. It seems to be the price to pay to obtain the maximum principle. However 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 [START_REF] Barth | Finite Volume methods: foundation and analysis[END_REF]. 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 L 1 and L 2 norms.

Simple test

The initial data is u 0 (x, y) = cos 2π(x + 2y). The velocity vector is defined by p = q = 1/2. The final time of observation is T = 2, so the final solution is equal to the initial solution, that is u(T ) = u 0 . The CFL is ν = 0.25. The asymptotic order of accuracy is approximated with the formula order ≈ log(E(∆x) -log(E(∆x/2) log 2 .

where E(∆x) and log(E(∆x/2) are evaluated with the two finest meshes. 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. 1. Error measurements for u 0 (x, y) = cos 2π(x + 2y).

LW L 1 LW L 2 LW L ∞ NL L 1 NL L 2 NL L ∞ ∆x =1/

Stationary solution

The velocity vector is still p = q = 1/2. Here u 0 (x, y) = cos 2π(x -y) which is a stationary solution u(t) = u 0 for all t ≥ 0 since (p∂ x + q∂ y )u 0 ≡ 0. Such non trivial stationary solution cannot exist in 1D. The final time of observation is T = 2. The CFL is ν = 0.25. 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. 2. Error measurements for u 0 (x, y) = cos 2π(x -y).

LW L 1 LW L 2 LW L ∞ NL L 1 NL L 2 NL L ∞ ∆x =1/

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 p = 1/3 and q = 2/3. The CFL is ν = 0.5. The initial data is u 0 (x, y) = exp(-100r 2 ) with r 2 = (x -0.5) 2 + (y -0.5) 2 . A plot of the solutions is displayed in Figure 1 on the coarse mesh at T = 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. 3. Error measurements for the Gaussian.

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 u 0 (x, y) = 1 if max(|x -0.5|, |y -0.5|) < 0.2 and u 0 (x, y) = 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 L ∞ because the solution is a weak solution with BV regularity only. The nonlinear scheme is more accurate than the Lax-Wendroff scheme in L 1 and L 2 norms. It seems to convergence at order 2 3 in L 1 and 1 3 in L 2 . This fact is easy to interpret by interpolation of L 2 between L 1 and L ∞ .

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 known if corner corrections are logically necessary for the obtention of quadratic stability of flux limitation techniques. Hereafter we evoke two possibilities to extend the modified equation to more interesting problems.

Let us firstly consider the 2D equation with smooth variable coefficients

∂ t u + ∂ x (a(y)u) + ∂ y (b(x)u) = ∂ t u + (a(y)∂ x + b(x)∂ y b(x))u = 0.
This equation is relevant in plasma physics [START_REF] Bernier | Splitting methods for rotations: application to Vlasov equations[END_REF][START_REF] Einkemmer | Convergence analysis of a discontinuous Galerkin/Strang splitting approximation for the Vlasov-Poisson equations[END_REF][START_REF] Filbet | Comparison of Eulerian Vlasov solvers[END_REF][START_REF] Latu | Field-aligned interpolation for semi-Lagrangian gyrokinetic simulations[END_REF]. 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 (34). One can for exemple discretize p and q at half index positions. More complicated coefficients a(x, y) and b(x, y) are possible also.

The other situation concerns a 3D equation with constant coefficients on a cartesian grid ∂ t u + p∂ x u + q∂ y u + r∂ z u = 0 where p, q, r ≥ 0 and p + q + r = 1. The modified equation of the upwind scheme is The development of corner correction techniques in 3D should be possible from this formulation. We also think of using the recent theory [START_REF] Bonnans | Second order monotone finite differences discretization of linear anisotropic differential operators[END_REF] in order to get more insights into the tensorial nature of the anisotropic diffusion operator.

∂ t u + (
The methods and proofs were developed on a cartesian grid, but for the sole 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.

Figure 2 .

 2 Figure2. Advection of square indicatrix function. On the left calculated with the linear Lax-Wendroff scheme. On the right calculated with limitation of the fluxes. Th maximum principle is not respected with the linear scheme. It is satisfied with the limitation procedure.

∂

  t u + (a(y)∂ x + b(x)∂ y )u -∆x 2 ∂ x (|a(y)|∂ x u) -∆x 2 ∂ y (|b(x)|∂ y u) + ∆t 2 (∂ x a(y) + ∂ y b(x)) 2 = 0.

  Figure1. Advection of Gaussian. On the left calculated with the linear Lax-Wendroff scheme. On the right calculated with limitation of the fluxes. Th maximum principle is not respected with the linear scheme. It is satisfied with the limitation procedure.LW L 1 LW L 2 LW L ∞ NL L 1 NL L 2 NL L ∞

	0.8				0.6
	0.5 0.6 0.7				0.4 0.5
	0.2 0.3 0.4				0.2 0.3
	-0.1 0 0.1				0 0.1
	-0.2				-0.1
	0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x	0	0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 y		0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x	0	0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 y
	∆x =1/40 0.823	0.499	0.438	0.461	0.386	0.454
	∆x =1/80 0.221	0.182	0.182	0.195	0.161	0.186
	∆x =1/160 0.0557 0.0492 0.0535	0.0594 0.0517 0.0731
	∆x =1/320 0.0139 0.0124 0.0136	0.0170 0.0164 0.0312
	order ≈	2		1.98	1.97	1.80	1.65	1.22
	Table

Table 4 .

 4 LW L 1 LW L 2 LW L ∞ NL L 1 NL L 2 NL L ∞Error measurements for the indicatrix function.

	∆x =1/20 0.856	0.487	0.773	0.718	0.479	0.767
	∆x =1/40 0.713	0.470	0.857	0.454	0.381	0.791
	∆x =1/80 0.508	0.381	0.860	0.280	0.297	0.801
	∆x =1/160 0.346	0.305	0.869	0.174	0.233	0.808
	∆x =1/320 0.234	0.245	0.876	0.109	0.184	0.812
	∆x =1/640 0.157	0.197	0.882	0.0683 0.145	0.814
	order ≈	0.57	0.31	0	0.67	0.34	0

  p∂ x + q∂ y + r∂ z )u -∆x 2 (p∂ xx + q∂ yy + r∂ zz ) u + ∆t 2 (p∂ x + q∂ y + r∂ z ) 2 u = 0.It can be rewritten as∂ t u + (p∂ x + q∂ y + r∂ z )u -∆x 2 (1 -ν)(p∂ x + q∂ y + r∂ z ) 2 u

	-	∆x 2	pq(∂ x -∂ y ) 2 u -	∆x 2	pr(∂ x -∂ z ) 2 u -	∆x 2	qr(∂

y -∂ z ) 2 u = 0.

Assume the numerical fluxes satisfies two family of inequalities which are, firstly the incoming conditions (44) for all i, j, and secondly the following outgoing conditions for all i, j U i,j + 1-ν ν (U i,j -M i,j ) ≤ pU i+ 1 2 ,j + qU i,j+ 1 2 , U i,j + 1-ν ν (U i,j -m i,j ) ≥ pU i+ 1 2 ,j + qU i,j+ 1 2 .

(46)

Inequalities (44) are called incoming because they concern incoming fluxes U i-1 2 ,j and U i,j-1 2 in a given cell with indices i and j. Inequalities (46) are called outgoing because they concern outgoing fluxes U i+ 1 2 ,j and U i,j+ 1 2 . The first inequalities can also be interpreted as a kind of consistency requirement of the incoming flux. Proof of Lemma 15. From (45), ( 44) and ( 46), one obtains the upper bound

that is U i,j ≤ M i,j after simplifications. Similarly one has

that is U i,j ≥ m i,j after simplifications.

Lemma 16. Assume 0 < ν ≤ 1. Then the upwind fluxes satisfy all inequalities (44) and (46).

Proof. The bounds (44) are already verified in Lemma 14. The remaining inequalities (46) are evident because the upwind fluxes yield pU i+ 1 2 ,j + qU i,j+ 1 2 = U i,j . 5.1.2. Explicit formula for ϕ i,j

Next we consider the scheme (34) defined by the minmod type limiters

where

and by

The latter means that the correction terms at corners are not taken into account. They will be reintroduced later.

Lemma 17. The scheme (34-35) with (47-49) satisfies the special form (41) of the maximum principle, is LVD and is stable in quadratic norm. It is formally only first order accurate.

Proof. By definition, the requirements (36) are satisfied so quadratic stability holds. The scheme is only firstorder accurate because of (49), that is the corner contributions miss with respect to the formulation (31-33) of the second-order accurate Lax-Wendroff scheme.

To show that it is LVD and satisfies the maximum principle under the form (41), one relies on the verification of of (44-46). One has

so the double inequality (44) holds for the total incoming flux. Moreover one has by construction (47-48) that

So one can write

.

Then one has the lower bound

from which the first inequality of (46) is deduced. Similar manipulations yield the second inequality of (46). Therefore the scheme is LVD.

Next we desire to show the asymptotic value ϕ i,j ≈ 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 ϕ i+ 1 2 ,j+ 1 2 ≈ 1, contrary to (49) which is analyzed in this section. We consider the smooth solution u issued from an initial data (2) and the discrete data obtained by interpolation U i,j = u(i∆x, j∆x) for all (i, j).

(50)

We will say that the solution is locally non characteristic if

This is a 2D extension of a point which is not an extremum in the classical theory of limiters [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics-A Practical Introduction[END_REF].

Lemma 18. Consider the scheme (34-35) with the limiters (47-49) and with the initialization (2-50). Assume the solution is locally non characteristic. Then ϕ i,j = 1 + O(∆x).

Proof. We note x i = i∆x and y j = j∆x. From (48), a Taylor expansion shows U i,j -pU i-1,j -qU i,j-1 = ∆x (p∂ x u + q∂ y u) (x i , y j ) + O(∆x 2 ).

Of course one also has U i-1,j -pU i-2,j -qU i-1,j-1 = ∆x (p∂ x u + q∂ y u) (x i , y j )+O(∆x 2 ). The non characteristic condition yields that r i-1 2 ,j = 1 + O(∆x). Similarly r i,j-1 2 = 1 + O(∆x). The claim is obtained by (47).


Saikou Kolley

Keywords: Detection and Prevention Phishing attacks: Detection and prevention Phishing, Detection, Prevention, Cyber Criminals, Email, URL, Domain Name, Fraud, Online, Internet

come

I. INTRODUCTION

The advent of the internet in the early 70s and 80s has brought about enormous advancement in terms of technology, knowledge sharing and material and financial prosperity. Today online activities dominate almost every aspect of the society. These activities range from education, banking, research, sports, shopping, entertainment, social interactions to Zoom meetings [START_REF] Graham | Protect.Comply.Thrive IT Governance Blog[END_REF].

In the vast majority of cases the internet is mostly safe and secure. But the success of the internet is not lost on the criminal fraternity. Criminals have found ways and means to disrupt the productivity, safety, and order of the internet. They do this by using their technical knowhow to attack or defraud organizations businesses and individuals for financial, political or other kinds of personal gains [2,[START_REF] Brinkmann | Scan files for hidden executable content[END_REF]. As a result of these attacks, the internet has now become a minefield to be careful of when using it. A simple seemingly helpful email from the HMRC could turn an individual's life into a nightmare. One of the most popular tools cyber criminals use to carry out their fraudulent activitiesamong other tools -is phishing [START_REF] Ghafir | Malicious SSL Certificate Detection: A Step Towards Advanced Persistent Threat Defence[END_REF].

This report aims to demonstrate and explain what phishing is, how it is carried out, its financial cost, its impact on society and how to detect and mitigate against it

II. WHAT IS PHISHING?

Phishing is the attempt to steal personal, confidential and sensitive informationusually through electronic communicationsby fraudulent means. Phishing attacks are usually carried out by using fake emails, fake URL links, popups or fake websites. It is better for a cyber criminal to make a victim click a malicious link in an email or a spoofed website than actually attempt to break into a computer system with robust defences. That is why phishing has become the number one choice for cyber crime. Malicious emails, URLs and websites can easily be created by a criminal method called spoofing [5]. Spoofing is the method of disguising emails, websites or links in an attempt to hide their true origin and or intent. Phishing falls under a category of attacks called advanced persistent threat (APT) [START_REF] Buber | Phishing URL Detection with ML[END_REF]. These are well organized attacks mostly motivated by data theft or financial gain and mostly directed towards individuals, businesses, and politicians.

III. COMMON TYPES OF PHISHING ATTACKS

Hackers sometimes engage in social engineering in order to get information. Social engineering is the attempt of manipulating people in order to gain confidential information [START_REF] Raza | An IoT and Business Processes Based Approach for the Monitoring and Control of High Value-Added Manufacturing Processes[END_REF]. Phishing is a social engineering vector. There are several ways by which a phishing attack is carried out: Email Phishing: This is where a fraudulent person registers a fake domain name that looks very similar to a real organisation's domain. The cyber criminal then uses the domain to spoof genuine email addresses which would then be used to fraudulently send random emails to unsuspecting victims. To disguise the emails several techniques can be used. One of these techniques is character substitution. This is where characters in the domain name or URL are substituted to mimic a genuine URL or domain name. Example an 'r' and 'n' written together 'rn' to mimic the letter 'm'. In this example 'rnicrosoft' could easily be mistaken for 'microsoft' [START_REF] Ghafir | A Survey on Intrusion Detection and Prevention Systems[END_REF]. The Cyber criminal can then create a fake 'rnicrosoft' account from which they can send fraudulent emails to unsuspecting victims.

Another technique is the use of the local part of a genuine email address and a fake domain, example 'microsoft@refundme.com'. In this example the local part (Microsoft) is genuine but the domain (refundme.com) is not. A victim seeing 'microsoft' in the local part of the email address might never suspect that the email was not from the genuine Microsoft [START_REF]:1611402804108&source=lnms&tbm=isch&sa=X&ve d=2ahUKEwig--md_7HuAhWRi1wKHVbbA4cQ_AUoAXoECBwQAw &biw=1536&bih=722#imgrc=4O3KuvYOC8XWnM[END_REF]. In order to carry out a successful email phishing campaign, cyber criminals need a lot of email addresses and sometimes they do not have to social engineer to get them. There are a few online tools that help them harvest emails. One of them is Jigsaw, the other is theHarvester. theHarvester is a python script that pulls out email addresses from major search engines [START_REF] Ghafir | Social Engineering Attack Strategies and Defence Approaches[END_REF]. With these tools, a cyber criminal can harvest as many emails as they like, clone them and use them for phishing and other malicious purposes and because they would look as if they are genuine, victims will find it difficult to figure out the messages are actually fake and the codes in them are malicious.

Vulnerabilities of Domains:

The vulnerabilities of domain names and URLs are also fully exploited in email phishing [START_REF] Hmrc | Examples of HMRC related phishing emails, suspecious phone calls and texts[END_REF]. A URL typically consists of a protocol, a subdomain, a root domain, a top domain, a subdirectory, and a file name. The combination of the subdomain, the root domain and the top level domain makes up the domain name. Domain names are registered and therefore constrained to the registered owner. The subdomain and the URL path however have a vulnerability because they can be manipulated by a cyber criminal and used for phishing [START_REF] Ghafir | Detection of Advanced Persistent Threat Using Machine-Learning Correlation Analysis[END_REF]. The vulnerability of domains makes it possible for hacker to manipulate and clone them and there are many websites that enable anyone to own and operate a domain name. An example of such a website is GoDaddy [START_REF] Irwin | IT Governance European Blog[END_REF].

Cyber criminals can register fraudulent domain names and email addresses that look like legitimate ones at GoDaddy.com. Criminals can also clone websites using a tool in Kali called Settoolkit. Those websites can then be used for criminal purposes including phishing. IV. REDIRECTS Cyber criminals do not only have to clone fake emails or websites, they can also take advantage of the vulnerability in HTTP request smuggling where some applications allow redirects from one URL to another [START_REF] Kaspersky | What is a whaling attack[END_REF]. They could use this technique to send you an email with a link that redirects a victim to a fake destination.

For example, a request is made for a file that does not have a trailing slash, but gets redirected to the same file but this time the slash is included [START_REF] Ghafir | Malicious File Hash Detection and Drive-by Download Attacks[END_REF].

GET /home HTTP/ 1.1 Host: normal-website.com HTTP/ 1.1 301 Moved Permanently Location: HTTPS: // norman-website.com/home/ This is harmless under normal circumstances but a cyber criminal can exploit it to redirect a victim to a different domain as seen in the example above.

The code below could also be used for redirects. POST / HTTP /1.1 Host: vulnerable-website.com Content-Length: 54 Transfer-Encoding: chunked 0 GET /home HTTP/1.1 Host: attacker-website.com Foo: x This will redirect a victim to the attacker's website and will affect subsequent requests processed by the back-end servers: GET: /home HTTP/1.1 Host: attacker-website.com Foo: XGET / scripts/include.js HTTP/1.1 Host: vulnerable-website.com HTTP/1.1 301 Moved Permanently Location: https//attacker-website.com/home/ (Portswigger, 2021).

V. SPEAR PHISHING

When phishing is targeted towards an individual, it is called spear phishing. For this to work, a criminal has to have personal detailsname, email, address, job title or similar personal informationof the victim.

VI. WHALE PHISHING

Whale phishing is targeted towards the wealthy, people with power or senior executives of an organization. In this method a cyber criminal would masquerade as someone with authority and sends fraudulent email, letter or memo to a target victim who complies with a demand [START_REF] Pentestgeek | How do I phish -Advanced Email Phishng Tactics[END_REF][START_REF]Anti-phishing toolbar: A tool to combat phishing attacks[END_REF]. An example of this is, in 2016 Snapchat became a whaling attack victim when their payroll department received an email that supposedly came from the CEO asking for employee payroll information [START_REF] Ghafir | Security Threats to Critical Infrastructure: The Human Factor[END_REF].

VII. CLONE PHISHING

This is where cyber criminals obtain legitimate emailswith links and attachmentsthat have been delivered in the past. The address of the sender and recipient will then be used to create a clone of the original but this time, the content, the links and attachments would contain malicious code. The emails would then be resent using spoofed original email addresses. The unsuspecting recipient will receive the email as a resend of the original email and may not suspect or question its authenticity, thus falling victim to scam [START_REF] Eltanani | Coverage Optimisation for Aerial Wireless Networks[END_REF].

VIII. VOICE PHISHING

Organizations, businesses or individuals do not have to receive a fraudulent email, visit a fake website or click a malicious link to become a phishing victim. They can easily and equally become a scamming victim by a simple telephone Call [START_REF] Portswigger | Exploiting HTTP request smuggling vulnerabilities[END_REF]. Voice phishing is where criminals directly call organizations, businesses or individuals pretending to represent some legitimate authority. They would then dupe them into handing over personal or financial information. Voice phishing is extremely common and usually targets the elderly [START_REF] Ghafir | Disguised Executable Files in Spear-Phishing Emails: Detecting the Point of Entry in Advanced Persistent Threat[END_REF].

IX. SYMPTOMS OF COMPROMISED SYSTEMS

Apart from the obvious demand for ransom and the presence of fraudulent links and emails on a system, the following can be symptoms of a compromised system.

•

Slow network devices: If you noticed an unexplained slow down in network performance, that could be a sign that your system has been compromised.

•

Presence of programs not installed by the organization: If programs that were not previously installed begin appearing on your system, that could be a sign of a comprised system. • Random pop-ups: Random unsolicited pop-ups are a sign of an infected system. • Redirects: Redirects are when a web page changes URLs and redirects a user to a totally different location when clicked. This is a good sign that it might be a malicious link.

X. EXAMPLES OF PHISHING ATTACKS

The following is an HMRC email scam that claims that customers can ask for a refund to help mitigate against the negative financial effects of Covid-19. They can do so by clicking the link in the email:

XII. SINGATURE BASED DETECTION

Cyber criminals use legitimate looking URLs, emails, and websites to carryout their malicious activities. One of the ways to detect this is by using digital signatures. The signature-based detection method analyses files and assigns a hash to them. These signatures are then stored in databases for future comparison with known malware. All files are then checked against theses databases to detect malicious files. Signature based detection methods are very useful tools but sometimes they are not good enough, so new detection methods have to be invented including the following:

XIII. HEURISTIC ANALYSIS Heuristic Analysis: Heuristic analysis first determines patterns for normal software and system behavior and then flogs up anything that is out of order as suspicious [START_REF] Ghafir | Proposed Approach for Targeted Attacks Detection[END_REF].

XIV. SANDBOXING Sandboxing: In this method, suspicious files are isolated and analyzed for malicious code before being allowed to run in the system [START_REF]pHEKHddXAaIQ_A UoAXoECBEQAw&biw=1536&bih=722#imgrc=dNXz iaZB2Zr5hM[END_REF].

XV. AL-BASED ANALYSIS

This is where machine artificial intelligence is used to work out what is malicious code and what is not. This involves machine learning where AI is trained to recognize malicious code.

XVI. DOMIAN-BASED DETECTION

A bit of detective work can throw up a lot about a registered domain. Sitechecker can be used to find out about the credentials of a site. Some of the things to look out for are: • Was the domain blacklisted.

•

Is the registrant hidden [START_REF] Ghafir | BotDet: A System for Real Time Botnet Command and Control Traffic Detection[END_REF].

XVII. FILE EXTENSION BASED DETECTION

File extensions can be faked. The Unitrix exploit uses special Unicode characters to reverse the true name and sequence of a file enabling it to hide the extension in the middle. Unitrix can reverse [u+202E]3PM.SCRwhich is the malicious file nameinto RCS.MP3. In this example, the extension is not a harmless MP3 file but a malicious SCR file masquerading as a harmless music file [START_REF] Hoffman | How hackers can hide malicious programs with a fake malicious extension[END_REF].

To detect some of these malicious extensions, use the MZReveal software. MZReveal can help detect malicious file extensions. It does this by detecting and logging hidden files that are executable [START_REF] Svoboda | Network Monitoring Approaches: An Overview[END_REF].

XVIII. URL BASED DETECTION

One of the ways to detect whether a link is what it says it is, is to hover over the link. This would tell the real origin of the link.

XIX. ANTI PHISHING TOOLBARS

Anti-phishing toolbars are not perfect but they do help detect and stop some malware. Some examples are:

• Netcraft: Netcraft is a toolbar with anti-phishing capabilities, it shows registration dates, locations of domains and other vital information about the reputation of a domain [START_REF]Social engineering (security)[END_REF].

• SpoofGuard: SpoofGuard is another toolbar that helps detect phishing attacks. It works by calculating a spoof score of a visited page using heuristics from previous phishing attacks [START_REF]Phishing[END_REF].

XX. PREVENTION

Phishing scams are increasing year-on -year because cyber criminals are making huge profits. But scams can be prevented by common sense, online awareness and following some simple basic rules [START_REF] Ghafir | Advanced Persistent Threat and Spear Phishing Emails[END_REF]. This section will be divided into two. Section one will discuss what the ordinary online user can do to prevent phishing and section two will discuss what an IT manager can do.

XXI. ORDINARY ONLINE USER

Ordinary online users must be vigilant at all times and should not fall for all the things they came across at face value, so keeping to the following rules will help.

• Too good to be true emails: Online users should be wary of emails that promise refunds, discounts or account updates. If in doubt call the senders of the email directly to verify their claims.

•

URLs: URLs can be easily disguised so think carefully before clicking on a URL. A good practice is to hover over the URL, if the result matches what the sender claims it to be, then it is a genuine URL. • Browsers: More often than not, browsers can provide minimal protection from cyber criminals, so users should always make sure that browsers are up to date.

•

Pop-ups: Pop-ups should be avoided at all cost if possible because they are one to the ways malware is propagated. So browsers should be configured to block popups.

• Firewalls: Firewalls play a great role in blocking malware. So firewalls should always be installed, up to date and activated. • Antivirus Software: Using antivirus software is the best way to block malware so investing in some good antivirus software is highly recommended. • Anti-Phishing Toolbars: Anti-phishing toolbars can detect block and report malware in some cases. So installing one is a good idea.

XXII. NETWROK IT MANAGER

Organizations and businesses are every-day targets for cyber criminals, that is why IT / network managers should be up to date with trending cyber attacks and advise their staff accordingly. The following are some of the steps that need to be taken to protect their businesses and organizations: Phishing: £960,000 and took 20 days to resolve. Ransomware: £520,000 and 23.1 days to resolve. Botnets: £260,000 and 2.5 days to resolve [START_REF] Graham | Protect.Comply.Thrive IT Governance Blog[END_REF].

XXIV. CONCLUSION

This report started by introducing the advent of the internet and its good intentions then highlighted the fact that the internet will be used for good and evil. It went on to describe phishing which is the main topic of the report.

The report then went on to discuss the different types of phishingemail phishing, spear phishing, whale phishing, clone phishing and voice phishingand the methods used in these attacks. It then discussed the symptoms of a compromised system other examples of phishing scams. The report concluded by discussing how phishing can be detected and prevented and the social and financial cost of phishing.

Figure 1 :

 1 Figure 1: Advanced persistent threat life cycle (Concept draw, 2020).

Figure 2 :

 2 Figure 2: Python theHarvester script (Pentestgeek, n.d.).

Figure 3 :

 3 Figure 3: URL structure (Cult of web, n.d.).

Figure 4 :

 4 Figure 4: GoDaddy.com (Pentestgeek, n.d.).

Figure 5 :

 5 Figure 5: Setoolkit in Kali.

Figure 6 :

 6 Figure 6: HMRC email scam (HMRC, 2021).

Figure 7 :

 7 Figure 7: HMRC text scam (HMRC, 2021).

  Apart from the embarrassment of being the victim of a phishing attack, there is a huge financial cost as well. But that is not all, the reputational cost cyber criminals can inflict on a business or organization is immense. No one would like to trade or do business with a company with lax network security. Lloyd's of London has said that the cost of cyber attacks is in the region of $120 billion[START_REF] Graham | Protect.Comply.Thrive IT Governance Blog[END_REF]. The following figures show the cost and the number of days it took to resolve an incident: Malware: £1.57 million and took 6.4 days to resolve. Web-based attacks: £1.52 million and took 22.4 days to resolve. DOS attacks: £1.31 million and took 16.8 days to resolve. Malicious code: £960,000 and took 55.2 days to resolve.

	•	Spam Filters: Install spam filters and always keep
	them up to date.
	•	Web Filters: Install a web filter to weed out all
	known malicious websites.
	•	Antivirus Software: Deploy good antivirus software
	for maximum protection.
		XXIII.	SOCIAL COST OF PHISHING
			•	Education and Training: IT managers should
			conduct regular education and training to sensitise staff about
			phishing, what to look out for and how to prevent the attacks.
			•	Security policy: A good security policy will go a
			long way in protecting businesses and organizations form
			cyber attacks. It must include password expiry dates and two-
			time authentication among other things.