Andrus Giraldo

Bernd Krauskopf

Neil G R Broderick

Juan A Levenson

Alejandro M Yacomotti

The driven-dissipative Bose-Hubbard dimer: phase diagram and chaos

We present the phase diagram of the mean-field driven-dissipative Bose-Hubbard dimer model. For a dimer with repulsive on-site interactions (U > 0) and coherent driving we prove that Z2symmetry breaking, via pitchfork bifurcations with sizable extensions of the asymmetric solutions, require a negative tunneling parameter (J < 0). In addition, we show that the model exhibits deterministic dissipative chaos. The chaotic attractor emerges from a Shilnikov mechanism of a periodic orbit born in a Hopf bifurcation and, depending on its symmetry properties, it is either localized or not.

The Bose-Hubbard model is a celebrated fundamental quantum mechanical model that accounts for boson dynamics in a lattice [1]. It successfully describes the interplay between the hopping of particles between neighboring sites of the lattice (with rate J) and on-site interactions. Such interactions appear as multi-boson terms in the Hamiltonian with interaction energy U . Importantly, this model accurately explains the superfluid to Mott insulator phase transitions, that has been experimentally demonstrated in ultracold atomic lattices [2]. A minimal building block in this context is the so-called Bose Hubbard dimer, consisting of only two interacting sites, also known as the bosonic Josephson junction [3,4]. Futhermore, the Bose-Hubbard dimer lies at the basis of a number of striking phenomena such as the Josephson effect, self-trapping and symmetry breaking [5][6][7].

In recent years there has been a growing interest in understanding open quantum systems, where the bosons can be added and destroyed by means of external driving and dissipation mechanisms [8][9][10][11]. In this context, photonic systems have attracted much attention since photons in optical cavities can be injected through an external driving laser, and dissipation comes in as a natural consequence of optical cavity losses. The drivendissipative Bose-Hubbard dimer has been realized in a number of experimental systems, such as semiconductor microcavities [12,13], superconducting circuits [14] and photonic crystals [15], where the interactions take the form of Kerr-type optical non-linearities. In many regards, these optical systems constitute outstanding platforms for studying many-body phenomena in open quantum systems [16]. Among them, dissipative phase transitions [17,18] are an especially exciting open topic, because they provide a conceptual basis for the understanding and the prediction of new collective states, both steady and dynamical ones, with the latter accounting for collective coherent oscillations.

As is well known, phase transitions are characterized by critical phenomena, which can emerge in the thermodynamic limit, i.e., with a large photon number in optical cavities. Remarkably, even single-mode cavities -i.e., with no spatial degrees of freedom -may display phase transitions including optical bistability (which is of first order) [18,19] and the emergence of an oscillation threshold in, e.g., two-photon pumped Kerr resonators [20] or laser devices [21] (which is a second-order phase transition). As pointed out in Ref. [20], driven-dissipative phase transitions are strongly related to mean-field semiclassical solutions in such a thermodynamic limit: these are known as bifurcations of a classical vector field from a non-linear dynamical point of view; see, e.g., [START_REF] Guckenheimer | Nonlinear oscillations, dynamical systems, and bifurcation of vector fields[END_REF][START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF].

Recent examples of complex dynamics found in phase diagrams of driven-dissipative nonlinear cavities include oscillating phases [START_REF] Jin | [END_REF] and exotic attractors [25]. Moreover, Lorenz-type chaos has already been found in the Gross-Pitaevsky equation accounting for incoherent pump of a double potential well Bose-Einstein condensate [26]. However, to our knowledge, dissipative chaos has not been predicted so far in the coherently driven regime, which is the one that is well described by the driven-dissipative Bose-Hubbard dimer. In this work, we show that dissipative chaos exist in this model and that it is intimately related to Shilnikov homoclinic bifurcations. In addition, we present the phase diagram as a result of the comprehensive analysis of local bifurcations of steady states and their symmetry properties.

The mean-field approximation of the Bose-Hubbard model is (see, e.g., [27])

i dα 1 dτ = -∆ -i γ 2 + 2U |α 1 | 2 α 1 -Jα 2 + F, i dα 2 dτ = -∆ -i γ 2 + 2U |α 2 | 2 α 2 -Jα 1 + F, (1)
where α 1 , α 2 ∈ C are the electric field envelopes in each optical cavity, which are linearly coupled by the tunneling parameter J. Both cavities are coherently driven by a field with amplitude F and detuning ∆ = ω p -ω c , where ω p and ω c are the driving and cavity frequencies, respectively. The cavities have a Kerr-type nonlinearity characterised by U , which is assumed to be fast compared to the dissipation rate γ. We apply to system (1) the coordinate transformation

(α 1 , α 2 ) → (A, B) = (-2iα * 1 |U |/γ, -2iα * 2 |U |/γ),
where α * i is the complex conjugate of α i . With time rescaled to τ = γt/2, we obtain the non-dimensionalized Bose-Hubbard model as the vector field

dA dt = -A + i δ + sign(U )|A| 2 A + κB + f, dB dt = -B + i δ + sign(U )|B| 2 B + κA + f, (2)
for the (rescaled) envelopes A, B ∈ C. Here

κ = - 2J γ , f = 2F |U | γ 3/2 and δ = - 2∆ γ
are the (rescaled) coupling, drive and detuning parameters, respectively.

Since system (2) is invariant under the transformation

(A, B, U, δ, κ) → (A * , B * , -U, -δ, -κ),
all results for positive U directly translate to those for negative U ; hence, we consider U > 0 and sign(U ) = 1 in (2) from now. Note that, although U is fixed by the material of the cavity, both the magnitude and sign of κ can be changed by means of photonic design in some particular geometries, for example, with potential barrier engineering for the case of photonic crystal nanocavities [15,28]. Because δ can be altered at will during the experiment by changing the frequency of the driving laser, states of negative interaction energy in an otherwise positive-U system can be assessed via the parameter transformations κ → -κ and δ → -δ. In addition, system (2) is invariant under the phasespace mirror symmetry of exchanging the two cavities, that is, of swapping A and B, which gives rise to Z 2equivariance of solutions [START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF]. Hence, solutions (equilibria, periodic orbits, trajectories) can be split into two major groups: symmetric and asymmetric ones. For a symmetric solution the field intensities and phases in both cavities are the same. The set of all symmetric solutions forms the symmetry subspace given by A = B, which is an invariant subspace of system (2). An asymmetric solution, on the other hand, is one where the intensity and/or the phase in both cavities differ, that is, A = B. Note that asymmetric solutions often come in pairs, one being the mirror image of the other under swapping A and B.

We start our analysis of system (2) by characterizing the equilibria and their bifurcations as the drive parameter f increases. against f for three different values of δ for the two cases of a negative and positive coupling parameter κ in the left and right columns, respectively. The corresponding phase diagrams of |B| 2 against f are the same due to Z 2 -equivariance. All bifurcations and branches of solutions have been computed with the numerical continuation package Auto07p [START_REF] Doedel | AUTO-07p: Continuation and Bifurcation Software for Ordinary Differential Equations[END_REF]. When f = 0 there is the stable equilibrium given by A = B = 0, which gives rise to a monotone branch of stable symmetric equilibria for any f provided the detuning δ is sufficiently large. However, when δ = 0 and for increasingly negative δ, one finds bifurcations on the branch of symmetric equilibria.

For negative κ, as in Fig. 1a-c for κ = -3.5, we observe for δ = 0 an interval of bistability that is delimited by two saddle-node bifurcations S of symmetric equilibria [Fig. 1a]. In this interval there are three symmetric equilibria that form a hysteresis loop: a stable equilibrium corresponding to low intensity in both cavities, an intermediate saddle equilibrium, and a stable equilibrium with higher intensity in both cavities. As δ decreases, the interval of bistability increases in size. While bistability can be encountered in one-cavity systems, spontaneous Z 2 -symmetry breaking requires coupled cavities with mirror symmetry; this is the case for the Bose Hubbard model ( 2) and spontaneous symmetry breaking phase transitions are known to exist in the form of pitchfork bifurcations [27,[START_REF] Cao | [END_REF]. As Fig. 1b for δ = -6 shows, we find two such symmetry breakings at the points labelled P on the low intensity branch. The left pitchfork bifurcation point P is subcritical and gives rise to a pair of unstable asymmetric equilibria, while the right point P is supercritical and gives rise to a pair of stable asymmetric equilibria (one with |A| 2 < |B| 2 and the other with |A| 2 > |B| 2 ); these two branches come together at a pair of saddle-node bifurcations of asymmetric equilibria denoted S * . These additional bifurcations create a small f -interval, between S * and the left point P, with four stable equilibria: higher and lower intensity symmetric equilibria, and a pair of asymmetric (lowerintensity) equilibria [Fig. 1b]. Notice that between the two points P there are still three stable objects: the two asymmetric equilibria and the higher-intensity symmetric equilibrium. For an even lower detuning, as in Fig. 1c for δ = -6.2, we find two pairs of supercritical Hopf bifurcations denoted H on the stable branch of asymmetric equilibria. They give rise to a pair of stable asymmetric periodic orbits, which are represented in the phase diagram by their maxima in |A| 2 . In the corresponding f -interval the asymmetric equilibria are now unstable. Note that the two asymmetric periodic orbits are new attractors that coexist with the stable higher intensity equilibrium.

The situation for κ > 0 is notably different: bistability of the symmetric state is, in fact, absent for all panels Fig. 1d-f. Instead, for δ = 0 we find symmetry breaking at two supercritical pitchfork bifurcations, which give rise to a pair of stable asymmetric equilibria in an quite large f -interval where the symmetric equilibrium is unstable [Fig. 1d]. This makes this set of abnormal coupling conditions well suited for experimental implementations. Other driving configurations might also allow one to fulfill these conditions, but then the driving phase of the cavities needs to be adjusted properly such that the antisymmetric state can be linearly excited [27].

As δ is decreased, as in Fig. 1e for δ = -3.0, the pair of asymmetric stable equilibria that exists in this interval of symmetry breaking also exhibit saddle-node bifurcations S * that create bistability of asymmetric states and a pair of hysteresis loops. Perhaps surprisingly, there is an frange where the intensity in one of the cavities is nearly zero. As Fig. 1f for δ = -4.5 shows, a high negative detuning generates a pair of supercritical Hopf bifurcation points labelled H, which create an f -interval with pair of asymmetrical stable periodic solutions, while the asymmetric equilibria are unstable.

The bifurcations of equilibria identified in Fig. 1 can be followed in the additional parameter δ, yielding a phase diagram in the (f, δ)-parameter plane consisting of bifurcation curves that bound regions with qualitatively different behavior. Figure 2 shows the two phase diagram for both κ = -3.5 and κ = 3.5 showing the curves S, S * , P and H of saddle-node, pitchfork and Hopf bifurcations, respectively. The dashed horizontal lines correspond to the respective phase diagrams in just f from Fig. 1a-f alytical expressions obtained by solving for the ria and the respective bifurcation condition [31]. These calculations also show that the saddle-node bifurcation S in the symmetry subspace (brown line) only occurs if δ < -√ 3 -κ and 1.24081 < f . At the points (f, δ) ≈ 1.24081, -√ 3 -κ there is a cusp bifurcation CP on the curve S that corresponds to a sharp bound for bistability in the (f, δ)-plane [START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF]. The region of bistable symmetric equilibria bounded by S is shaded in brown in both panels of Fig. 2; note the quite different position of this region for κ = -3.5 and for κ = 3.5. Additionally, we find that the curve of pitchfork bifurcation P (purple line), which bounds the shaded region with asymmetric equilibria, has a maximum at δ = -√ 3 + κ; hence, symmetry breaking only occurs if δ < -√ 3 + κ, which is equivalent to the bound found in [START_REF] Brunstein | Nonlinear Dynamics in III-V Semiconductor Photonic Crystal Nano-cavities[END_REF] and also consistent with [START_REF] Maes | [END_REF]. As an unexpected consequence, the condition for symmetry breaking at P to be possible is that the detuning ω p -(ω c -κγ/2), with respect to the antisymmetric mode, must exceed √ 3/2 times the cavity linewidth. This is similar to the bistability condition, but the detuning is measured from the antisymmetric mode of the system, i.e., the one that cannot be linearly excited from the outside world. These bounds explain why the first phenomenon observed, as δ decreased, is bistability when κ < 0, while it is symmetry breaking when κ > 0; compare the two panels of Fig. 2 and see also Fig. 1. Notice that the curve P is tangent to the curve S at the saddle-node pitchfork point SP, to the right of the cusp point CP for κ < 0 and to the left of CP for κ > 0; moreover, for κ < 0 the curve P changes criticality at the point DP, from which the curve S * of saddle-node bifurcations of asymmetric equilibria emerges.

The curves SN * and H in Fig. 2 of saddle-node and Hopf bifurcations or asymmetric equilibria are computed numerically by continuation techniques. They bound regions of bistability and of stable asymmetric periodic solutions, respectively. Notice that for κ < 0 the curve H emanates from a saddle-node pitchfork point SP and is superctitical throughout, while for κ > 0 we find a change of criticality of H at a generalised Hopf bifurcation (labelled GH) [START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF]. The considerable amount of multistability of the Bose-Hubbard model ( 2) is represented by the overlap between the different shaded region of Fig. 2. For even lower negative values of δ than shown in Fig. 2, one can find up to nine equilibria, as pointed out in [START_REF] Cao | [END_REF].

As δ is decreased for both signs of κ, the asymmetric stable periodic orbits created at the supercritical Hopf bifurcation H exhibit bifurcations scenarios to more complex and chaotic dynamics; the region of complex dynamics is indicated by grey shading in Fig. 2. Its exact boundary is formed by an intricate arrangement of accumulating bifurcation curves, period-doubling and different types of homoclinic and heteroclinic bifurcations [31]. Notice that for κ < 0 periodic solutions become more complex almost immediately, while for κ > 0 we find a quite large region of stable periodic solutions. While a full discussion of the scenarios is beyond the scope of this paper, we conclude by showing the existence of two types of Shilnikov-type chaos in the Bose-Hubbard model: asymmetric chaotic dynamics localized in one of the cavities, and symmetric chaotic dynamics with irregular switching between the two cavities.

At a homoclinic bifurcation a system of differential equations posesses a special trajectory, called a homoclinic orbit, that converges both in backward and forward time to a saddle equilibrium. For the special case that the equilibrium is a saddle focus (with a pair of complex conjugate eigenvalues), a celebrated result by Shilnikov [34] states that (under a certain eigenvalue conditions) infinitely many periodic orbits and chaotic behavior can be found nearby. This phenomenon is now referred to as a Shilnikov bifurcation. In the Z 2 -equivariant system (2) we find Shilnikov bifurcations both to asymmetric saddle foci and to symmetric saddle focus, as well as two types of nearby chaotic attractors.

Figure 3 shows both cases of an asymmetric and of a symmetric Shilnikov bifurcation and associated chaotic attractors. For κ < 0 in Fig. 3a, there exist asymmetric Shilnikov homoclinic orbits to each of the pair of asym- metric equilibria, one either side of the symmetry line |A| 2 = |B| 2 . For nearby f as in Fig. 3b, we find a pair of chaotic asymmetric and localized attractors. The time traces in panels (c1,c2) show that this type of localized chaotic dynamics is characterized by a dominance in intensity of one of the two cavities for all time. For κ > 0, on the other hand, we find a pair of Shilnikov homoclinic orbits to a single symmetric saddle equilibrium [Fig. 3d]. This has the consequence that the associated chaotic attractor in Fig. 3e is much larger and no longer localized as it crosses the line |A| 2 = |B| 2 . As the time traces in panels (f1,f2) clearly show, this type of chaotic dynamics is characterized by episodes of dominance in intensity of one of the two cavities for some time, with rapid and irregular switching events of the intensity to the other cavity. We remark that a local analysis of equilibria alone is not able to predict and explain this global switching behavior of the system.

In conclusion, our systematic bifurcation analysis enabled a complete phase diagram description of both equilibria and time-dependent solutions of the drivendissipative Bose-Hubbard dimer model. In particular we were able to identify bifurcations to limit cycles, eventually leading to localized and non-localized chaotic behavior near homoclinic bifurcations of Shilnikov type. Crucial to the dynamics is the fact that the two cavities are coherently driven, which means that the system remains four-dimensional since there is no phase freedom as in ac-tive cavities. Interestingly, in a system with positive onsite interaction energy (U > 0), nonlocalized Shilnikov chaos can be observed for negative tunneling parameter (J < 0), in a bifurcation cascade whose details will be reported elsewhere. Our mean-field dynamical results thus pave the way to study complex nonequilibrium orbits -including quantum dissipative chaos-in the quantum master equation of the open Bose-Hubbard dimer.

This work has been partially funded by the "Investissements d'Avenir" program (Labex NanoSaclay, Grant No. ANR-10-LABX-0035) and the ANR UNIQ DS078.

FIG. 1 .

 1 FIG.1. Phase diagram of system (2) of intensities |A| 2 against the drive parameter f at different δ values with κ = -3.5 (left column) and κ = 3.5 (right column); shown are branches of stable equilibria (blue), of saddle equilibria with one unstable eigenvalue (cyan) and with two unstable eigenvalues (orange), and of stable periodic solutions (dark green).

FIG. 2 .

 2 FIG. 2. Phase diagrams of system(2) in the (f, δ)-plane for κ = -3.5 (top panel) and κ = 3.5 (bottom panel). Shown are the curve of saddle-node bifurcations of symmetry equilibria S (brown) and asymmetric equilibria S * (orange), of pitchfork bifurcations P (purple), and of Hopf bifurcations H (green); also shown are points CP and CP * of cusp, SP of saddle-node pitchfork, DP of degenerate pitchfork, and GH of generalised Hopf bifurcation. Regions are shaded: brown for symmetric bistability; purple for symmetry broken equilibria; orange, yellow and red for multistability of equilibria; green for coexisting stable periodic solutions; and gray for more complex dynamics.

FIG. 3 .

 3 FIG. 3. Shilnikov orbits (a,d), trajectories on chaotic attractors (b,e) in the (|A| 2 , |B| 2 )-plane, and temporal traces of |A| 2 (c1,f1) and of |B| 2 (c2,f2) for the black trajectories. Here κ = -3.5, δ = -8 with f ≈ 6.943 for (a) and f = 6.93 for (b,c); and κ = 3.5, δ = -8 with f ≈ 2.379 for (d) and f = 3.0 for (e,f).


Thierry Goudon

Julie Llobell

Sebastian Minjeaud

An asymptotic preserving scheme on staggered grids for the barotropic Euler system in low Mach regimes

We present a new scheme for the simulation of the barotropic Euler equation in low Mach regimes. The method uses two main ingredients. First, the system is treated with a suitable time splitting strategy, directly inspired from [J. Haack, S. Jin, J.-G. Liu, Comm. Comput. Phys., 12 (2012) 955-980], that separates low and fast waves. Second, we adapt a numerical scheme where the discrete densities and velocities are stored on staggered grids, in the spirit of MAC methods, and with numerical fluxes derived form the kinetic approach. We bring out the main properties of the scheme in terms of consistency, stability, and asymptotic behaviour, and we present a series of numerical experiments to validate the method.

MSC. 65M08, 76M12, 76NXX, 35Q31.

Introduction

Euler equations and low Mach regimes

This work is concerned with the numerical resolution of the barotropic Euler system, which describes the evolution of a compressible fluid, in low Mach regimes. The unknowns ρ and u stand respectively for the local density and velocity of the fluid. They both depend on the time and space variables, t 0 and x ∈ Ω ⊂ R N , where Ω is a bounded domain. Working with dimensionless quantities, the evolution of the fluid is governed by the PDE system

     ∂ t ρ + ∇ • ρu = 0, ∂ t ρu + ∇ • ρu ⊗ u + 1 2 ∇ p(ρ) = 0.
(1

)
The system has obviously to be completed with initial and boundary conditions (periodic or vanishing flux for instance). Here and below, we suppose that the pressure law ρ → p(ρ) is the classical power-law p(ρ) = ρ γ with γ > 1. In [START_REF] Alazard | A minicourse on the low Mach number limit[END_REF], stands for the Mach number, that is the ratio of the typical velocity of the fluid to the typical sound speed

( 2 = |u ref | 2 c 2 ref = ρ ref |u ref | 2 γp ref )
. We are interested in the numerical approximations of the solution of system (1) for 0 < 1 varying from very small values to unity.

Low Mach regimes

In order to understand the behaviour for small 's, let us expand the solution as a power series Q = Q (0) + Q (1) + 2 Q (2) +..., with Q = (ρ, p, u), and identify the terms that arise with the same order of magnitude with respect to in [START_REF] Alazard | A minicourse on the low Mach number limit[END_REF]. At the leading order O 1/ 2 in the momentum equation, we obtain ∇p (0) = 0 from what we deduce p (0) (x, t) = p (0) (t) and consequently ρ (0) (x, t) = ρ (0) (t). Similarly, at order O (1/ ) we get ∇p (1) = 0, thus ρ (1) (x, t) = ρ (1) (t) depends only on the time variable too. Given this information, system (1) at order O (1) becomes ∂ t ρ (0) + ρ (0) ∇ • u (0) = 0, ∂ t (ρ (0) u (0) ) + ρ (0) ∇ • (u (0) ⊗ u (0) ) + ∇p (2) = 0.

(

) 2
Integrating the mass equation in [START_REF] Barsukow | A numerical scheme for the compressible low-Mach number regime of ideal fluid dynamics[END_REF] over the domain Ω, yields

Ω ∂ t ρ (0) + ρ (0) ∇ • u (0) dx = |Ω|∂ t ρ (0) + ρ (0) Ω ∇ • u (0) dx = 0.
Applying the Green theorem leads to

∂ t ρ (0) = - 1 |Ω| ρ (0) ∂Ω u (0) • n dσ(x)
where n is the outward-pointing normal vector on ∂Ω and dσ(x) is the Lebesgue measure. Let us restrict to the case of zero-flux boundary conditions, u • n ∂Ω = 0, or the case of periodic boundary conditions, so that in both cases we have ∂Ω u (0) • n dσ(x) = 0. It follows that ∂ t ρ (0) = 0. Going back to the mass conservation equation in [START_REF] Barsukow | A numerical scheme for the compressible low-Mach number regime of ideal fluid dynamics[END_REF] we arrive at ∇ • u (0) (x, t) = 0 and ρ (0) (x, t) = ρ (0) on Ω × [0, ∞).

The momentum equation in [START_REF] Barsukow | A numerical scheme for the compressible low-Mach number regime of ideal fluid dynamics[END_REF] can now be cast as ∂ t u (0) + ∇ • (u (0) ⊗ u (0) ) + 1 ρ (0) ∇p (2) = 0.

In this relation the pressure term p (2) can be seen as the Lagrange multiplier associated to the incompressibility constraint ∇ • u (0) = 0. Applying the divergence operator we are led to ∆p (2) = -ρ (0) ∇ • ∇ • (u (0) ⊗ u (0) ) .

Finally, the mass equation in [START_REF] Alazard | A minicourse on the low Mach number limit[END_REF] at order O ( ) reads ∂ t ρ (1) + ρ (0) ∇ • u (1) = 0.

With a similar reasoning, we deduce that ρ (1) (x, t) = ρ (1) is constant in time and space and that ∇ • u (1) = 0. Definition 1.1. The initial data (ρ, u)(x, t = 0) for ( 1) is said to be well-prepared when ρ(x, t = 0) = ρ (0) + ρ (1) + O 2 (x) with ρ (0) , ρ (1) real constants, u(x, t = 0) = u (0) (x, t = 0) + O ( ) (x) with ∇ • u (0) (x, t = 0) = 0. Therefore, as goes to 0, we expect that solutions (ρ, u) of (1) behave like ρ(x, t) = ρ (0) + ρ (1) + O 2 (x, t), u(x, t) = u (0) (x, t) + O ( ) (x, t) with u (0) solution of the Incompressible Euler system

∂ t u (0) + ∇ • (u (0) ⊗ u (0) ) + 1 ρ (0) ∇p (2) = 0, ∇ • u (0) = 0. (3)
Further details on the low Mach asymptotics can be found in A. Majda's book [START_REF] Majda | Compressible fluid flow and systems of conservation laws in several space variables[END_REF]. Rigorous justification dates back to [START_REF] Klainerman | Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids[END_REF][START_REF] Klainerman | Compressible and incompressible fluids[END_REF]. When the initial data are not well-prepared, initial layer might occur [START_REF] Ukai | The incompressible limit and the initial layer of the compressible Euler equation[END_REF]. We refer the interested reader to [START_REF] Alazard | A minicourse on the low Mach number limit[END_REF] for an up-to-date review detailing further results and indicating further relevant references on the analysis of the low Mach regimes.

Numerical issues of the low Mach regimes

There are several numerical issues that make the simulation of (1) challenging as becomes small:

• The stiffness of the pressure term induces severe stability conditions for standard timeexplicit schemes. Indeed, the characteristic speeds of the system behave like |u| ± 1 p (ρ) and thus become very large as → 0. It imposes to constrain the time step δt to be proportional to δx with δx the mesh size, which leads to non affordable computational costs.

• Moreover, the usual conservative finite-volume type schemes are also subject to a dramatic loss of accuracy when becomes small, so that it can be very difficult to capture the correct solution. This phenomena has been brought out in the seminal papers [START_REF] Guillard | On the behavior of upwind schemes in the low Mach number limit. II: Godunov type schemes[END_REF][START_REF] Guillard | On the behavior of upwind schemes in the low Mach number limit[END_REF] which explain how Godunov type schemes create spurious pressure waves that prevent the discrete solutions to be close to a discrete incompressible flow. This has been further analyzed in [START_REF] Barsukow | A numerical scheme for the compressible low-Mach number regime of ideal fluid dynamics[END_REF][START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF][START_REF] Dellacherie | Etude et discrétisation de modèles cinétiques et de modèles fluides à bas nombre de Mach[END_REF][START_REF] Dellacherie | Construction of modified Godunovtype schemes accurate at any Mach number for the compressible Euler system[END_REF][START_REF] Dellacherie | The influence of cell geometry on the Godunov scheme applied to the linear wave equation[END_REF]. It turns out that the loss of accuracy in fact highly depends on the geometry of the mesh [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit: a review[END_REF]. In particular UpWind schemes on Cartesian meshes (that we shall consider here) are inaccurate at low Mach numbers. These analysis have led to design suitable preconditioning techniques [START_REF] Guillard | On the behavior of upwind schemes in the low Mach number limit. II: Godunov type schemes[END_REF][START_REF] Guillard | On the behavior of upwind schemes in the low Mach number limit[END_REF][START_REF] Turkel | Preconditioned methods for solving the incompressible and low speed compressible equations[END_REF].

We are thus faced with the difficulty of designing numerical methods that preserve the asymptotic limit, without inducing too high computational costs. A possible avenue to cope with these difficulties is to try to develop "Asymptotic Preserving" methods, a term coined by S. Jin [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF][START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF].

Definition 1.2. Asymptotic Preserving (AP) Method.

A method is said to be AP when it preserves at the discrete level the asymptotic passage from a model to another. To be more specific, if the space and time steps δx and δt are kept fixed, the method automatically transforms to a stable discretization of the limit model when the small scale parameter tends to zero. Such methods are intended to be robust in the sense that the asymptotic regimes can be captured without resolving the small time and spatial scales. The efficiency of an AP scheme relies on careful time and space discretizations. Since we expect to make use of time steps much larger than the one imposed by the standard CFL conditions, it is clear that a part of the system should be treated implicitly to overcome the stiffness. However, the numerical cost of the implicit solver should remain as reduced as possible. We can mention [START_REF] Degond | All speed scheme for the low Mach number limit of the isentropic Euler equations[END_REF][START_REF] Dimarco | Study of a new asymptotic preserving scheme for the Euler system in the low Mach number limit[END_REF][START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF][START_REF] Klein | Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics. I. One-dimensional flow[END_REF][START_REF] Noelle | A weakly asymptotic preserving low Mach number scheme for the Euler equations of gas dynamics[END_REF] for attempts in that direction.

A new numerical strategy

In this work, we wish to develop an AP procedure by following the ideas introduced by J. Haack, S. Jin and J.-G. Liu in [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF]. The stiff Euler system is split into two parts: a nonlinear and hyperbolic system of conservation laws that involves only waves propagating with O(1) speeds, and a linear system that contains the fast acoustic dynamics. The former equations are treated explicitly, the latter are treated implicitly. Next, in this framework, we propose a space discretization that relies on strategies introduced in [START_REF] Berthelin | Kinetic schemes on staggered grids for barotropic Euler models: entropy-stability analysis[END_REF] and that presents the following originalities:

• we work on staggered grids, with densities (and thus pressures) and velocities stored on dual locations. This is not usual for hyperbolic conservation laws. However, this approach is well suited to avoid spurious oscillations in regimes close to incompressible: on Cartesian grids, we can expect to recover a final scheme based on the principles of MAC schemes [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF]. This approach is also developed, possibly with a different definition of numerical fluxes though, in e. g. [START_REF] Berthelin | Multifluid flows: a kinetic approach[END_REF][START_REF] Gastaldo | Staggered discretizations, pressure correction schemes and all speed barotropic flows[END_REF][START_REF] Herbin | Staggered schemes for all speed flows[END_REF][START_REF] Herbin | Explicit staggered schemes for the compressible Euler equations[END_REF][START_REF] Zaza | Contribution à la résolution numérique d'écoulements à tout nombre de Mach et au couplage fluide-poreux en vue de la simulation d'écoulements diphasiques homogénéisés dans les composants nucléaires[END_REF].

• the numerical fluxes are defined with an inspiration from kinetic schemes, see [START_REF] Coron | Numerical passage from kinetic to fluid equations[END_REF][START_REF] Deshpande | Kinetic theory based new upwind methods for inviscid compressible flows[END_REF] and the references therein. It provides simple formulae, which permit us to justify the positivity of the density and the entropy-stability [START_REF] Berthelin | Kinetic schemes on staggered grids for barotropic Euler models: entropy-stability analysis[END_REF]. We also refer the reader to [START_REF] Llobell | Staggered discretization for conservation laws of gas dynamics[END_REF] for further developments of the method, in particular its extension to second order (using muscl-like techniques and a RK2 scheme) and to the full Euler system.

Remark 1.3. We point out that the inspiration from the kinetic viewpoint has been recently used to develop new methods for solving the Incompressible Navier-Stokes system [START_REF] Bouchut | Second-order entropy satisfying BGK-FVS schemes for incompressible Navier-Stokes equations[END_REF], with interesting theoretical properties, in particular with respect to energy/entropy dissipation. Note however that, despite the derivation of this method relies on asymptotic analysis and the introduction of vanishing Mach numbers, it is not designed for performing well in low Mach regimes.

The paper is organized as follows. In Section 2.1 we describe the time splitting strategy of [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF]. In Section 2.2 and 2.3 we detail the construction of our scheme. Section 2.4 investigates the properties of the method: we show how the scheme becomes a (simple) incompressible solver as → 0 and we check that well-prepared data do not produce spurious oscillations, according to the analysis of [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF][START_REF] Dellacherie | Etude et discrétisation de modèles cinétiques et de modèles fluides à bas nombre de Mach[END_REF][START_REF] Dellacherie | The influence of cell geometry on the Godunov scheme applied to the linear wave equation[END_REF]. Finally, we discuss a set of numerical experiments in Section 3.

Numerical resolution

Splitting of the compressible Euler system

The compressible system (1) is split with two scales: the fast acoustics wave scale and the slow convection scale that contains the underlying incompressible dynamics. Following [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF], we rewrite the system in the following way

   ∂ t ρ + α∇ • (ρu) + (1 -α)∇ • (ρu) = 0, ∂ t (ρu) + ∇ • (ρu ⊗ u) + ∇ p(ρ)-a(t)ρ 2 + a(t) 2 ∇ρ = 0. (4)
The parameter 0 α < 1 is used in order to keep some momentum in the system for the fast dynamics, see (FD) below, which is necessary to enforce incompressibility. The splitting of the stiff pressure term also involves a new parameter a(t) > 0.

For the numerical resolution of (4) we adopt the following time-splitting strategy that makes the two distinct scales appear. We start with the Slow Dynamics system, the hyperbolic system of conservation laws that reads

   ∂ t ρ + α∇ • (ρu) = 0, ∂ t (ρu) + ∇ • (ρu ⊗ u) + ∇ p(ρ) -a(t)ρ 2 = 0. (SD)
Then, we consider the following Fast Dynamics system

   ∂ t ρ + (1 -α)∇ • (ρu) = 0, ∂ t (ρu) + a(t) 2 ∇ρ = 0. (FD)
Let us focus on the two-dimensional framework. Denoting u and v the components of the velocity field u = (u, v), the wave speeds of the (SD) in the horizontal (resp. vertical) direction are given by u (resp. v) and λ± (ρ, u) (resp. λ± (ρ, v)) where

λ± (ρ, ν) = ν ± c(ρ, ν) with c(ρ, ν) = (1 -α)ν 2 + α p (ρ) -a(t) 2 . ( 5
)
The wave speeds λ± (ρ, ν) are real for all ν ∈ R and ρ > 0 as soon as

a(t) min x∈Ω p ρ(x, t) .
The system is thus hyperbolic under this condition. It may seem natural to simply set

a(t) = min x∈Ω p ρ(x, t) . ( 6
)
But, with this choice, spurious oscillations are observed in some test cases for large values of . They appear in regions where the density is nearly uniform and the material velocity vanishes. Indeed, in these regions, the corresponding sound speed vanishes and the spurious oscillations are probably due to a lack of numerical diffusion in the Slow Dynamic part of the splitting. To overcome this difficulty, in these test cases (see Section 3.2), we set

a(t) = min x∈Ω p ρ(x, t) -ι(t) 2 > 0, (7)
where ι(t) is a constant smaller than the minimum of p (so that a(t) ι(t)(1 -2 ) > 0 for < 1). In the numerical tests of Section 3.2, we simply choose ι(t) = 0 or ι(t) = 1 when additional numerical diffusion is needed for large .

The choice of the splitting parameter α might impact the stability of the scheme: [START_REF] Zakerzadeh | A note on the stability of implicit-explicit flux-splittings for stiff systems of hyperbolic conservation laws[END_REF] prove that, for a linearized version of the system when 0 < 1, choosing α ∈ (0, 1 2 ) ensures the stability (in a sense defined in [START_REF] Majda | Stable viscosity matrices for systems of conservation laws[END_REF], see also [START_REF] Schütz | Flux splitting for stiff equations: A notion on stability[END_REF][START_REF] Zakerzadeh | A note on the stability of implicit-explicit flux-splittings for stiff systems of hyperbolic conservation laws[END_REF]) of the above splitting when it is combined with Rusanov-type numerical fluxes. In [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF], the numerical experiments are performed with α = 2 . As detailed below, we also find advantages in adopting this choice. It leads to

λ± (ρ, ν) = ν ± (1 -2 )ν 2 + p (ρ) -a(t).
and consequently since a(t) > 0, we have λ± (ρ, ν)

± 2ν ± + p (ρ), (8)
where the notation • ± is defined by x ± = |x| ± x /2, ∀x ∈ R. Whereas the wave speeds of the original system (1) behave like ν ± p (ρ)/ (as recalled in the introduction), we observe that the wave speeds of the system (SD) remain bounded when goes to 0 and thus ensures unconstrained CFL conditions with respect to (see Proposition 2.4).

Remark 2.1. As emphasized in [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF], the choice (6) ensures also that wave speeds of the system (SD) remain bounded when goes to 0 whatever the value of α. Indeed, if the initial conditions are well-prepared, we have

α 2 p (ρ) -min x∈Ω p ρ(x, t) = α ρ (2) (t, x) -min x∈Ω ρ (2) (t, x) p ρ (0) + O ( ) .
We now turn to the numerical resolution of those two systems. The characteristic speeds of the system (SD) are no longer stiff -in contrast to the wave speed of (1) -and it does not induce severe stability conditions. We shall make use of the space discretization set up in [START_REF] Berthelin | Kinetic schemes on staggered grids for barotropic Euler models: entropy-stability analysis[END_REF] for this system. The system (FD) is a simple linear hyperbolic system for the variable (ρ, q = ρu). However, it is not clear whether a direct resolution of the system can preserve the positivity of the density, a property which is absolutely crucial for the stability of the whole procedure, and for physical purposes. For this reason, we shall treat (FD) as a non linear system (for variables ρ and u) to be solved implicitly.

Time discretization and stability

From now on we use the following notation: given a quantity Q at the current time t we denote by Q its update at the next time t + δt. Let us detail the scheme by using the explicit Euler scheme for (SD) and the implicit Euler scheme for (FD), which can be summarized in the semi-discrete form:

       ρ * -ρ δt + α∇ • (ρu) = 0, ρ * u * -ρu δt + ∇ • (ρu ⊗ u) + ∇ p(ρ) -a(t)ρ 2 = 0,
and

     ρ -ρ * δt + (1 -α)∇ • (ρ u) = 0, ρ u -ρ * u * δt + a(t) 2 ∇ρ = 0.
Remark 2.2. In [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF], see also [START_REF] Noelle | A weakly asymptotic preserving low Mach number scheme for the Euler equations of gas dynamics[END_REF], the time discretization is constructed to reach second order accuracy: the two-level Adams-Bashforth scheme is used for the explicit slow system and the Crank-Nicolson scheme is used for the fast implicit system. Of course, it is straightforward to adapt such a time-discretization, but the analysis would be much more involved.

The construction proposed in [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF] relies on the following stability statement for split systems, see [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF]Lemma 3.1].

Lemma 2.3. If both methods

U n+1 -U n δt + AU n = 0 and U n+1 -U n δt + BU n+1 = 0,
are L 2 -stable, then the following method is also L 2 -stable

U n+1 -U n δt + AU n + BU n+1 = 0 in the sense that U n L 2 ≤ e Cnδt U 0 L 2
holds for a certain positive constant C.

Note however that this is a rough notion of stability. It gives an intuition on the relevance of the splitting. However, as said in introduction, it turns out that certain fast-slow splittings, that might look physically relevant at first sight with formal stability conditions independent of , give rise in practice to numerical instabilities. This is where the interplay between time and space discretizations arises: the origin of these instabilities is due to the checkerboard pressure profiles, a well-known difficulty in incompressible computational fluid mechanics. Since introducing stabilization terms in the scheme impacts the stability conditions and often leads to restriction depending on again, see in particular [34, sp. Section 3.5, Theorem 4.1 and Remark 4.2], we take advantages, in this article, of the use of staggered grids and its wellknown stability properties in the incompressible regime. The stability properties of the full discrete scheme are further analysed in Section 2.4.

Space discretization

We detail the scheme in 2D for a Cartesian grid. The computational domain is the square

Ω = [a x , b x ] × [a y , b y ] ⊂ R 2 .
The discretization strategy adopts the principles of the MAC schemes [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF], which are precisely suited to avoid the spurious modes in incompressible regimes. We proceed as follows:

• the primal mesh is given by subdividing each direction

a x = x 1 < x 2 < ... < x i-1 < x i < x i+1 < ... < x N < x N +1 = b x , a y = y 1 < y 2 < ... < y j-1 < y j < y j+1 < ... < y M < y M +1 = b y .
• Then, we define the midpoints

x i+1/2 = x i + x i+1 2
, ∀i ∈ {1, ..., N }, and y j+1/2 = y j + y j+1 2 , ∀j ∈ {1, ..., M }.

• We denote by δx i+1/2 , δy j+1/2 ,δx i and δy j the length of

[x i , x i+1 ], [y j , y j+1 ], [x i-1/2 , x i+1/2 ]
and [y j-1/2 , y j+1/2 ] respectively. We also denote by h = max(max i δx i+1/2 , max j δy j+1/2 ) the size of the mesh.

The discrete densities, the horizontal and the vertical velocities are stored at different locations, see Fig. 1:

• the density ρ is evaluated at the centers (x i+1/2 , y j+1/2 ) of the primal cells: the corresponding numerical unknowns are denoted by ρ i+1/2,j+1/2 ,

• the horizontal velocity u is evaluated on the vertical edges of the primal cells, at point (x i , y j+1/2 ): the numerical unknowns thus read u i,j+1/2

• the vertical velocity v is evaluated on the horizontal edges of the primal cells, at point (x i+1/2 , y j ): the numerical unknowns thus read v i+1/2,j .

ρ i+1/2,j+1/2 u i+1,j+1/2 × u i,j+1/2 × v i+1/2,j+1 × v i+1/2,j × i i + 1/2 i + 1 j j + 1/2 j + 1
Figure 1: Position of the unknowns on a MAC grid

The density can be naturally defined at the edges of the primal mesh by averaging as follows

ρ i,j+ 1 2 = δx i+ 1 2 ρ i+ 1 2 ,j+ 1 2 + δx i-1 2 ρ i-1 2 ,j+ 1 2 2δx i , ρ i+ 1 2 ,j = δy j+ 1 2 ρ i+ 1 2 ,j+ 1 2 + δy j-1 2 ρ i+ 1 2 ,j-1 2 2δy j .
Indeed, the different sets of unknows can be thought of as piecewise constant functions on different meshes: the primal unknowns

ρ i+ 1 2 ,j+ 1 2 can be considered as constant on primal cells [x i , x i+1 ] × [y j , y j+1 ] whereas the dual unknowns ρ i,j+ 1 2 , resp. ρ i+ 1 2 ,j , can considered as constant on dual cells [x i-1 2 , x i+ 1 2 ] × [y j , y j+1 ], resp. [x i , x i+1 ] × [y j-1 2 , y j+ 1 2 ]. The above definitions of ρ i,j+ 1 2 , resp. ρ i+ 1 2
,j , ensure that the integral of the discrete densities over the dual cell

[x i-1 2 , x i+ 1 2 ] × [y j , y j+1 ], resp. [x i , x i+1 ] × [y j-1 2 , y j+ 1 2
], is the same whatever the considered set of unknowns (that is the primal unknowns ρ i+ 1 2 ,j+ 1 2 or the dual unknowns ρ i,j+ 1 2 , resp. ρ i+ 1 2 ,j ). These expressions are different from linear interpolation but bearing in mind that

x i+ 1 2 is the midpoint of [x i , x i+1 ], we have δx i-1 2 + δx i+ 1 2 = 2δx i (and similarly δy j-1 2 + δy j+ 1 2 = 2δy j ) so that ρ i,j+ 1 2
and ρ i+ 1 2 ,j are convex combinations of the densities stored at the centers of the primal cells.

We solve (SD) in order to produce an intermediate solution (

ρ * i+1/2,j+1/2 , u * i,j+1/2 , v * i+1/2,j )
with the following scheme

                                                         ρ * i+ 1 2 ,j+ 1 2 -ρ i+ 1 2 ,j+ 1 2 δt + α   F x i+1,j+ 1 2 -F x i,j+ 1 2 δx i+ 1 2 + F y i+ 1 2 ,j+1 -F y i+ 1 2 ,j δy j+ 1 2   = 0, ρ * i,j+ 1 2 u * i,j+ 1 2 -ρ i,j+ 1 2 u i,j+ 1 2 δt + G u,x i+ 1 2 ,j+ 1 2 -G u,x i-1 2 ,j+ 1 2 δx i + G u,y i,j+1 -G u,y i,j δy j+ 1 2 + Π i+ 1 2 ,j+ 1 2 -Π i-1 2 ,j+ 1 2 δx i = 0, ρ * i+ 1 2 ,j v * i+ 1 2 ,j -ρ i+ 1 2 ,j v i+ 1 2 ,j δt + G v,x i+1,j -G v,x i,j δx i+ 1 2 + G v,y i+ 1 2 ,j+ 1 2 -G v,y i+ 1 2 ,j-1 2 δy j + Π i+ 1 2 ,j+ 1 2 -Π i+ 1 2 ,j-1 2 δy j = 0, (SDd)
where we are going to define the numerical fluxes in these formulae. The system (SD) differs from the usual Euler system. Nevertheless, we can readily adapt the scheme proposed in [START_REF] Berthelin | Kinetic schemes on staggered grids for barotropic Euler models: entropy-stability analysis[END_REF][START_REF] Llobell | Staggered discretization for conservation laws of gas dynamics[END_REF].

Bearing in mind that the sound speed c(ρ, ν) defined in ( 5) now depends on ρ and ν a velocity component u or v, we get the following definition for the fluxes. We introduce

F + (ρ, ν) =        0 if ν -c(ρ, ν), ρ 4 c(ρ, ν) ν + c(ρ, ν) 2 if |ν| c(ρ, ν), ρν if ν c(ρ, ν),
and

F -(ρ, ν) =        ρν if ν -c(ρ, ν), - ρ 4 c(ρ, ν) (ν -c(ρ, ν)) 2 if |ν| c(ρ, ν), 0 if ν c(ρ, ν).
Note that, with this convention, we have ±F ± (ρ, ν) 0. Then we set

F x i,j+ 1 2 = F x,+ i,j+ 1 2 + F x,- i,j+ 1 2 and F y i+ 1 2 ,j = F y,+ i+ 1 2 ,j + F y,- i+ 1 2 ,j , with F x,+ i,j+ 1 2 = F + (ρ i-1 2 ,j+ 1 2 , u i,j+ 1 2 ) and F x,- i,j+ 1 2 = F -(ρ i+ 1 2 ,j+ 1 2 , u i,j+ 1 2
),

F y,+ i+ 1 2 ,j = F + (ρ i+ 1 2 ,j-1 2 , v i+ 1 2 ,j ) and F y,- i+ 1 2 ,j = F -(ρ i+ 1 2 ,j+ 1 2 , v i+ 1 2 ,j ).
The definition of the convection fluxes in the momentum equation is given by

G u,x i+ 1 2 ,j+ 1 2 = u i,j+ 1 2 F x,+ i+ 1 2 ,j+ 1 2 + u i+1,j+ 1 2 F x,- i+ 1 2 ,j+ 1 2 and G u,y i,j = u i,j-1 2 F y,+ i,j + u i,j+ 1 2 F y,- i,j
and a similar definition for G v,x i,j and G u,y i+ 1 2 ,j+ 1

2

. It uses the averaged mass fluxes

F x,± i+ 1 2 ,j+ 1 2 = 1 2 F x,± i,j+ 1 2 + F x,± i+1,j+ 1 2 and F y,± i,j = δx i+ 1 2 F y,± i+ 1 2 ,j + δx i-1 2 F y,± i-1 2 ,j 2δx i .
The discrete analog of the parameter a in ( 7) is the discrete quantity

a d = min i,j p (ρ i+ 1 2 ,j+ 1 2 ) -ι 2 , ( 9
)
and thus the pressure term Π i+ 1 2 ,j+ 1 2 -which differs from the usual Euler case dealt with in [3] -is defined by

Π i+ 1 2 ,j+ 1 2 = p(ρ i+ 1 2 ,j+ 1 2 ) -a d ρ i+ 1 2 ,j+ 1 2 2
.

We now turn to the resolution of the system (FD) to find (ρ

i+ 1 2 ,j+ 1 2 , u i,j+ 1 2
, v i+ 1 2 ,j ), using the solution of (SDd) we have just obtained as initial data. Working on staggered grids offers us the possibility to use a centered scheme for the gradient of the pressure term, which is the key to overcome the difficulties recorded with methods on colocalized meshes, see [START_REF] Dellacherie | Etude et discrétisation de modèles cinétiques et de modèles fluides à bas nombre de Mach[END_REF][START_REF] Noelle | A weakly asymptotic preserving low Mach number scheme for the Euler equations of gas dynamics[END_REF]. We set

                         ρ i+ 1 2 ,j+ 1 2 -ρ * i+ 1 2 ,j+ 1 2 δt + (1 -α)    F U p,x i+1,j+ 1 2 -F U p,x i,j+ 1 2 δx i+ 1 2 + F U p,y i+ 1 2 ,j+1 -F U p,y i+ 1 2 ,j δy j+ 1 2    = 0, ρ i,j+ 1 2 u i,j+ 1 2 -ρ * i,j+ 1 2 u * i,j+ 1 2 δt + a d 2 ρ i+ 1 2 ,j+ 1 2 -ρ i-1 2 ,j+ 1 2 δx i = 0, ρ i+ 1 2 ,j v i+ 1 2 ,j -ρ * i+ 1 2 ,j v * i+ 1 2 ,j δt + a d 2 ρ i+ 1 2 ,j+ 1 2 -ρ i+ 1 2 ,j-1 2 δy j = 0. (FDd)
The upwind fluxes F U p,x and F U p,y are obtained by applying the upwinding principle based on the sign of the velocity. Precisely, we use in (FDd) the following implicit fluxes

F U p,x i,j+ 1 2 = ρ i-1 2 ,j+ 1 2 [u i,j+ 1 2 ] + -ρ i+ 1 2 ,j+ 1 2 [u i,j+ 1 2 ] - and F U p,y i+ 1 2 ,j = ρ i+ 1 2 ,j-1 2 [v i+ 1 2 ,j ] + -ρ i+ 1 2 ,j+ 1 2 [v i+ 1 2 ,j ] -.
Here, we are still using the definition x ± = |x| ± x /2, ∀x ∈ R. In contrast to [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF], we are working here on staggered grids and the resolution does not involve the momentum as a numerical unknown but instead the velocity. Despite the fact that (FD) is a linear system for the conservative quantities (ρ, q = ρu), the corresponding discrete equation (FDd) becomes a nonlinear problem for the non-conservative variables. Accordingly, the update does not come from the resolution of a linear system, and, instead, we have to use a root-finding algorithm as the Newton-Raphson method. This is absolutely crucial in order to preserve the positivity of the density (stated below, see Proposition 2.4). Nevertheless, it is important to note that, even if the system looks fully implicit, the expressions for u i,j+1/2 and v i+1/2,j as a function of ρ i+1/2,j+1/2 can be readily obtained from the last two equations of (FDd) so that the system can be solved as a non linear system of scalar equations involving only the unknowns (ρ i+1/2,j+1/2 ) with i ∈ {1, ..., N }, j ∈ {1, ..., M }. This is of great importance to save computational cost.

In terms of stability let us give the following statement, which justifies that the scheme produces physically relevant quantities with a stability condition that does not degenerate when becomes small. Proposition 2.4. We set α = 2 . Suppose that the data (ρ, u) is well-prepared, in the sense of Definition 2.6 and such that ρ j+ 1 2 0 for all j. We assume the following CFL-like condition

δt δx i+ 1 2 u i,j+ 1 2 - + u i+1,j+ 1 2 + + p (ρ i+ 1 2 ,j+ 1 2 ) + δt δy j+ 1 2 v i+ 1 2 ,j - + v i+ 1 2 ,j+1 + + p (ρ i+ 1 2 ,j+ 1 2 ) 1 2 , ( 10
)
then the scheme preserves the non-negativity of the density: ρ j+ 1 2 0 for all j.

Proof. We adapt readily the arguments in [START_REF] Llobell | Staggered discretization for conservation laws of gas dynamics[END_REF] to justify that (SDd) preserves positive densities, under the following CFL-like condition

δt δx i+ 1 2 λ-(ρ i+ 1 2 ,j+ 1 2 , u i,j+ 1 2 ) - + λ+ (ρ i+ 1 2 ,j+ 1 2 , u i+1,j+ 1 2 ) + + δt δy j+ 1 2 λ-(ρ i+ 1 2 ,j+ 1 2 , v i+ 1 2 ,j ) - + λ+ (ρ i+ 1 2 ,j+ 1 2 , v i+ 1 2 ,j+1 ) + 1.
Owing to (8) which holds since α = 2 , this condition is ensured when ( 10) is satisfied. Thus, it proves that ρ *

i+ 1 2 ,j+ 1 2
> 0 for all i, j. We now turn to the resolution of the Fast Dynamics (FDd). The monotonicity of the implicit Upwind scheme [START_REF] Boyer | Analysis of the upwind finite volume method for general initial-and boundaryvalue transport problems[END_REF] ensures that

ρ i+ 1 2 ,j+ 1 2 > 0 as soon as ρ * i+ 1 2 ,j+ 1 2 > 0. We conclude that ρ i+ 1 2 ,j+ 1 2
> 0 at every time step under the CFL-like condition [START_REF] Dellacherie | Etude et discrétisation de modèles cinétiques et de modèles fluides à bas nombre de Mach[END_REF].

Low Mach regimes and AP-scheme

We wish now to show that the scheme performs well in low Mach regimes and exhibits the AP features. Let us suppose that the initial data of the system are well prepared in the sense of Definition 1.1. As said above, the solutions of the compressible Euler system (1) converge to a constant density and a velocity field solution of the Incompressible Euler system (3). Our purpose is two-fold:

1) to show that our solver automatically becomes an incompressible solver for (3) in the low Mach number regime 0 < 1. It makes the AP character of the scheme precise.

2) to justify that the scheme preserves well prepared data, in the discrete sense given in Definition 2.6 below. This crucial property justifies that the proposed discretization does not produce spurious oscillations that can prevent to capture the correct solution in low Mach regimes, see [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF][START_REF] Dellacherie | Etude et discrétisation de modèles cinétiques et de modèles fluides à bas nombre de Mach[END_REF][START_REF] Dellacherie | The influence of cell geometry on the Godunov scheme applied to the linear wave equation[END_REF][START_REF] Guillard | On the behavior of upwind schemes in the low Mach number limit[END_REF][START_REF] Guillard | On the behavior of upwind schemes in the low Mach number limit. II: Godunov type schemes[END_REF][START_REF] Noelle | A weakly asymptotic preserving low Mach number scheme for the Euler equations of gas dynamics[END_REF].

In order to bring out that the behaviour of the continuous case is well reproduced, we introduce an ansatz of the discrete solution:

• ρ i+ 1 2 ,j+ 1 2 = ρ (0) i+ 1 2 ,j+ 1 2 + ρ (1) i+ 1 2 ,j+ 1 2 + 2 ρ (2) i+ 1 2 ,j+ 1 2 + ... • u i,j+ 1 2 = u (0) i,j+ 1 2 + u (1) i,j+ 1 2 + 2 u (2) i,j+ 1 2 + ... and v i+ 1 2 ,j = v (0) i+ 1 2 ,j + v (1) i+ 1 2 ,j + 2 v (2) i+ 1 2 ,j + ... • p i+ 1 2 ,j+ 1 2 = p (0) i+ 1 2 ,j+ 1 2 + p (1) i+ 1 2 ,j+ 1 2 + 2 p (2) i+ 1 2 ,j+ 1 2 + ...
The investigation of the asymptotic regime → 0 will make the following discrete operators appear.

Definition 2.5. We define the discrete divergence operator ∇ d • of a vector u = (u, v) as

∇ d • u i+ 1 2 ,j+ 1 2 = u i+1,j+ 1 2 -u i,j+ 1 2 δx i+ 1 2 + v i+ 1 2 ,j+1 -v i+ 1 2 ,j δy j+ 1 2 ,
and the discrete Laplacian operator ∆ d of a scalar quantity p as

∆ d p i+ 1 2 ,j+ 1 2 = 1 δx i+ 1 2 p i+ 3 2 ,j+ 1 2 -p i+ 1 2 ,j+ 1 2 δx i+1 - p i+ 1 2 ,j+ 1 2 -p i-1 2 ,j+ 1 2 δx i + 1 δy j+ 1 2 p i+ 1 2 ,j+ 3 2 -p i+ 1 2 ,j+ 1 2 δy j+1 - p i+ 1 2 ,j+ 1 2 -p i+ 1 2 ,j-1 2 δy j .
Since the preparation of data is an important assumption for the analysis of the asymptotic behaviour, it is convenient to introduce an equivalent definition at the discrete level. Definition 2.6. Discrete data (ρ, u = (u, v)) are said to be well-prepared if they satisfy, for all i, j, (1) + O( 2 ) with ρ (0) , ρ (1) real constants,

• ρ i+ 1 2 ,j+ 1 2 = ρ (0) + ρ
• ∇ d • u (0) i+1/2,j+1/2 = 0, and thus ∇ d • u i+1/2,j+1/2 = O ( ) .
This definition appeared in [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF][START_REF] Dellacherie | Etude et discrétisation de modèles cinétiques et de modèles fluides à bas nombre de Mach[END_REF][START_REF] Dellacherie | The influence of cell geometry on the Godunov scheme applied to the linear wave equation[END_REF]; it plays a crucial role in the analysis of the sensitivity of numerical solvers to the grid geometry. It is the discrete counterpart of Definition 1.1.

Let us suppose that the discrete solution is well-prepared in the sense of Definition 2.6 from time t = 0 up to time t > 0. We wish to show that the updated physical quantities still satisfy the conditions of Definition 2.6.

The first thing to notice is that the explicit pressure term drops out, as the modified explicit pressure has a simple expansion for small 's.

Lemma 2.7. Let (ρ, u) be well-prepared. The modified pressure p(ρ) -a d ρ in (SDd) admits the expansion

Π i+ 1 2 ,j+ 1 2 := p ρ i+ 1 2 ,j+ 1 2 -a d ρ i+ 1 2 ,j+ 1 2 2 = C -2 2 + C -1 + C 0 + O ( )
where the coefficients C -2 , C -1 , C 0 do not depend on the considered cell.

Proof. Since the data are well-prepared, the pressure

p i+ 1 2 ,j+ 1 2 = p(ρ i+ 1 2 ,j+ 1 2
) expands as

p i+ 1 2 ,j+ 1 2 = p ρ (0) + ρ (1) p ρ (0) + 2 ρ (2) i+ 1 2 ,j+ 1 2 p ρ (0) + ρ (1) 2 2 p ρ (0) + O 3 . ( 11
)
Hence we obtain

p i+ 1 2 ,j+ 1 2 -a d ρ i+ 1 2 ,j+ 1 2 2 = p ρ (0) -a (0) d ρ (0) 2 + ρ (1) p ρ (0) -a (0) d ρ (1) -a (1) d ρ (0) + ρ (2) i+ 1 2 ,j+ 1 2 p ρ (0) + ρ (1) 2 2 p ρ (0) -a (0) d ρ (2) i+ 1 2 ,j+ 1 2 -a (1)
d ρ (1) -a

(2)

d ρ (0) + O ( ) .

An expansion of the coefficient a

d = min i,j p (ρ i+ 1 2 ,j+ 1 2 ) -ι 2 : a (0) d = p ρ (0) , a (1)
d = ρ (1) p ρ (0) , a

(2)

d = min i,j ρ (2) i+ 1 2 ,j+ 1 2 p ρ (0) + ρ (1) 2 2 p ρ (0) -ι, (12)
leads to the conclusion with the following expressions of the constants

C -2 , C -1 , C 0 C 0 = - ρ (1) 2 2 p ρ (0) + ρ (0) p ρ (0) -min i,j ρ (2) i+ 1 2 ,j+ 1 2 ρ (0) p ρ (0) + ιρ (0) , C -1 = -ρ (0) ρ (1) p ρ (0) , C -2 = p ρ (0) -ρ (0) p ρ (0) .
As we did for the continuous system, we now insert the expansion in the scheme obtained by taking the sum of (SDd) and (FDd)

ρ i+ 1 2 ,j+ 1 2 -ρ i+ 1 2 ,j+ 1 2 δt + α   F x i+1,j+ 1 2 -F x i,j+ 1 2 δx i+ 1 2 + F y i+ 1 2 ,j+1 -F y i+ 1 2 ,j δy j+ 1 2   +(1 -α)    F U p,x i+1,j+ 1 2 -F U p,x i,j+ 1 2 δx i+ 1 2 + F U p,y i+ 1 2 ,j+1 -F U p,y i+ 1 2 ,j δy j+ 1 2    = 0, ( 13
)
ρ i,j+ 1 2 u i,j+ 1 2 -ρ i,j+ 1 2 u i,j+ 1 2 δt + G u,x i+ 1 2 ,j+ 1 2 -G u,x i-1 2 ,j+ 1 2 δx i + G u,y i,j+1 -G u,y i,j δy j+ 1 2 + Π i+ 1 2 ,j+ 1 2 -Π i-1 2 ,j+ 1 2 δx i + a d 2 ρ i+ 1 2 ,j+ 1 2 -ρ i- ,j+ 1 2 δx i = 0, (14)
ρ i+ 1 2 ,j v i+ 1 2 ,j -ρ i+ 1 2 ,j v i+ 1 2 ,j δt + G v,x i+1,j -G v,x i,j δx i+ 1 2 + G v,y i+ 1 2 ,j+ 1 2 -G v,y i+ 1 2 ,j-1 2 δy j + Π i+ 1 2 ,j+ 1 2 -Π i+ 1 2 ,j-1 2 δy j + a d 2 ρ i+ 1 2 ,j+ 1 2 -ρ i+ ,j-1 2 δy j = 0, (15)
and we identify terms arising with the same power of .

1) The leading order term in the momentum equations ( 14) and ( 15) scales as O 1 2 . It comes from the pressure gradient. However, owing to Lemma 2.7, the explicit pressure gradient defined from Π i+ 1 2 ,j+ 1 2 contains only O( ) terms. It leads to

ρ (0) i+ 1 2 ,j+ 1 2 -ρ (0) i-1 2 ,j+ 1 2 δx i = 0 and ρ (0) i+ 1 2 ,j+ 1 2 -ρ (0) i+ 1 2 ,j-1 2 δy j = 0.
From this, we deduce that, for all (i, j), ρ

(0) i-1 2 ,j+ 1 2 = ρ (0) i+ 1 2 ,j+ 1 2 = ρ (0) i+ 1 2 ,j-1 2
. In other words, at the discrete level the leading term in the updated density ρ (0) is constant in space.

2) Then, at order O 1 , we similarly get

ρ (1) i+ 1 2 ,j+ 1 2 -ρ (1) i-1 2 ,j+ 1 2 δx i = 0 and ρ (1) i+ 1 2 ,j+ 1 2 -ρ (1) i+ 1 2 ,j-1 2 δy j = 0.
We deduce that, for all (i, j), ρ

(1) i-1 2 ,j+ 1 2 = ρ (1) i+ 1 2 ,j+ 1 2 = ρ (1) i+ 1 2 ,j-1 2
: at the discrete level ρ (1) is constant in space, too.

3) Let us now consider terms of order O (1). Since the data (ρ, u) is assumed to be wellprepared, by multiplying the mass equation ( 13) by δx i+ 1 2 δy j+ 1 2 and summing over i and j (with periodic or wall boundary conditions and denoting by N, M the number of grid points in the horizontal and vertical direction, respectively), we get

M j=1 N i=1 δx i+ 1 2 δy j+ 1 2 ρ (0) -ρ (0) δt = 0.
We remind the reader that

M j=1 N i=1 δx i+ 1 2 δy j+ 1 2 = |Ω|. Therefore ρ (0) = ρ (0)
: the discrete density at leading order is constant in space and time.

At the order O (1), since the density ρ (0) and the updated density ρ (0) at leading order are constant in space and equal, the mass fluxes become

F x,(0) i,j+ 1 2 = F + (ρ (0) , u (0) i,j+ 1 2 ) + F -(ρ (0) , u (0) i,j+ 1 2 ) = ρ (0) u (0) i,j+ 1 2 , similarly F y,(0) i+ 1 2 ,j = ρ (0) v (0) i+ 1 2 ,j , F U p,x i,j+ 1 2 = ρ (0) u (0) i,j+ 1 2 and F U p,y i+ 1 2 ,j = ρ (0) v (0) i+ 1
2 ,j . Going back to the discrete mass balance equation ( 13) it becomes at the order O (1)

αρ (0) ∇ d • u (0) i+ 1 2 ,j+ 1 2 + (1 -α)ρ (0) ∇ d • u (0) i+ 1 2 ,j+ 1 2 = 0. Since ∇ d • u (0) i+ 1 2 ,j+ 1 2
vanishes for all (i, j), the update also satisfies ∇ d • u (0) = 0 and the second property in Definition 2.6 is preserved in the low Mach regime.

Remark 2.8. When the initial data u(x, 0) = u init (x) satisfies ∇ • u init,(0) (x) = 0 at the continuous level, the simple evaluation on the grid

∇ d • u init,(0) i+ 1 2 ,j+ 1 2 = u init,(0) (x i+1 , y j+ 1 2 ) -u init,(0) (x i , y j+ 1 2 ) δx i+ 1 2 + v init,(0) (x i+ 1 2 , y j+1 ) -v init,(0) (x i+ 1 2 , y j ) δy j+ 1 2
is of order O(h). This consistency error can propagate, and the obtained discrete velocity field is not discrete-divergence-free. In particular, if the mesh does not resolve the small scale (δx ), the discrete initial data does not fulfill the criterion in Definition 2.6. Obviously, it is possible to adopt more involved definitions of u init i+ 1 2 ,j+ 1 2 , for instance based on a discrete Helmholtz decomposition to get rid of this issue but we can also take advantages to set α = 2 .

This choice, which already appears in [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF], modifies the expansion above and the identification of the O(1) and O( 2 ) terms in the discrete mass balance equation. In particular, with this choice, the method automatically projects the velocity onto the space of discretedivergence-free fields when → 0: indeed, in this case, at the order O (1), the discrete mass balance equation becomes

∇ d • u (0) i+ 1 2 ,j+ 1 2 = 0.
4) At the order O ( ), we multiply the mass equation ( 13) by δx i+ 1 2 δy j+ 1 2 and sum for i, j and, since ρ (1) is a constant in space, we obtain that ρ (1) = ρ (1) . The first property in Definition 2.6 is satisfied in the low Mach regime.

Remark 2.9. Using that ρ (k) = ρ (k) are constant for k ∈ {0, 1}, we can expand the mass flux F x i,j+ 1 2 as

F x i,j+ 1 2 = F + ρ (0) + ρ (1) , u (0) i,j+ 1 2 + u (1) i,j+ 1 2 + F -ρ (0) + ρ (1) , u (0) i,j+ 1 2 + u (1) i,j+ 1 2 + O( 2 ) = ρ (0) + ρ (1) u (0) i,j+ 1 2 + ρ (0) u (1) i,j+ 1 2 + O( 2 )
, so that we obtain the following expression

F x,(1) i,j+ 1 2 = ρ (1) u (0) i,j+ 1 2 + ρ (0) u (1) i,j+ 1 2 .
Similarly, we have

F y,(1) i+ 1 2 ,j = ρ (1) v (0) i+ 1 2 ,j + ρ (0) v (1) i+ 1 2 ,j , F U p,x,(1) i,j+ 1 2 = ρ (1) u (0) i,j+ 1 2 + ρ (0) u (1) i,j+ 1 2 , F U p,y,(1) i+ 1 2 ,j = ρ (1) v (0) i+ 1 2 ,j + ρ (0) v (1) i+ 1
2 ,j . Thus, when α = 2 , the first order term O ( ) in the mass equation ( 13) reads:

ρ (1) ∇ d • u (0) i+1/2,j+1/2 + ρ (0) ∇ d • u (1)
i+1/2,j+1/2 = 0.

Since, as said above in Remark 2.8, in this case we have ∇ d • u (0) = 0, we conclude that ∇ d • u (1) = 0, in analogy to the continuous case. 5) At the order O (1), using that ρ (0) and ρ (1) are constant, the momentum equation [START_REF] Dimarco | Study of a new asymptotic preserving scheme for the Euler system in the low Mach number limit[END_REF] gives us:

ρ (0) u (0) i,j+ 1 2 -u (0) i,j+ 1 2 δt + G u,x,(0) i+ 1 2 ,j+ 1 2 -G u,x,(0) i-1 2 ,j+ 1 2 δx i + G u,y,(0) i,j+1 -G u,y,(0) i,j δy j+ 1 2 + p (2) i+ 1 2 ,j+ 1 2 -p (2) i-1 2 ,j+ 1 2 δx i = 0,
where p

(2)

i+ 1 2 ,j+ 1 2 is the O( 2 ) term in the expansion of p(ρ i+ 1 2 ,j+ 1 2
) (which is similar to [START_REF] Dellacherie | Construction of modified Godunovtype schemes accurate at any Mach number for the compressible Euler system[END_REF]). Indeed, bearing in mind the expansion (12) of a d , we have p

(2) i+ 1 2 ,j+ 1 2 -p (2) i-1 2 ,j+ 1 2 = a (0) d ρ (2) i+ 1 2 ,j+ 1 2 -ρ (2) i-1 2 ,j+ 1 2 .
We now turn to give a simple expression of the fluxes G Going back to the definition of F ± we thus observe that

F +,(0) (ρ, ν) = ρ (0) [ν (0) ] + and F -,(0) (ρ, ν) = -ρ (0) [ν (0) ] -.
Hence we get G u,x,(0) = ρ (0) G u,x and G u,y,(0) = ρ (0) G u,y with

G u,x i+ 1 2 ,j+ 1 2 = {u (0) } x i+ 1 2 ,j+ 1 2 • {u (0) } x i+ 1 2 ,j+ 1 2 -{|u (0) |} x i+ 1 2 ,j+ 1 2 • u (0) x i+ 1 2 ,j+ 1 2 (16)
and

G u,y i,j = {δxv (0) } x i,j δx i • {u (0) } y i,j - {|δxv (0) |} x i,j δx i • u (0) y i,j , (17)
where we have introduced the following notations for the average and the jump of a quantity Q, for k and l in N or N + 1 2 N,

{Q} x k,l = Q k-1 2 ,l + Q k+ 1 2 ,l 2 and {Q} y k,l = Q k,l-1 2 + Q k,l+ 1 2 2 , Q x k,l = Q k+ 1 2 ,l -Q k-1 2 ,l 2 and Q y k,l = Q k,l+ 1 2 -Q k,l-1 2 2 .
The fluxes G u,x i+ 1 2 ,j+ 1 2 (resp. G u,y i,j ) can be interpreted as fluxes in the horizontal (resp. vertical) direction of the advected quantity u (0) at velocity u (0) (resp. v (0) ). We obtain a similar expression of the fluxes for the second component of the velocity v:

G v,y i+ 1 2 ,j+ 1 2 = {v (0) } y i+ 1 2 ,j+ 1 2 • {v (0) } y i+ 1 2 ,j+ 1 2 -{|v (0) |} y i+ 1 2 ,j+ 1 2 • v (0) y i+ 1 2 ,j+ 1 2 (18)
and

G v,x i,j = {δyu (0) } y i,j δy j • {v (0) } x i,j - {|δyu (0) |} y i,j δy j • v (0) x i,j . ( 19
)
Thus, with similar calculation for the vertical component, we arrive at the following result Proposition 2.10. When setting α = 2 , formally, as tends to zero, the compressible Euler solver (SDd)-(FDd) behaves like the following scheme for solving the Incompressible Euler system (3)

                                                   u (0) i,j+ 1 2 -u (0) i,j+ 1 2 δt + G u,x i+ 1 2 ,j+ 1 2 -G u,x i-1 2 ,j+ 1 2 δx i + G u,y i,j+1 -G ,u,y i,j δy j+ 1 2 + 1 ρ (0) p (2) i+ 1 2 ,j+ 1 2 -p (2) i-1 2 ,j+ 1 2 δx i = 0, v (0) i+ 1 2 ,j -v (0) i+ 1 2 ,j δt + G v,y i+ 1 2 ,j+ 1 2 -G v,y i+ 1 2 ,j-1 2 δy j + G v,x i+1,j -G ,v,x i,j δx i+ 1 2 + 1 ρ (0) p (2) i+ 1 2 ,j+ 1 2 -p (2) i+ 1 2 ,j-1 2 δy j = 0, ∇ d • u (0) i+ 1 2 ,j+ 1 2 = 0, ( 20
)
with the definitions ( 16), ( 17), ( 18), ( 19) for fluxes G.

Remark 2.11. The 1D version of the fluxes G for the linear transport equation with constant velocity is exactly the upwind fluxes. In the case of the 1D Burgers equation

∂ t u + ∂ x u 2 2
= 0, the analogue of the fluxes G would be the following fluxes:

G i+ 1 2 = 1 4          (u i + u i+1 )u i if u i > 0 and u i+1 > 0, (u i + u i+1 )u i+1 if u i < 0 and u i+1 < 0, 2u i u i+1 if u i < 0 < u i+1 , u 2 i+1 + u 2 i if u i > 0 > u i+1 .
Numerical experiments show that these fluxes do not produce non entropic shocks, in contrast to the upwind scheme which is known to be non entropic since it captures wrong stationary shock solutions [START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF]Sections 13 & 14].

In the limit scheme [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF], the divergence free constraint on the velocity can also be written as an elliptic problem for the pressure which reads as follows:

∆ d p (2) i+ 1 2 ,j+ 1 2 = -ρ (0)   U i+1,j+ 1 2 -U i,j+ 1 2 δx i+ 1 2 + V i+ 1 2 ,j+1 -V i+ 1 2 ,j δy j+ 1 2   , ( 21
)
where U i,j+ 1 2 and V i+ 1 2 ,j are defined by

U i,j+ 1 2 = G u,x i+ 1 2 ,j+ 1 2 -G u,x i-1 2 ,j+ 1 2 δx i + G u,y i,j+1 -G u,y i,j δy j+ 1 2 , V i+ 1 2 ,j = G v,y i+ 1 2 ,j+ 1 2 -G v,y i+ 1 2 ,j-1 2 δy j + G v,x i+1,j -G v,x i,j δx i+ 1 2 .
Note that ∆ d is the standard 5-points discretization of the Laplace operator, in contrast to the operator that arises in the original scheme [20, eq. ( 28)] which decouples odd/even points of the discretization. This is a well-known drawback of colocalized methods which can lead to spurious oscillations due to the nontrivial checkerboard modes that belong to the kernel of the operator. Using staggered grids avoids this difficulty; it can be related to the MAC discretization, which gets rid of these spurious modes [START_REF] Dellacherie | The influence of cell geometry on the Godunov scheme applied to the linear wave equation[END_REF].

Note that, as mentioned in Remark 2.8, the limit scheme (written here only in the case α = 2 ) automatically projects the initial velocity to discrete-divergence-free velocity fields. Actually, for the first time step the elliptic problem ( 21) should be written

∆ d p (2) i+ 1 2 ,j+ 1 2 = -ρ (0)   U i+1,j+ 1 2 -U i,j+ 1 2 δx i+ 1 2 + V i+ 1 2 ,j+1 -V i+ 1 2 ,j δy j+ 1 2   + ∇ d • u (0) i+ 1 2 ,j+ 1 2 δt ,
which enforces ∇ d • u (0) = 0 as required in the limit scheme [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF]

even if ∇ d • u (0) = 0.

Numerical simulations

In this section we present some numerical test cases, in 1D and 2D, to illustrate the performances of the proposed scheme. In particular, we point out in the numerical experiments the stability and consistency issues.

As suggested above, we use α = 2 for all numerical experiments and otherwise specified, we adopt in the numerical tests the definition [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF] with ι = 0.

Simulation of 1D Riemann problems in low Mach regime

We start with a simple 1D Riemann problem, issued from [START_REF] Dimarco | Study of a new asymptotic preserving scheme for the Euler system in the low Mach number limit[END_REF]. The solution is made of two moving rarefaction waves, a left and a right one, separated by a constant state. The pressure law is given by p(ρ) = ρ 2 and the initial conditions are

ρ (0) (x) = 1 + 2 if x < 0.5, 1 if x > 0.5, and
u (0) (x) = 1 - if x < 0.5, 1 + if x > 0.5.
We compute the solution on the interval (0, 1) for three different values of ( , T ) where T is the final time. For each case, the number of grid points is N = 200 and the time step is defined by δt = βδx with β = 0.2, 0.1 or 0.01. We compare the results produced by the AP scheme described in this paper with the results given by the first order explicit scheme of [START_REF] Berthelin | Kinetic schemes on staggered grids for barotropic Euler models: entropy-stability analysis[END_REF].

We first choose close to one ( = √ 0.99) and thus the Fast Dynamics part of the scheme does not play a significant role in the resolution. The results obtained with the two schemes are very close and all the values of β can be used. The results obtained with β = 0.2 at time T = 0.1 are reported in Figure 2.

We then choose smaller ( = √ 0.1) and for β = 0.2 the fully explicit scheme returns a negative density (which means that the CFL condition for this scheme is violated and the time step should be reduced) whereas the AP scheme still returns a relevant result. A time step divided by 2 (that is β = 0.1) ensures the positivity of the density for both schemes. The results obtained with β = 0.2 (for the AP scheme) and β = 0.1 at time T = 0.05 are reported in Figure 3.

Finally we choose close to zero ( = √ 0.001) and the explicit scheme returns a negative density for β = 0.2 and β = 0.1. We recover a relevant result with β = 0.01. The AP scheme provides a relevant results for all the values of β. The results obtained with β = 0.2, 0.1 (for the AP scheme) and β = 0.01 at time T = 0.007 are reported in Figure 4. This test confirms that the stability condition for the explicit method is sensitive to the value of the Mach number and it shows that the AP strategy allows us to keep significantly larger time steps. It leads to significant computational cost savings (especially when 1), but obviously at the price of a loss of accuracy.

Low Mach regime simulations in 2D

We discuss several test cases in 2D, all with periodic boundary conditions. We adopt the following notation. For discrete quantities ρ, u or v, depending on their location, we introduce the L 1 and L ∞ discrete norms

ρ 1 = i,j δx i+1/2 δy j+1/2 |ρ i+1/2,j+1/2 |, ρ ∞ = max i,j |ρ i+1/2,j+1/2 |, u 1 = i,j δx i δy j+1/2 |u i,j+1/2 |, u ∞ = max i,j |u i,j+1/2 |, v 1 = i,j δx i+1/2 δy j |v i+1/2,j |, v ∞ = max i,j |v i+1/2,j |.
To evaluate the behaviour of the scheme when the mesh size goes to zero, associated to a quantity q (which can stand for ρ, u or v), we consider the following errors in

L 1 or L ∞ discrete norm e (N )
1,q = q (N ) -q exact 1 and e (N ) ∞,q = q (N ) -q exact ∞ ,

where q (N ) stands for the approximate solution obtained on a N × N grid, whereas q exact is the exact solution. Similarly, to evaluate the behaviour of the scheme when goes to zero, we introduce the following quantities e ( ) 1,q = q ( ) -q (0) 1 and e ( ) ∞,q = q ( ) -q (0) ∞ , where q ( ) stands for the approximate solution obtained for a fixed value of , whereas q (0) is the exact solution of the limit problem when → 0. In the computation of the norms, the expression of q exact and q (0) are simply evaluated at the corresponding grid points.

Test 1: traveling vortex

The traveling vortex test case is extracted from [START_REF] Kaiser | A new stable splitting for the isentropic Euler equations[END_REF]. The exact solution is known and a direct comparison is possible. We consider the domain [0, 1] × [0, 1]. The pressure law is given by p(ρ) = 1 2 ρ 2 and the initial density and velocity are given by

ρ(x, y, 0) = 110 + 2 (4π) 2 f (r c ), u(x, y, 0) = v 0 v 1 + g(r c ) 0.5 -y x -0.5 ,
where r c = 4π (x -0.5) 2 + (y -0.5) 2 1 2 and f , g are two functions such that f is continuous and piecewise C 1 , g is continuous and

f (r) = rg(r) 2 , for almost all r.
This problem reduces to a pure transport, the vortex is traveling at speed (v 0 , v 1 ): the exact solution is simply

ρ(x, y, t) = ρ(x -v 0 t, y -v 1 t, 0), u(x, y, t) = u(x -v 0 t, y -v 1 t, 0).
As in [START_REF] Kaiser | A new stable splitting for the isentropic Euler equations[END_REF], we first choose

f (r) = (1.5) 2 δ(r c )(k(r c ) -k(π)), g(r) = 1.5(1 + cos(r c ))δ(r c ), v 0 = 0.6, v 1 = 0,
where

k(r) = 2 cos(r) + 2r sin(r) + 1 8 cos(2r) + 1 4 r sin(2r) + 3 4 r 2 , δ(r) = 1 r<π .
We observe in numerical simulations for this test case some spurious oscillations for large values of (for instance = 0.8 and = 0.5) when we set ι = 0 in the definition (9) of the coefficient a d . We observe that setting ι = 1 is sufficient to eliminate these oscillations so that we adopt this definition for all the simulations performed for this test case. In Figures 5, 6 and7 we show the density, the horizontal and vertical velocity respectively at initial time (left) and final time T = 0.5 (right) for = 0.8, δx = δy = 1/32 and δt = 5 × 10 -4 . Cuts are available in Figure 8. Numerical solutions at time T = 0 (blue solid lines) and T = 0.5 (black dashed lines) are represented. The density (top) and the vertical velocity (bottom right) are plotted along the line y = 0.5 as a function of x -v 0 T . It allows us to compare the shape of the initial solution and the shape of the numerical solution at time T = 0.5. The horizontal velocity (bottom left) is plotted along the line x = 0.5 + v 0 T as a function of y. The results are oscillations-free but we observe an important numerical diffusion inherent to the use of a first order scheme. At T = 0.5, min i,j v i+ 1 2 ,j = -0.069, max i,j v i+ 1 2 ,j = 0.064. We then provide a convergence study in Fig 9 . For different Mach numbers ∈ {0.01, 0.1, 0.4, 0.8}, we perform computations with δt = 5 × 10 -4 on different grids (10 × 10, 14 × 14, 20 × 20, 28 × 28, 40 × 40 and 56 × 56). We plot in logarithmic scales the discrete L 1 error norm at time T = 0.01 between the exact solution and the discrete one (for the density on the left and for the first component of the velocity on the right) as a function of the mesh size. We clearly observe a first order convergence. Similar results are obtained for the second component of the velocity. Finally, to justify the AP character of the scheme, using a 32 × 32 grid and a time step δt = 5 × 10 -4 , we study for different values of the Mach number ∈ {0.0125, 0.025, 0.05, 0.1, 0.2, 0.4, 0.8}, the discrete L 1 and L ∞ -norm at time T = 0.01 of the difference between the discrete solution and ρ (0) = 110. The results are reported in Table 1. It confirms that the density tends to be constant, equal to ρ (0) when tends to zero and that the norm of the difference behaves like c 2 (with c a constant independent of ).

∞,ρ = ρ ( ) -ρ (0)
∞ at T = 0.01 for different on a 32 × 32 grid and δt = 5 × 10 -4 .

Next, concerning the velocity, we provide in Table 2, for different values of ∈ {0.0125, 0.025, 0.05, 0.1, 0.2, 0.4, 0.8} and at T = 0.01, the L 1 norm of the discrete divergence and the L 1 norm of the difference between the approximate velocity u ( ) and the exact limit velocity u (0) obtained when → 0. Unfortunately, this test case does not illustrate the fact that the approximate velocity converges towards the correct solution when → 0. Indeed, in Table 2, the difference in L 1 norm seems to remain close to 2 × 10 -4 when goes to 0. This is probably due to the fact that, in this test case, contrarily to the density, the exact velocity does not depend on or, in other words, this is due to the fact that, even for = 0, the exact velocity is equal to the limit velocity. Therefore, the observed difference is certainly a truncation error due to the mesh rather than a difference due to a positive value of . Nevertheless, even if this difference does not vanish when → 0, we observe that the discrete divergence goes to 0, in accordance with Remark 2.8. For small values of , we recover that the norm of the discrete divergence behaves like c 2 as expected. 1,u = u ( ) -u (0) 1 at T = 0.01 for different on a 32 × 32 grid and δt = 5 × 10 -4 .

In Table 3, we present the results of the same study on a finer mesh (a 64 × 64 grid). The behaviour of the discrete divergence is similar and, as previously, we also observe that the difference in L 1 norm between the approximate velocity u ( ) and the exact limit velocity u (0) does not vanish. However it now remains close to 9.30 × 10 -5 . It confirms the hypothesis that this almost constant error could be due to the discretization error since it becomes (approximatively) twice smaller on a twice finer grid. To observe the behaviour of the discrete velocity when goes to zero, we propose to study a perturbed version of this test case, more relevant for our purpose, by setting

f (r) = (1.5 + ) 2 δ(r c )(k(r c ) -k(π)), g(r) = (1.5 + )(1 + cos(r c ))δ(r c ).
The results are reported in Table 4 (for the density) and in Table 5 (for the velocity and the divergence of the velocity). The behaviour of ρ ( ) and ∇ d • u ( ) remains unchanged. But we now clearly observe that the discrete velocity tends to the exact limit velocity when tends to zero and, as expected, that the norm of the difference behaves like c (with c a constant independent of ). 1,u = u ( ) -u (0) 1 at T = 0.01 for different on a 32 × 32 grid and δt = 5 × 10 -4 .

Test 2: cylindrical explosion problem

We consider a 2D isentropic cylindrical explosion problem extracted from [START_REF] Dimarco | Study of a new asymptotic preserving scheme for the Euler system in the low Mach number limit[END_REF]. The computational domain is Ω = [-1, 1] × [-1, 1], discretized with a 50 × 50 grid. The pressure is given by p(ρ) = ρ and the initial density and velocity are given by ρ(x, y, 0)

= 1 + 2 1 r(x,y) 1/4 , u(x, y, 0) = -xβ(x, y) r(x, y)ρ(x, y, 0) 1 r(x,y)>10 -15 , v(x, y, 0) = -yβ(x, y) r(x, y)ρ(x, y, 0) 1 r(x,y)>10 -15 ,
where r is the distance to the center of the domain, r(x, y) = x 2 + y 2 and β(x, y) = max{0, 1 -r(x, y)} × 1 -e -16r(x,y) 2 .

The initial density and velocity field are shown in Fig. 10. Like in [START_REF] Dimarco | Study of a new asymptotic preserving scheme for the Euler system in the low Mach number limit[END_REF] we display the numerical solution obtained at time t = 0.1, t = 0.25 and t = 0.5 for = 1 in Fig. 11 to Fig. 13. The simulation is performed on a 50 × 50 grid with a time step δt = 5 × 10 -4 .

The results obtained for = 1 at time t = 0.1 and t = 0.5 are very close to the results obtained with a different AP-scheme in [START_REF] Dimarco | Study of a new asymptotic preserving scheme for the Euler system in the low Mach number limit[END_REF]. Note that, as confirmed by its authors, the result reported in [START_REF] Dimarco | Study of a new asymptotic preserving scheme for the Euler system in the low Mach number limit[END_REF] for t = 0.25 does not correspond to this simulation time. We warmly thanks the authors for taking the time to re-consider this test case and for many useful hints. We present in Fig. 14 the results obtained at t = 0.05 using the AP-scheme with = 5×10 -3 . We use again a 50 × 50 grid and compare the velocity fields obtained with a small time step, δt = 10 -4 , at left and a large one, δt = 5 × 10 -3 , at right. As in [START_REF] Dimarco | Study of a new asymptotic preserving scheme for the Euler system in the low Mach number limit[END_REF], we observe that the AP scheme resolves the small waves dynamics when using a small time step but also that the AP-scheme can handle large time step and, in this case, it captures the limit solution.

Test 3

The third test case is extracted from [START_REF] Degond | All speed scheme for the low Mach number limit of the isentropic Euler equations[END_REF]. The computational domain is Ω = [0, 1] × [0, 1]. The final time is T = 1 and the time step is given by δt = 2.5 × 10 -2 . The pressure law is given by p(ρ) = ρ 2 and the initial data reads ρ = 1 + 2 sin 2 (2π(x + y)), q u = sin(2π(x -y)) + 2 sin(2π(x + y)) and u = q u ρ , q v = sin(2π(x -y)) + 2 cos(2π(x + y)) and v = q v ρ .

As in [START_REF] Degond | All speed scheme for the low Mach number limit of the isentropic Euler equations[END_REF], we compute the solution for = 0.8 and = 0.05. For = 0.8, the result in [START_REF] Degond | All speed scheme for the low Mach number limit of the isentropic Euler equations[END_REF] is presented on a 20 × 20 grid; we obtain a similar figure with a 12 × 12 grid, but the results are different when the mesh is refined. With a 20 × 20 grid, the result given by the AP method is close to the one obtained with a 64 × 64 grid. Results for = 0.8 using a 64 × 64 grid are presented in Figure 15. The result obtained with = 0.05 using a 64 × 64 grid are presented in Figure 16 for the AP-Scheme and in Figure 17 with the limit scheme [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF]. The density obtained with the AP-scheme is almost equal to the constant state ρ (0) = 1: ρ -ρ (0) ∞ ∼ 10 -3 , the horizontal and vertical velocities have converged to the solution of the limit scheme.

Incompressible Euler simulations using the limit scheme

We now turn to the resolution of the incompressible Euler system in order to check the validity of the limit scheme [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF]. We consider two test cases from [START_REF] Shu | A numerical resolution study of hight order eseentially non-oscillatory schemes applied to incompressible flow[END_REF]. The computational domain is Ω = [0, 2π] × [0, 2π] for both test cases. The density is assumed to be constant equal to 1.

Test A: analytical solution

We test the capability of the limit scheme [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF] to (approximately) preserve a stationary solution of the incompressible Euler system: u(x, y, t) = -sin(x) cos(y), v(x, y, t) = cos(x) sin(y) and p(x, y, t) = 1 2 (cos(x) 2 + cos(y) 2 ).

The exact vorticity equal to w(x, y, t) = -2 sin(x) sin(y).

The initial data is chosen as the evaluation on the grid of the stationary solution and using the limit scheme (20), we perform different numerical simulations up to time T = 2 on different meshes with a time step δt = 2 × 10 -4 . The result obtained using a 160 × 160 mesh is reported in Figure 18 (vorticity) and in Figure 19 (velocity and pressure), the discrete vorticity w being defined by

w i,j = v i+ 1 2 ,j -v i-1 2 ,j δx i - u i,j+ 1 2
-u i,j-1 2 δy j .

We next provide in Table 6 the L ∞ error norms e ∞,w , N ∈ {40, 60, 80, 100, 120, 140, 160}, for the horizontal velocity u (the results for the vertical velocity v are similar) and for the vorticity w. As expected, we observe a first order convergence. ∞,u and e (N ) ∞,w using the limit scheme [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF] for different meshes at time T = 2 with δt = 2 × 10 -4

Test B: double shear layer

The initial data of the last test case are given by u(x, y, 0) = tanh y -π 2 π/15 1 y π + tanh The exact solution is not known but we can compare with the results given in [START_REF] Shu | A numerical resolution study of hight order eseentially non-oscillatory schemes applied to incompressible flow[END_REF], based on computations with the ENO scheme on 64 × 64 and 128 × 128 grids. We perform a simulation grid with δt = 10 -3 using the limit scheme (20).

using the limit scheme (20) on a 128 × 128 grids. We display the vorticity in Figure 20 with 60 equally spaced contours from -4.9 to 4.9 at t = 4 (left up), t = 6 (right up), t = 8 (left bottom) and t = 10 (right bottom). These results are qualitatively very close to the ones obtained with the ENO scheme on the 64 × 64 grid in [START_REF] Shu | A numerical resolution study of hight order eseentially non-oscillatory schemes applied to incompressible flow[END_REF]. The ENO scheme used in [START_REF] Shu | A numerical resolution study of hight order eseentially non-oscillatory schemes applied to incompressible flow[END_REF] is fourth-order accurate (in the L 1 sense), while our scheme is only first order accurate. Nevertheless, the vortices developed by the solution with smaller and smaller scales are well reproduced.

  where α = 2 . We first remark that in this case c(0) (ρ, ν) = |ν (0) |.

Figure 2 :

 2 Figure 2: 1D Riemann problem, comparison between the AP scheme and an explicit scheme for ( , T ) = ( √ 0.99, 0.1) and δt = βδx.

Figure 3 :

 3 Figure 3: 1D Riemann problem, comparison between the AP scheme and an explicit scheme for ( , T ) = ( √ 0.1, 0.05) and δt = βδx.

Figure 4 :

 4 Figure 4: 1D Riemann problem, comparison between the AP scheme and an explicit scheme for ( , T ) = ( √ 0.001, 0.007) and δt = βδx.

Figure 5 : 1 2 ,j+ 1 2 =Figure 6 : 2 = 0 2 =Figure 7 :

 51262027 Figure 5: Test 1 (traveling vortex), density at time T = 0 (left) and T = 0.5 (right). At T = 0.5, min i,j ρ i+ 1 2 ,j+ 1 2

Figure 8 :

 8 Figure 8: Test 1 (traveling vortex), translated cuts at time T = 0 (blue solid lines) and T = 0.5 (black dashed lines). Density along y = 0.5 as a fonction of x -v 0 T (top). Horizontal velocity along x = 0.5+v 0 T as a function of y (bottom left). Vertical velocity along y = 0.5 as a fonction of x -v 0 T (bottom right).

Figure 9 :

 9 Figure 9: Test 1 (traveling vortex) L 1 error norm at T = 0.01 as a function of the mesh size for different values of and δt = 5 × 10 -4 .

  e

  e

  = u ( ) -u (0) 1 at T = 0.01 for several on a 64 × 64 grid and δt = 5 × 10 -4 .

Figure 10 :

 10 Figure 10: Test 2 (explosion problem), time t = 0 and = 1, velocity field (left) and density profile (middle and right).

Figure 11 :

 11 Figure 11: Test 2 (explosion problem), time t = 0.1 and = 1, velocity field (left) and density profile (middle and right).

Figure 12 :

 12 Figure 12: Test 2 (explosion problem), time t = 0.25 and = 1, velocity field (left) and density profile (middle and right).

Figure 13 :

 13 Figure 13: Test 2 (explosion problem), time t = 0.5 and = 1, velocity field (left) and density profile (middle and right).

  (a) δt = 10 -4 (b) δt = 5 × 10 -3

Figure 14 :

 14 Figure 14: Test 2 (explosion problem), velocity field at time t = 0.05, using the AP-scheme with = 5 × 10 -3 , comparison between the results obtained for δt = 10 -4 (left) and δt = 5 × 10 -3 (right).

Figure 15 :

 15 Figure 15: Test 3, = 0.8, horizontal and vertical velocity (left and middle) and density (right).

Figure 16 :

 16 Figure 16: Test 3, = 0.05, horizontal (left) and vertical velocity (right).

Figure 17 :

 17 Figure 17: Test 3, horizontal (left) and vertical velocity (right) using the limit scheme (20).

Figure 18 :Figure 19 :

 1819 Figure 18: Test A, vorticity w using the limit scheme (20) at T = 2 on a 160 × 160 grid with δt = 2 × 10 -4 .

3π 2 -

 2 y π/15 1 y>π , v(x, y, 0) = 0.05 sin(x).

Figure 20 :

 20 Figure 20: Test B (double shear layer), contour of the vorticity w at different times on a 128×128grid with δt = 10 -3 using the limit scheme[START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF].

Table 1 :

 1 Test 1 (traveling vortex) Errors e

			( ) 1,ρ	log e log( +1 / ) ( +1 ) 1,ρ /e ( ) 1,ρ	( ) e ∞,ρ	log e log( +1 / ) ( +1 ) ∞ ( ) /e ∞
	1	0.8	1.06 × 10 -3	1.62	2.54 × 10 -2	1.84
	2	0.4	3.44 × 10 -4	2.24	7.11 × 10 -3	2.00
	3	0.2	7.30 × 10 -5	2.32	1.78 × 10 -3	2.00
	4	0.1	1.46 × 10 -5	2.00	4.46 × 10 -4	1.99
	5	0.05	3.64 × 10 -6	2.00	1.12 × 10 -4	2.01
	6 0.025 9.13 × 10 -7	2.00	2.79 × 10 -5	2.00
	7 0.0125 2.28 × 10 -7	--	6.97 × 10 -6	--
				( ) 1,ρ = ρ ( ) -ρ (0)	1 and e	( )

Table 2 :

 2 Test 1 (traveling vortex) Discrete L 1 norm of the discrete divergence e

	( ) 1,div =

Table 3 :

 3 Test 1 (traveling vortex) Discrete L 1 norm of the discrete divergence e

Table 4 :

 4 Test 1 (traveling vortex, perturbed version) Discrete L 1 norm of e

	( +1 ) 1,ρ	/e ( ) 1,ρ )
	log( +1 / )

( ) 1,ρ = ρ ( ) -ρ (0) 1 at T = 0.01 for different on a 32 × 32 grid and δt = 5 × 10 -4 .

Table 5 :

 5 Test 1 (traveling vortex, perturbed version) Discrete L 1 norm of the discrete divergence e

	( ) 1,div = ∇ d • u ( )	1 and e ( )

Table 6 :

 6 Test A, error norms e

	1,w )
	log(N /N +1 )