domain: phys sdv

text:
Daniele Funaro
email: daniele.funaro@unimore.it

A Dynamic Representation of mRNA Nucleotides Clarifies the Conundrum of Codon Redundancy

Keywords: RNA, codon degeneracy, electric circuit, biophoton

The deciphering of the genetic code takes place through the reading of the nitrogenous bases, which are 4 in number. In most cases, the bases are taken three by three, thus generating 64 possible combinations with repetition. Each combination (codon) allows for the synthesis of a specific amino acid. Since the latter are only 21 in number, the codon-amino acid conversion table shows a strong redundancy. Countless efforts have been made to understand the true encryption mechanism. Here we want to add our version, which consists of associating a periodic sound based on 3 notes to each codon. RNA now becomes a dynamic object and not just a list of static instructions. In addition to a different interpretation of the genetic code, there is also a considerable reduction in redundancy, given that the number of periodic sounds that can be produced with 3 notes drops to 20 (with the addition of 4 pure frequencies). Finally we discuss the possibility of how these sounds can be generated and travel inside the double helix, and possibly emitted as biophotons.

Introduction

It is well known that RNA is a string of nitrogenous bases containing the sequence of instructions for constructing amino acids, as well as other interesting information. The organization suggests a strictly deterministic procedure. Once having established the beginning of a sequence, the nucleotides are considered three by three, thus forming a series of codons. Each codon is therefore identified by a triplet of symbols, taken from the four available, i.e.: U (Uracil), C (Cytosine), A (Adenine), G (Guanine). There are 64 possible combinations. Each of them allows us to identify a single amino acid according to the famous 1 conversion chart (see table 1), obtained thanks to the pioneering efforts of M. Nirenberg and coworkers [START_REF] Nirenberg | The genetic code[END_REF].

CUU → Leu CCU → Pro CAU → His CGU → Arg CUC → Leu CCC → Pro CAC → His CGC → Arg CUA → Leu CCA → Pro CAA → Gin CGA → Arg CUG → Leu CCG → Pro CAG → Gin CGG → Arg AUU → Ile ACU → Thr AAU → Asn AGU → Ser AUC → Ile ACC → Thr AAC → Asn AGC → Ser AUA → Ile ACA → Thr AAA → Lys AGA → Arg AUG → Met ACG → Thr AAG → Lys AGG → Arg GUU → Val GCU → Ala GAU → Asp GGU → Gly GUC → Val GCC → Ala GAC → Asp GGC → Gly GUA → Val GCA → Ala GAA → Glu GGA → Gly GUG → Val GCG → Ala GAG → Glu GGG → Gly
The result is extraordinary but not perfectly convincing from the theoretical viewpoint. In fact, there are only 21 amino acids so there are several codons that lead to the same result. The table therefore displays some redundancy (called codon degeneracy) and there does not seem to exist a clear law that allows the combinations to be thinned out in a logical way. In particular, there are 8 amino acids having 4 or more codon representations (Ser, Leu, Pro, Arg, Thr, Val, Ala, Gly) and 12 amino acids having 1 to 3 codon representations (Phe, Leu, Tyr, Cys, Trp, His, Gin, Ile, Asn, Lys, Asp, Glu). Finally, there are 4 special codons. One of them corresponds to the amino acid Met, that determines the beginning of a protein. Each of the last 3 codons represent the command stop that marks instead the end of a protein. Many attempts have been made in the past to try to better formalize the table (see, e.g. [START_REF] Grosjean | An integrated, structure-and energy-based view of the genetic code[END_REF]). Given that nature makes its resources available in a very parsimonious way, the aim is to understand if there are hidden relationships that justify the existence of 64 combinations when in the end only a third would be necessary.

The idea we propose here is to interpret the genetic code in a dynamic way. RNA is not trivially a score to be read and performed by other, more or less complex, molecules. Instead, it is itself the entire musical instrument, which decrypts and at the same time plays appropriate chords associated with each codon. The trick is to assume that the message carried by a codon can be repeated periodically so as to become a sort of sound. Different amino acids will correspond to different sounds, but there will be however less sounds than codons. Let's do some calculations.

Suppose to have a list of four symbols Γ, ∆, Ξ, Ω. We can consider all possible triplets that can be formed using the above characters. We know that all combinations amount to 64. We then take the periodic signal generated by the repetition of a given triplet for an indefinite number of times. For example, if we choose ΩΓΞ, the corresponding "sound" is

• • • ΩΓΞΩΓΞΩΓΞΩΓΞΩΓΞΩ • • • .
The question is: how many different sequences can be built? First of all, there exist 4 triplets ΓΓΓ, ∆∆∆, ΞΞΞ, ΩΩΩ, that generate trivial constant monochromatic sequences. We are left with other 60 more interesting cases. Let us now observe that the sequence originating for example from ΩΓΞ coincides with those generated by ΓΞΩ or ΞΩΓ. Indeed, the above three combinations give rise to the same sound, therefore they are redundant. This means that it is possible to group within a single class the 3 triplets that produce a specific sound. Dividing 60 by 3, we get 20 classes, i.e. 20 different possible sounds. Incidentally, with the exception of Met, 20 is also the number of amino acids. The coincidence is somewhat surprising. But we haven't finished yet. There are 8 sounds where the symbols are all different. It may not be a just coincidence that exactly 8 amino acids admit representations with 4 or 6 redundant codons. There are 12 sequences where a given symbol appears twice, such as for example: • • • ∆∆Γ∆∆Γ∆∆Γ∆∆Γ∆∆ • • • . Do they correspond to the other 12 amino acids? There are finally 4 monochromatic sounds. Do they correspond to the start and stop signals?

In section 2, we will show how to generate 20 + 4 different sounds starting from certain assigned frequencies. Let's say right away that we will not be able to fully solve the problem of rewriting with a minimal number of entries the conversion chart of Table 1. In practice we will not be able to associate in automatic way the generic symbols Γ, ∆, Ξ, Ω with the classic ones representing the nitrogenous bases, i.e.: U, C, A, G. Indeed, it will be soon evident that a straight relation does not hold. We think that, physically, a set of three ordered codons constitutes a unique oscillating electric circuit which elaborate three resonating frequencies (described by some of the symbols Γ, ∆, Ξ, Ω) which in turn produce a characteristic sound. We do not yet know exactly how these circuits work. We provide some ideas in section 3. The purpose of this paper is mainly to interpret the genetic code encrypted in the DNA in an active way, which is realized through the emission of electromagnetic messages at a distance. For the reasons that will be better expressed in section 3, we believe that DNA is not only a data container but also transmits its data into the environment. At the moment this conjecture allows us to understand why the symbols U, C, A, G produce redundant combinations, while this is no longer verified if we use symbols Γ, ∆, Ξ, Ω. The connection between the two sets of symbols is yet to be found. Note that we are not just sidestepping the problem with semantic speculation, we are offering a new way to understand how DNA works.

Sound generation

We propose the following values for our parameters:

Γ = 1 ∆ = 2 Ξ = 4 Ω = 5 (1)
Such a choice is arbitrary, meaning that other choices can work in a similar way. We denote by (ω a , ω b , ω c ) the generic triplet of values obtained by choosing randomly (with possible repetition) three members from the set {Γ, ∆, Ξ, Ω}.

According to what has been said in section 1, we want the triplet (ω a , ω b , ω c ) to be equivalent to (ω c , ω a , ω b ) and (ω b , ω c , ω a ) (in fact, they all generate the same sound). Instead, (ω a , ω b , ω c ) will not be equivalent in general to (ω a , ω c , ω b ), corresponding to another sound in our construction. Successively, we set δ = 2π/3 and define the functions of the variable x:

f 1,a (x) = sin[ω a (x -δ)] sin(ω a x) sin(2ω a δ) sin(ω a δ) (2)
f 2,b (x) = - sin[ω b (x -δ)] sin[ω b (x + δ)] [sin(ω b δ)] 2 (3)
f 3,c (x) = sin[ω c (x + δ)] sin(ω c x) sin(2ω c δ) sin(ω c δ) (4)
f 4,a,b,c (x) = 1 2 sin[ω a (x + δ)] sin(ω b x) sin[ω c (x -δ)] (5)
f a,b,c (x) = f 1,a (x) + f 2,b (x) + f 3,c (x) + f 4,a,b,c (x) (6)
The above are periodic functions with period equal to 3δ = 2π. For a fixed a, we have the relations, ∀x:

f 1,a (x) + f 2,a (x) + f 3,a (x) = 1 f a,a,a (x) = 1 -1 8 sin(3ω a x) (7)
f 1,a (x + δ) = f 3,a (x) f 1,a (x -δ) = f 2,a (x) (8)
We finally have:

f a,b,c (x + δ) = f b,c,a (x) = f c,a,b (x -δ) (9)
This last property says that an even permutation of the indices produces the same wave, up to a phase change. According to the second equation in [START_REF] Bixon | Long-range and very long-range charge transport in DNA[END_REF], when all indices are equal to a we get a monochromatic wave of frequency ω a /δ. Otherwise, using (1), we obtain 20 different sounds which are plotted in Figs. 1 and2. As we said, the idea is to associate each sound with an amino acid, and vice versa. In this way there is no more redundancy. The problem remains of how a sound can be generated from a single codon, respecting the rules of the conversion table 1.

DNA circuitry

Electric circuits displaying a phase shift equal to 2π/3 can be easily built (see Fig. 3). Signals such as those studied in the previous section can be generated by collecting the outputs of three devices of this type resonating at different frequencies. Each circuit has its own RC network and the sum of the frequencies is done by imposing suitable shifting parameters so to recover periodic composite waves.

Figure 3: Typical oscillating electronic circuit with a 2π/3 phase shift.

It has been ascertained under various circumstances that the DNA double helix can be associated with some electric circuit. The results are both experimental and computational. They confirm for instance the conductivity of the periodic phosphate backbone [START_REF] Brillouin | Giant molecules and semiconductors[END_REF][START_REF] Ikeura-Sekiguchi | Attosecond electron delocalization in the conduction band through the phosphate backbone of genomic DNA[END_REF], thanks to its one-dimensional periodic structure. Electronic properties of DNA and its nucleotides have been studied for instance in [START_REF] Bakhshi | On the electronic structure and conduction properties of aperiodic DNA and proteins, II: Electronic structure of aperiodic DNA[END_REF][START_REF] Lewis | Direct measurement of hole transport dynamics in DNA[END_REF][START_REF] Bixon | Long-range and very long-range charge transport in DNA[END_REF][START_REF] Giese | Electron transfer in DNA[END_REF]. We observe that the presence of nitrogen can give compounds a tendency to flip from one configuration to another, activating spontaneous oscillations. Concatenations of logic gates using basic molecules as switches have been taken into consideration [START_REF] Andréasson | Molecules with a sense of logic: a progress report[END_REF][START_REF] Zhou | DNA-based constitutional dynamic networks as functional modules for logic gates and computing circuit operations[END_REF]. A model of DNA composed by multiple oscillatory RLC circuits is esamined in [START_REF] Xu | Stepwise oscillatory circuits of a DNA molecule[END_REF]. These micro circuits are joined by the phosphate backbones that collect the data. All these results reinforce the idea that DNA should be interpreted in a dynamic way. Therefore, together with the intrinsically static message written in the codons, we also must consider a high probability that this should be assigned a periodic sound, which brings us back to the content of this paper.

We have no element at this time to continue our investigation. A way should be proposed to assign to each nitrogenous base an electric circuit. The global circuit associated with each codon should then produce a periodic sound (among the 20+4 possible). In principle this is not difficult. However, it is difficult to respect the hidden rules of the codon table. For example, the generic sequence GGX (where X is any of the symbols U, C, A, G) corresponds to Glycine (last box in Table 1). For some unknown reason, the electric circuit associated to this codon resonates with the couple GG and inhibits the presence of the third component. We undertake to study this issue further in a future paper. Meanwhile let's add some other considerations in favor of a vision in which DNA is an instrument actually played, and not just a static collection of data.

The analysis of vibration frequency maps carried out for nucleobases [START_REF] Jiang | Development of vibrational frequency maps for nucleobases[END_REF][START_REF] Prokhorenko | New insights into the photophysics of DNA nucleobases[END_REF] shows strong DNA activity in producing or absorbing biophotons. Such an intense electromagnetic turmoil regulates the most fundamental biological processes. This aspect strengthens the hypothesis of signals originating from DNA, which remotely guide the main functions of a single cell, but which can extend further, transferring information between groups of cells [START_REF] Morales | Photon dissipation as the origin of information encoding in RNA and DNA[END_REF][START_REF] Cifra | Electromagnetic cellular interactions[END_REF][START_REF] Scholkmann | Non-chemical and non-contact cell-to-cell communication: a short review[END_REF]. These aspects have also been taken into consideration in [START_REF] Funaro | From Photons to Atoms: The Electromagnetic Nature of Matter[END_REF] (section 3.7) and [START_REF] Funaro | The impact of a pervasive electrodynamical background on biological interactions[END_REF] in a more general context, where the chemical elements themselves arise from an agglomeration of atomic nuclei and electromagnetic manifestations that connect them in a peculiar way characterizing the compound.

Finally, we must not forget that DNA is a twisted strand that is difficult to access from a physical viewpoint. Moreover, the cell is a system overcrowded with molecules [START_REF] Zimmerman | Macromolecular crowding: biochemical, biophysical, and physiological consequences[END_REF][START_REF] Ellis | Macromolecular crowding: obvious but underappreciated[END_REF], in which it is practically impossible to create connections except through communications that take place with the exchange of electromagnetic impulses. Of course, it is necessary to assume that suitable receiving stations exist. The RNA polymerase could be one of them, thus eliminating direct contact with DNA, which appears unrealistic in in vivo situations.

Conclusions

Recent research shows more and more that living organisms are an extended set of molecules (variably complex) that communicate through the exchange of electromagnetic messages. DNA itself is a source of biophotons. In this context it is natural to hypothesize that the genetic code should mainly be interpreted in a dynamic way, which occurs for example through the emission of periodic signals (sounds in our context). The conjecture is strengthened if one considers that the periodic signals obtainable from three (out of four) independent frequencies are 20+4 in number. This suggests that each codon is a transmitting machine, which assigns each sound the creation of a specific amino acid. Although from a circuital point of view the codon still remains to be investigated, the representation through sounds allows us to take a step forward in understanding why codon degeneracy occurs in the static version of

4 2Figure 1 :

 41 Figure 1: Sounds produced by the triplets specified at the top of each plot.

Figure 2 :

 2 Figure 2: Sounds produced by the triplets specified at the top of each plot.

6

 6

Table 1 :

 1 Codon table

	UUU → Phe	UCU → Ser	UAU → Tyr	UGU → Cys
	UUC → Phe	UCC → Ser	UAC → Tyr	UGC → Cys
	UUA → Leu	UCA → Ser	UAA → Stop	UGA → Stop
	UUG → Leu	UCG → Ser	UAG → Stop	UGG → Trp

DNA, but not in the dynamic one.


Evelyne Akopyan Tésa

Riadh Dhaou
email: riadh.dhaou@toulouse-inp.fr

Emmanuel Lochin
email: emmanuel.lochin@enac.fr

Bernard Pontet
email: bernard.pontet@cnes.fr

Jacques Sombrin
email: jacques.sombrin@tesa.prd.fr

On the Network Characterization of Nano-Satellite Swarms

Keywords: Swarm of Nano-Satellites, Network Characterization, Graph Theory

Low-frequency radio interferometry is crucial to understanding the universe and its very early days. Unfortunately, most of the current instruments are ground-based and thus impacted by the interferences massively produced by the Earth. To alleviate this issue, scientific missions aim at using Moonorbiting nano-satellite swarms as distributed radio-telescopes in outer space, keeping them out of Earth interference range. However, swarms of nano-satellites are systems with complex dynamics and need to be appropriately characterized to achieve their scientific mission. This paper presents a methodology based on graph theory for characterizing the swarm network system by computing graph theory metrics around three properties: the node density, network connectivity and ISL availability. We show that these properties are well-suited for highlighting a possible heterogeneity in a network and adapt a routing strategy accordingly. This work is the first milestone in defining the best-suited routing strategy within the swarm from the derived network properties.

I. INTRODUCTION

The study of the low-frequency range (below 100 MHz) is essential for many scientific fields, such as astrophysics for sky mapping and monitoring or the observation of the Dark Ages signals, which are signatures of the very early universe [START_REF] Brown | Non Thermal Galactic Emission Below 10MHz[END_REF], [START_REF] Manning | The Galactic background radiation from 0.2 to 13.8 MHz[END_REF]. Until now, the majority of low-frequency radio interferometry instruments are ground-based and provide high-quality observations of space. However, the capabilities of these instruments are limited by ionospheric distortions, terrestrial Radio Frequency Interferences (RFI), as well as a complete reflection of radio waves below 10-30 MHz [START_REF] Bentum | A novel astronomical application for formation flying small satellites[END_REF]. Therefore, the cosmic radio signals that are weak in terms of power can be easily altered if not completely masked.

One solution to this problem is to create an interferometer directly in space: the Nano-satellites for a Radio Interferometer Observatory in Space (NOIRE) [START_REF] Cecconi | NOIRE study report: towards a low frequency radio interferometer in space[END_REF] study investigates and proves the feasibility of using a swarm of nano-satellites for low frequency radio observation in space. Indeed, a swarm of approximately 100 nano-satellites orbiting the Moon is highly protected from the Earth's RFI, and thus appears as a promising solution to the interference problem met by groundbased telescopes.

The setup and configuration of a swarm of nano-satellites as a space observatory is a challenging problem in terms of communication, as the network relies solely on wireless Inter-Satellite Links (ISL). It is possible to approximate this system as a Wireless Sensor Network (WSN), with the particularity that the average inter-satellite distance goes up to 50 km (interferometry instrument requirement), making it a very lowdensity network. Besides, the satellites move at an average speed of 1 to 10 km/s, making the swarm a highly mobile system and adding constraints to the communication model.

But most importantly, the system needs to operate in a distributed manner. Indeed, unlike traditional observation satellites that usually carry one or many measuring devices, a swarm of nano-satellites consists of a single instrument distributed over many satellites. Thus, each satellite collects observation data from space, then shares them among the swarm to compute one single image of cosmic rays. This image is a matrix of cross-correlations between all data collected by every satellite and shared among the swarm.

The main challenge is thus to disseminate significant amounts of data within the swarm, which raises major communication issues. First, the simultaneous propagation of several gigabytes of data over a radio channel by each satellite can lead to potential link congestion, delay, and packet loss. Then, the energy consumption of the swarm is proportional to the amount of data shared between the satellites. Therefore, a trade-off has to be made to limit the total energy consumption of the swarm to alleviate these issues and prevent too early energy depletion of the satellites. It becomes primordial to correctly measure the system: finding an optimal routing strategy to perform this data collection cannot be conceived without proper network characterization.

In this context, we focus on studying the specificities of a swarm of nano-satellites on the telecommunications level and, more precisely, on the characterization of the network layer to further determine a suitable routing scheme. We present a set of metrics based on graph theory that prove to perfectly grasp the network dynamics and heterogeneous nature when appropriate. The choice of routing strategy is only discussed and oriented at this stage, but will be thoroughly studied in a future work. We first present, in Section II, the state-of-the-art related to routing challenges. Section III describes the nanosatellite swarm model, then the network properties are derived in section IV. Finally, section V summarizes the essential results and perspectives of our work.

II. RELATED WORK

Similar case studies exist in the literature, such as the Orbiting Low-Frequency Array for Radioastronomy (OLFAR) project [START_REF] Budianu | OLFAR: Adaptive topology for satellite swarms[END_REF]. Their objective is to use a swarm of at least 50 nano-satellites as a distributed radio-telescope orbiting the Moon. In this paper, the authors focus on the relevance of dynamic node clustering to minimize the energy consumption related to data sharing. The authors demonstrate that classic clustering algorithms (such as global k-means or MOBIC) are not well-suited for highly mobile Wireless Sensor Networks (WSN). Thus, they provide a slave-master algorithm that proves to be a better solution for this topology. Regarding the same OLFAR project, in [START_REF] Budianu | Swarm-to-Earth communication in OLFAR[END_REF], the authors emphasize the critical nature of the downlink communication from the swarm to the Earth and present a resilient transmission mechanism.

To this day, routing algorithms are better studied on satellite constellations or mobile WSN. In [START_REF] Raveneau | A multi-level FREAK DTN: Taking care of disconnected nodes in the IoT[END_REF], the authors propose an IoT-compliant routing algorithm, named Multi-Level FREAK, which improves routing performances in a network where nodes can disconnect from their destination. Their solution requires very few memory and computation power, but their case study is different from ours: indeed, our objective is to share data with a large subset of nodes, whereas their objective is to send data to a precise destination. In [START_REF] Diana | DTN routing for quasi-deterministic networks with application to LEO constellations[END_REF], the authors introduce a simple yet efficient network characterization methodology for choosing an adapted routing scheme. They show that in the case of satellite constellations, which are quasi-deterministic mobile networks, an opportunistic routing scheme is better suited as the studied network architecture is closer to a Delay-Tolerant Network (DTN) than a Mobile Ad-Hoc Network (MANET).

However, in the case of a swarm of nano-satellites, abrupt changes in topology can occur. One solution to this problem is proposed in [START_REF] Whitbeck | HYMAD: Hybrid DTN-MANET routing for dense and highly dynamic wireless networks[END_REF], where the authors introduce the HYMAD algorithm. HYMAD is a hybrid DTN-MANET routing scheme that relies on topology information to adapt to the node density of the network, i.e., the algorithm behaves as a MANET routing scheme within highly dense zones, and a DTN scheme between these zones where the density is lower. The authors prove the efficiency of their algorithm in terms of delay and delivery ratio. Another solution is to rely on a routing scheme based on swarm intelligence: in [START_REF] Caro | Swarm intelligence for routing in mobile ad hoc networks[END_REF], the authors introduce an intelligent algorithm inspired by ants colonies, AntHocNet. Their algorithm is adapted for MANET and proves to be more efficient regarding delay jitter and routing overhead than a classic algorithm such as AODV.

Finally, a promising approach for determining a routing scheme for a swarm of nano-satellites is to study Flying Ad-Hoc Networks (FANET), and more specifically, swarms of drones. Indeed, the characteristics of FANET are very close to our study case: the node mobility is very high, the average density is low, and there are fast topology changes. The mobility model is also pre-determined in most cases. In [START_REF] Hyland | Simulation-based performance evaluation of mobile ad hoc routing protocols in a swarm of unmanned aerial vehicles[END_REF], the authors prove that it is possible to exploit the spatial multiplexing of the wireless channel by decreasing the connection range and implementing multi-hop routing. They evaluate the performances of MANET routing protocols such as OLSR and AODV, and Greedy Perimeter Stateless Routing (GPSR) [START_REF] Karp | GPSR: Greedy perimeter stateless routing for wireless networks[END_REF] in a swarm of drones. A survey on FANET communication protocols for physical, MAC, network, and transport layers is also proposed in [START_REF] Bekmezci | Flying ad-hoc networks (FANETs): A survey[END_REF].

III. MODEL DESCRIPTION

We present in this section the objectives and requirements of our study, then describe the modelling of the swarm network as well as the basic connection hypotheses.

A. Objectives and Requirements

The adaptability of swarms of nano-satellites is a strong advantage: such missions can offer a wide spectrum of study opportunities. In addition, a lunar orbit is particularly wellsuited for these missions, as the swarm is not disturbed by interferences coming from the Earth.

At each revolution around the Moon, scientific missions, like those described in the NOIRE study, operate in at least four steps:

1) Observation: each nano-satellite collects raw data from space. 2) Inter-satellite transmission: each nano-satellite broadcasts their data to the others and gathers the data from the other nano-satellites. This step is the focus of the paper. 3) Computation: each nano-satellite computes a part of the overall space image combining the data collected and received from the others. 4) Swarm-to-Earth transmission: the swarm downlinks the computed space image to a base station on Earth. This step is thoroughly studied in [START_REF] Budianu | Swarm-to-Earth communication in OLFAR[END_REF]. The objective of our study is network-based. The main goal is to determine a methodology for characterizing the properties of a swarm topology and thus derive the most suited routing scheme for data propagation within the swarm. Indeed, understanding the network dynamics is essential to efficiently share the data within the swarm while optimizing the energy consumption and the duration of the most interesting phase for the scientific community: the observation of space.

B. Description of the Network

Let us model the swarm network at each timestamp t ∈ T as a graph G t (N t , E t ) with a set of nodes N t , representing the nano-satellites, and a set of edges E t , representing the ISL. We assume the nodes operate in a faultless manner, i.e.:

∀(i, j) ∈ T 2 , |N i | = |N j | = |N |
This condition guarantees that the number of nodes in the graph remains constant with time. In our case, the studied network is a swarm composed of |N | = 100 nano-satellites orbiting the Moon, as depicted in Fig. 1a. The nano-satellites are all identical and use omnidirectional antennas to communicate with each other. The swarm system has the following properties:

1) Mobility: the nano-satellites are mobile within the swarm according to a quasi-deterministic mobility model, i.e., ∃(i, j) ∈ T 2 | E i ̸ = E j . This property highlights the dynamic nature of the ISL within the swarm network as depicted in Fig. 1b and Fig. 1c. The mobility model is more thoroughly studied in section IV.

2) Temporal Drift: the swarm mobility with respect to the Moon is pseudo-periodic, i.e., a small drift over time is randomly induced for each nano-satellite.

These properties imply a temporal evolution of the swarm topology. Intuitively, the routing strategy should be able to adapt to topology changes.

C. Connection Hypotheses

The connectivity within the swarm is exclusively based on ISL. Let e t (u, v) be the edge between two nodes u and v of the graph G at the time t, and d(u, v) the Euclidean distance between these nodes. The set of edges E t is defined as follows:

E t = {e t (u, v) | d(u, v) ≤ R G } ∀ (u, v) ∈ N 2 t (1)
where R G is the connection range of the model: it is independent of time and identical for all nodes at all times. We also assume that each ISL is a duplex link, i.e., e t (u, v) = e t (v, u). Additional connection hypotheses are listed below:

• All data messages are identical in size;

• every message is broadcasted. This hypothesis is based on the presence of the omnidirectional antenna on each nano-satellite; • each node of the network can forward messages in the manner of a router; • the throughput can be adaptive. A trade-off has to be made on the throughput, as a high throughput reduces the transmission duration and thus gives more time for the observation phase. Still, it also increases the energy consumption of the nano-satellite.

IV. CHARACTERIZATION OF NETWORK PROPERTIES

In this section, we first describe our working dataset and present the metrics used for characterization. Then, we present our experimental evaluation and discuss the results.

A. Description of the Dataset

To determine the properties and features of a swarm of nano-satellites, we use synthetic data generated in Matlab. Each nano-satellite orbits the Moon and follows Kepler's laws. The trajectory parameters of the first node are set manually on a given orbit. In contrast, the parameters of the others are generated dynamically with a random offset to the first one. The dataset has the following specificities:

• The data represent the coordinates of |N | = 100 nodes in the Moon-centered coordinate system, distributed in a sphere of 100 km in diameter in its initial state; • the (x,y,z) coordinates of the nodes are sampled every 10 seconds; • the simulation duration is T = 100, 000 seconds; • a nano-satellite performs a complete revolution in approximately 18,000 seconds.

Classic geolocalization is impossible to achieve because of the lack of GPS tracking in outer space. The Moon-centered coordinate system is convenient in this case, as the nanosatellites will process the coordinates to perform a distancebased peer localization. Relative distances and positions are indeed known because they are necessary for interferometry image computation.

B. Metrics of Interest

In this paper, we adopt a graph theory approach because graph theory metrics are well-suited for network characterization [START_REF] Hwang | Water distribution system classification using system characteristics and graph-theory metrics[END_REF]. We focus on characterizing three significant network properties: the node density, the network connectivity, and the ISL availability. Fig. 2: Evolution of the degree distribution with time through its average value and statistical range.

1) Node Density: the node density is generally characterized by the node degree, denoted for a node u as δ(u), which is the number of nodes within its connection range. By definition, the degree strongly depends on the chosen connection range. The knowledge of the neighborhood can be extended beyond the connection range by computing the k-vicinity [START_REF] Phe-Neau | The strength of vicinity annexation in opportunistic networking[END_REF]. The k-vicinity, denoted δ k (u), characterizes this extended neighborhood as the number of nodes reachable within at most k hops from the source node (note that δ 1 (u) = δ(u)). The analysis of the k-vicinity is of interest in our case because it shows that some nodes with few direct neighbors can reach a larger set by performing additional hops, making it a valuable metric for routing configuration. It can be useful in two cases:

Case 1: there is a limited number of message copies to transmit, and a choice has to be made for the next hop: a node with a larger k-vicinity would be preferred over an isolated node. This case is not relevant here because the messages are broadcasted and thus not limited in terms of copies.

Case 2: we need to optimize the energy consumption by minimizing the number of hops of a message: transmitting the message in priority to nodes with a high k-vicinity can decrease the total duration of the transmission. By doing so, it is possible to adapt the route of the messages and optimize the propagation. This is the case to be studied in the context of our mission.

2) Network Connectivity: the network connectivity can be characterized through various metrics, but we choose to focus only on the graph density and the connected components. The graph density (not to be confused with node density) gives an overview of the connectivity by evaluating the proportion of ISL in the graph, and is defined as follows:

ρ t = 2|E t | |N |(|N | -1)
The second part is the analysis of the connected components, which are defined as groups of nodes connected together. The higher the number of connected components, the weaker is the connectivity of the graph. We denote C t the set of connected components at a time t:

C t = {c i t |c i t ⊂ N t } and c 0 t ∪ c 1 t ... ∪ c n t = N t
3) ISL availability: we characterize the availability of the ISL in the swarm by analyzing the Inter-Contact Times (ICT) [START_REF] Conan | Characterizing pairwise inter-contact patterns in delay tolerant networks[END_REF]. The ICT is defined, for a pair of nodes u and v, as the proportion of time spent without a direct connection between them, and is denoted τ (u, v). As such, the analysis of the ICT distribution provides direct information on the availability and consistency of the ISL within the swarm. The ICT is usually computed on a given period of time (e.g. on a complete revolution around the Moon).

C. Experimental Evaluation and Results

We model the swarm system in Python3 by creating a simple yet adapted simulation tool called swarm sim1 . This simulator is a Python3 module that comprises Node and Swarm objects a well as many basic operation functions (such as trajectory simulation, metrics computation, etc.). We choose to implement our swarm simulator because the existing frameworks do not match our case study. No tool is efficient in both networking and spatial mechanics: indeed, in [START_REF] Raveneau | A multi-level FREAK DTN: Taking care of disconnected nodes in the IoT[END_REF], [START_REF] Whitbeck | HYMAD: Hybrid DTN-MANET routing for dense and highly dynamic wireless networks[END_REF], the authors use the DTN simulator The ONE presented in [START_REF] Keränen | The ONE simulator for DTN protocol evaluation[END_REF], which performs well in network-layer simulation from 2D data, but cannot simulate 3D data in the Moon-centered coordinate system. The Hypatia simulator is presented in [START_REF] Kassing | Exploring the" Internet from space" with Hypatia[END_REF] and displays many functionalities for constellation networks. However, in this model the ISL are static and not dynamic, making this simulator not adapted to swarms of nano-satellites. The global parameters of the simulation are listed in Tab. I. At each timestamp, we analyze the network connectivity by computing the number of connected components and the graph density. The evolution of these metrics is drawn in Fig. 4. We also compute the ICT for all pairs of nodes in N t : the ICT distribution is drawn in Fig. 5. In parallel, for each node, we compute the degree and k-vicinity to grasp the dynamics of the node density. Fig. 2 shows the evolution of the distribution of node degrees through its average value and statistical range. Fig. 3 shows the evolution of the average k-vicinity with respect to the node degree.

1) Results on the Node Density: because the average node degree is strongly linked to the inter-satellite distances, the evolution of the instantaneous average node degree shown on Fig. 2 highlights a periodicity in the mobility within the swarm and thus in the node density. The minimum node degree is mostly equal to 0, even when the average degree is high: it shows that the statistical range of node degrees within the swarm can be very wide, with highly and weakly connected nodes. On the other hand, the evolution of the average kvicinity with respect to the node degree as drawn in Fig. 3 was predictable: intuitively, we assumed that the extended number of reachable nodes would increase with additional hops, which is indeed the case in this graph. By performing only one additional hop (k = 2), it is possible to double the number of node neighbors, no matter the degree of the node.

2) Results on the Network Connectivity: Fig. 4 shows the evolution, on a complete revolution around the Moon, of the number of connected components in the graph as well as the graph density. The cross-evolution of these metrics is interesting because they are inverted: a high number of connected components implies fewer edges in the graph, hence a low graph density. On the contrary, a very high graph density implies that the nodes are connected within a single, bigger component. We can also notice the graph density follows the same evolution as the average node degree. Nonetheless, the highest graph density value merely equals 0.11, i.e., we can observe up to 11% of ISL out of all possibilities in the strongest connected state. 3) Results on the ISL availability: Fig. 5 shows the ICT distribution computed for each theoretical ISL on a complete revolution around the Moon. We call theoretical ISL all the possible pairs of nodes, while existing ISL refers to the edges actually observed in E t . In this figure, the average ICT on existing pairs is above 80%, and the average on theoretical ISL is higher. This means that an ISL is available less than 20% of the time on average. But the most interesting information are the extreme values. First, the large number of unavailable links, represented by the pairs at 100% ICT, shows that some pairs of nodes are too far away from each other and will never establish a direct contact (these pairs are not considered existing ISL). We can also notice that a significant number of pairs are at 0% ICT, which means they are always connected and thus act as a backbone in the network.

Although the analysis of these five metrics provides various information on the swarm network, the central characteristic is its important heterogeneity on all studied aspects.

D. Routing Perspectives

The characterization of the nano-satellite swarm highlights the peculiar nature of the network. Indeed, the overall connectivity of the network is weak (based on the average number of neighbors and the ISL availability) and evolves periodically in orbit around the Moon. These results encourage us to prefer a DTN routing scheme over a MANET one. However, an analysis based on average values and tendencies alone can hide the specific dynamics of the network. In our case, the network is actually very heterogeneous in terms of density, connectivity, and even ISL availability. Therefore, the routing scheme will perform better in some cases than in others. Thus, the best solution is to develop a hybrid approach to alleviate the heterogeneous nature of the network, e.g. with a MANET routing within the high-density zones and DTN routing within the low-density zones, as proposed in [START_REF] Whitbeck | HYMAD: Hybrid DTN-MANET routing for dense and highly dynamic wireless networks[END_REF] with the HYMAD protocol.

V. CONCLUSION

This paper has highlighted the importance of studying nanosatellite swarm systems more thoroughly for the scientific community. We have presented an efficient, methodology based on graph theory for the characterization of any swarm network operating in space, intending to grasp the dynamics of the network and derive its properties to determine the most adapted routing scheme. Indeed, a wireless mobile network may require either MANET or DTN routing schemes, but a complex system such as a swarm of nano-satellites orbiting the Moon needs a hybrid protocol to adapt to any situation in space.

In our case, we showed that the node density, network connectivity and ISL availability within the swarm oscillate periodically according to the position of the swarm around the Moon and are very heterogeneous. Hence, the routing strategy must be adaptive to be efficient at all times and for all nodes. Because we have high-and low-density zones within the swarm, a hybrid MANET/DTN protocol such as HYMAD seems to be the best-suited routing scheme based on our simulation results.

In our future work, we aim to study the energy consumption related to the data transmission within the swarm, and find the best way to minimize it. Our first objective will be to limit the number of transmissions of data packets from one node to the rest of the swarm by performing a fair graph division. Then, we will study the impact of different routing strategies on the energy consumption of the swarm, which will be a key criterion for the choice of the routing strategy.

  (a) Trajectory of the nodes on their orbit around the Moon. (b) Topology of G at snapshot (1). (c) Topology of G at snapshot (2).

Fig. 1 :

 1 Fig. 1: Dynamic topology of graph G. The orbit of G is defined as the average orbit of its |N | nodes. Each node position is derived with a random offset to the first node position, which is fixed (see Section IV-A).

Fig. 3 :

 3 Fig. 3: Average k-vicinity with respect to node degree, for three different values of k.

Fig. 4 :

 4 Fig. 4: Cross-evolution of the number of connected components and the graph density over a complete revolution around the Moon.

Fig. 5 :

 5 Fig. 5: ICT distribution for all pairs of nodes. Pairs on 100% ICT are considered theoretical only. The average value on existing ISL is over 80%.

TABLE I :

 I Simulation parameters.

Link to the source code available upon request.

ACKNOWLEDGMENTS

The authors would like to thank the CNES as well as TeSA laboratory for supporting their work.


Pierre Degond
email: pierre.degond@math.univ-toulouse.fr

Min Tang

All speed scheme for the low mach number limit of the Isentropic Euler equation

Keywords: AMS subject classification: 65M06,65Z05,76N99,76L05 Low Mach number, isentropic euler equation, compressible flow, incompressible limit, asymptotic preserving, Lax-Friedrich scheme

An all speed scheme for the Isentropic Euler equation is presented in this paper. When the Mach number tends to zero, the compressible Euler equation converges to its incompressible counterpart, in which the density becomes a constant. Increasing approximation errors and severe stability constraints are the main difficulty in the low Mach regime. The key idea of our all speed scheme is the special semi-implicit time discretization, in which the low Mach number stiff term is divided into two parts, one being treated explicitly and the other one implicitly. Moreover, the flux of the density equation is also treated implicitly and an elliptic type equation is derived to obtain the density. In this way, the correct limit can be captured without requesting the mesh size and time step to be smaller than the Mach number. Compared with previous semi-implicit methods [11, 13, 27], nonphysical oscillations can be suppressed. We develop this semi-implicit time discretization in the framework of a first order local Lax-Friedrich (LLF) scheme and numerical tests are displayed to demonstrate its performances.

Introduction

Singular limit problems in fluid mechanics have drawn great attentions in the past years, like low-Mach number flows, magneto-hydrodynamics at small Mach and Alfven numbers and multiple-scale atmospheric flows. As mentioned in [START_REF] Klein | Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics I: one-dimensional flow[END_REF], the singular limit regime induces severe stiffness and stability problems for standard computational techniques. In this paper we focus on the simplest Isentropic Euler equation and propose a numerical scheme that is uniformly applicable and efficient for all ranges of Mach numbers.

The problem under study is the Isentropic Euler equation

∂ t ρ ε + ∇ • (ρ ε u ε ) = 0, ∂ t (ρ ε u ε ) + ∇ ρ ε u ε ⊗ u ε + 1 ε 2 ∇p ε = 0. ( 1
)
where ρ ε , ρ ε u ε is the density and momentum of the fluid respectively and ε is the scaled Mach number. This is one of the most studied nonlinear hyperbolic systems. For standard applications, the equation of state takes the form

p(ρ) = Λρ γ , (2)
where Λ, γ are constants depending on the physical problem.

It is rigorously proved by Klainerman and Majda [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] that when ε → 0, i.e. when the fluid velocity is small compared with the speed of sound [START_REF] Constantin | On the Euler equations of incompressible fluids[END_REF], the solution of (1) converges to its incompressible counterpart. Formally, this can be obtained by inserting the expansion

ρ ε = ρ 0 + ε 2 ρ (2) + • • • , u ε = u 0 + ε 2 u (2) + • • • , (3)
into [START_REF] Bijl | A unified method for computing incompressible and compressible flows in boundary-fitted coordinates[END_REF] and equate the same order of ε. The limit reads as follows [START_REF] Klainerman | Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids[END_REF][START_REF] Klein | Asymptotic adaptive methods for multi-scale problems in fluid mechanics[END_REF]:

ρ = ρ 0 , (4a) ∇ • u 0 = 0, ( 4b
) ∂ t u 0 + ∇(u 0 ⊗ u 0 ) + ∇p 0 = 0. ( 4c
)
Here p 0 is a scalar pressure that can be viewed as the Lagrange multiplier which enforces the incompressibility constraint. Physically, this limit means that in slow flows (compared with speed of sound), the factor 1/ε 2 in the momentum equation in front of the pressure gradient generates fast pressure waves, which makes the pressure and therefore, the density, uniform in the domain [START_REF] Park | Multiple pressure variables methods for fluid flow at all Mach numbers[END_REF][START_REF] Munz | The extension of incompressible folw solvers to the weakly compressible regime[END_REF].

For atmosphere-ocean computing or fluid flows in engineering devices, when ε is small in [START_REF] Bijl | A unified method for computing incompressible and compressible flows in boundary-fitted coordinates[END_REF], standard numerical methods become unacceptably expensive. Indeed, [START_REF] Bijl | A unified method for computing incompressible and compressible flows in boundary-fitted coordinates[END_REF] has wave speeds of the form

λ = u ε ± 1 ε p ′ (ρ ε ),
where p ′ (ρ ε ) is the derivative with respect to ρ ε . If a standard hyperbolic solver is used, the CFL requirement is ∆t = O(ε∆x). Moreover in order to maintain stability, the numerical dissipation required by the hyperbolic solver is proportional to

|λ |. If |λ | = O( 1 ε
), in order to control the diffusion, we need to have ∆x = o(ε r ), where r is some appropriate constant. Thus the stability and accuracy highly depend on ε.

Our aim is to design a method whose stability and accuracy is independent of ε. The idea is to find an asymptotic preserving (AP) method, i.e. a method which gives a consistent discretization of the isentropic Euler equations (1) when ∆x, ∆t resolve ε, and a consistent discretization of the incompressible limit (4) when ε → 0 (∆x, ∆t being fixed). The efficiency of AP schemes at the low Mach number regime can be proved similarly as in [START_REF] Golse | The Convergence of Numerical Transfer Schemes in Diffusive Regimes I: The Discrete-Ordinate Method[END_REF]. The key idea of our all speed scheme is a specific semi-implicit time discretization, in which the low Mach number stiff term is divided into two parts, one part being treated explicitly and the other one implicitly. Moreover, the flux of the density equation is also treated implicitly. For the space discretization, when ε is O [START_REF] Bijl | A unified method for computing incompressible and compressible flows in boundary-fitted coordinates[END_REF], even if the initial condition is smooth, shocks will form due to the nonlinearity of the div ρ ε u ε ⊗ u ε term and shock capturing methods should be employed here.

In the literature, lots of efforts have been made to find numerical schemes for the compressible equation that can also capture the zero Mach number limit [START_REF] Bijl | A unified method for computing incompressible and compressible flows in boundary-fitted coordinates[END_REF][START_REF] Degond | Mach-number uniform asymptoticpreserving gauge schemes for compressible flows[END_REF][START_REF] Park | Multiple pressure variables methods for fluid flow at all Mach numbers[END_REF][START_REF] Munz | The extension of incompressible folw solvers to the weakly compressible regime[END_REF][START_REF] Haack | All speed asymptotic preserving schemes for compressible flows[END_REF]. In [START_REF] Bijl | A unified method for computing incompressible and compressible flows in boundary-fitted coordinates[END_REF], Bijl and Wesseling split the pressure into thermodynamic and hydrodynamic pressure terms and solve them separately. Similar to this approach, the multiple pressure variable (MPV) method was proposed by Munz et al. in [START_REF] Park | Multiple pressure variables methods for fluid flow at all Mach numbers[END_REF][START_REF] Munz | The extension of incompressible folw solvers to the weakly compressible regime[END_REF]. There is also some recent work by J. Hauck, J-G. Liu and S. Jin [START_REF] Haack | All speed asymptotic preserving schemes for compressible flows[END_REF]. Their approach involves specific splitting of the pressure term. We avoid using this splitting, the proper design of which seems very crucial in some cases.

Some similar ideas can be found in the ICE method, which is designed to adapt incompressible flow computation techniques using staggered meshes to the simulation compressible flows. The method was first introduced by Harlow and Amsdan in 1965 and 1971 [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF][START_REF] Harlow | A numerical fluid dynamics calculation method for all flow speeds[END_REF] and is called Implicit Continuous-fluid Eulerian (ICE) technique. It is used to simulate single phase fluid dynamic problems with all flow speeds. They introduce two parameters in the continuity equation and the momentum equation to combine information from both previous and forward time steps. However this method is not conservative, which leads to discrepancies in the shock speeds. Additionally it suffers from small wiggles when there are moving contact discontinuities. The first problem was solved by an iterative method, for example SIMPLE [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF], or PISO [START_REF] Issa | The computation of compressible and incompressible flow of fluid with a free surface[END_REF]. In some recent work, Heul and Wesseling also find a conservative pressure-correction method [START_REF] Van Der Heul | A conservative pressurecorrection method for flow at all speeds[END_REF]. All these methods are based on the so called MAC staggered mesh in order to be consistent with the staggered grid difference method for the incompressible Euler equations [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF]. Specifically, if we write the simplified ICE technique presented in [START_REF] Bonner | Compressible subsonic flow on a staggered grid[END_REF] in conservative form, we are led to the semi-discrete framework:

ρ * ε -ρ n ε ∆t + ∇ • (ρ ε u ε ) n = 0, (ρ ε u ε ) * -(ρ ε u ε ) n ∆t + ∇(ρ n ε u n ε ⊗ u n ε ) = 0, (5)
ρ n+1 ε -ρ * ε ∆t + ∇ • (ρ ε u ε ) n+1 -(ρ ε u ε ) * = 0, (ρ ε u ε ) n+1 -(ρ ε u ε ) * ∆t + 1 ε 2 ∇p(ρ n+1 ε ) = 0. (6)
By substituting the gradient of the second equation of ( 6) into its first equation and using the results of the first equation ( 5), ρ ε can be updated by solving an elliptic equation which does not degenerate when ε → 0. We use a similar idea in our method. However, we do not use the predictorcorrector procedure but we rather discretize the problem in a single step. We use standard shock capturing schemes which allows to guarantee the conservativity and the desired artificial viscosity. We only use implicit evaluations of the mass flux and pressure gradient terms to ensure stability and provide an extremely simple way to deal with the implicitness. Additionally, we propose a modification of the implicit treatment of the pressure equation. Indeed, using a similar idea as in [START_REF] Haack | All speed asymptotic preserving schemes for compressible flows[END_REF], we split the pressure into two parts and put α p(ρ ε ) into the hyperbolic system. This makes the first system no longer be weakly hyperbolic and much more stable. The numerical results show the advantage of our method in the following sense:

• The method is in conservative form and can capture the right shock speeds.

• The non-physical oscillations [START_REF] Haack | Oscillatory Behavior of Asymptotic-Preserving Splitting Methods for a Linear Model of Diffusive Relaxation[END_REF] can be suppressed by choosing the proper value of the parameter which determines the fraction of implicitness used in the evaluation of the pressure gradient term. The choice of this parameter depends on the time and space step and on the specific problem.

In this paper we only use the first order LLF scheme. Higher order space and time discretizations will be subject of future work. The main objective of this work is to show that the semi-discrete time discretization provides a framework for developing AP methods for singular limit problems. Similar ideas can be extended to the full Euler equation and more complicated fluid model and have also been used in other contexts such as quasineutrality limits [START_REF] Crispel | An asymptotic preserving scheme for the two-fluid Euler-Poisson model in the quasineutral limit[END_REF][START_REF] Degond | Analysis of an asymptotic preserving scheme for the Euler-Poison system in the quasineutral limit[END_REF] and magnetized fluids under stong magnetic fields [START_REF] Degond | An asymptotic preserving scheme for the Euler equations in a strong magnetic field[END_REF].

The organization of this paper is as follows. Section 2 exposes the semiimplicit scheme and its capability to capture the incompressible limit is proved. The detailed one dimensional and two dimensional fully discretized schemes and their AP property are presented in section 3 and 4 respectively. In section 5, how to choose the ad-hoc parameter is discussed and finally, some numerical tests are given in section 6 to discuss the stability and accuracy of our scheme. The efficiency at both the compressible and low mach number regime are displayed. Finally, we conclude in section 6 with some discussion.

Time Semi-discrete scheme

Let ∆t be the time step, t n = n∆t, n = 0, 1, • • • and let the 'n' superscript denote the approximations at t n . The semi-discrete scheme for the nth time step is

ρ n+1 ε -ρ n ε ∆t + ∇ • (ρ ε u ε ) n+1 = 0, (7)
(ρ ε u ε ) n+1 -(ρ ε u ε ) n ∆t + div ρ n ε u n ε ⊗ u n ε + α p(ρ n ε ) + 1 -αε 2 ε 2 ∇p(ρ n+1 ε ) = 0, ( 8
)
where α is an ad-hoc parameter which satisfies α ≤ 1/ε 2 . The choice of α depends on the space and time steps and on the fluid speed. When the shock is strong, α should be bigger, which means that the system should be more explicit to follow the discontinuity more closely. We discuss the choice of α for specific equations of state in this paper and test its effect numerically. It depends on the required accuracy, the small parameter ε and the shock amplitude in a sometimes quite complex way.

Rewriting the momentum equation ( 8) as

(ρ ε u ε ) n+1 = (ρ ε u ε ) n -∆t∇ ρ n ε u n ε ⊗ u n ε + α p(ρ n ε ) -∆t 1 -αε 2 ε 2 ∇P(ρ n+1 ε )
and substituting it into the density equation, one gets

ρ n+1 ε -∆t 2 1 -αε 2 ε 2 ∆P(ρ n+1 ε ) = φ (ρ n ε , u n ε ) (9)
which is an elliptic equation that can be solved relatively easily. Here

φ (ρ n ε , u n ε ) = ρ n ε -∆t∇ • (ρ n ε u n ε ) + ∆t 2 ∇ • ∇ ρ n ε u n ε ⊗ u n ε + α p(ρ n ε ) . ( 10
)
The Laplace operator in ( 9) can be approximated by

∇ P ′ (ρ n ε )∇ρ n+1 ε and (9) becomes ρ n+1 ε -∆t 2 1 -αε 2 ε 2 ∇ • P ′ (ρ n ε )∇ρ n+1 ε = φ (ρ n ε , u n ε ), (11)
Though shocks will form for the original system (7)(8), we always add some numerical diffusion terms so that φ (ρ n ε , u n ε ) is smooth. Then so is ρ n+1 ε . When we implement this method, ρ n+1 ε can be obtained from (9) first and u n ε is then updated by the momentum equation ( 8) afterwards. Therefore, apart from the resolution of the elliptic equation [START_REF] Haack | All speed asymptotic preserving schemes for compressible flows[END_REF], the scheme only involves explicit steps.

We now show that the scheme (7)( 8) is asymptotic preserving. We introduce the formal expansion

ρ n ε (x) = ρ n 0c + ερ n (1) (x) + ε 2 ρ n (2) (x) + • • • , u n ε = u n 0 (x) + εu n (1) (x) + • • • . (12)
In the sequel, the 'c' in the index means that the quantity is independent of space. When ∆x, ∆t are fixed and ε goes to 0 in (9), we formally have ∆P(ρ n+1 0 ) = 0, which implies that ρ n+1 0 is independent of space, where ρ n+1 0 is the limit of ρ n+1 ε when ε → 0. Thus we have

ρ n+1 0c -ρ n 0c ∆t + ∇ • (ρ 0c u 0 ) n+1 = 0 ( 13
)
by equating the O(1) terms in the density equation [START_REF] Degond | Analysis of an asymptotic preserving scheme for the Euler-Poison system in the quasineutral limit[END_REF]. Integrating (13) over the computational domain, one gets

|Ω| ρ n+1 0c -ρ n 0c ∆t = -ρ n+1 0c Ω ∇ • (u 0 ) n+1 = -ρ n+1 0c ∂ Ω n • u n+1 0 . (14)
As discussed in [START_REF] Haack | All speed asymptotic preserving schemes for compressible flows[END_REF], for wall boundary condition, periodic boundary condition and open boundary condition, [START_REF] Issa | The computation of compressible and incompressible flow of fluid with a free surface[END_REF] gives

ρ n+1 0c = ρ n 0c , (15)
that is ρ 0 is also independent of time. Thus (13) also implies

∇ • u n+1 0 = 0. ( 16
)
Then, by using the fact that the curl of the gradient of any scalar field is always zero, the curl of the O(1) terms of the momentum equation ( 8) becomes

∇ × u n+1 0 -u n 0 ∆t + ∇ × ∇ u n 0 ⊗ u n 0 = 0. ( 17
)
Thus

u n+1 0 -u n 0 ∆t + ∇ u n 0 ⊗ u n 0 + ∇p n (2) = 0, ( 18
)
where p n (2) is some scalar field. Equations ( 15), ( 16), [START_REF] Klein | Asymptotic adaptive methods for multi-scale problems in fluid mechanics[END_REF] are the semi-discretization in time of (4) and thus the scheme [START_REF] Degond | Analysis of an asymptotic preserving scheme for the Euler-Poison system in the quasineutral limit[END_REF], ( 8) is consistent with the low Mach number limit ε → 0 of the original compressible Euler equation. This statement is exactly saying that the scheme is AP. We can see that, in order to obtain the stability and AP properties, it is crucial to treat the flux in the density equation [START_REF] Degond | Analysis of an asymptotic preserving scheme for the Euler-Poison system in the quasineutral limit[END_REF] implicitly.

Letting U = (ρ ε , ρ ε u ε ) T , we can write [START_REF] Degond | Analysis of an asymptotic preserving scheme for the Euler-Poison system in the quasineutral limit[END_REF], ( 8) abstractly as

U n+1 -U n ∆t + ∇ • F(U n+ 1 2 ) + QU n+1 = 0, (19)
where

F(U n+1/2 ) = (ρ ε u ε ) n+1 ρ n ε u n ε ⊗ u n ε + α p(ρ n ε ) , Q = 0 0 1-αε 2 ε 2 ∇P 0 . ( 20
)
Here P is an operator on ρ ε and U n+1/2 reminds that the flux is partly implicit and partly explicit. This semi-discretization gives us a framework for developing AP schemes that can capture the incompressible limit. Now we are left with the problem of how discretizing the space variable. Because shocks can form, considerable literature has been devoted to the design of high resolution methods that can capture the correct shock speed. Upwind schemes and central schemes are among the most widely used Godunov type schemes [START_REF] Kurganov | New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations[END_REF][START_REF] Kurganov | Solution of two-dimensional Riemann problems for gas dynamics without Riemann problem solvers[END_REF][START_REF] Leveque | Numerical methods for conservation laws[END_REF].

In the present paper, the hyperbolic operator

U n+1 -U n ∆t + ∇ • F(U n+ 1 2 )
is approximated by an upwind hyperbolic solver and the stiff 1/ε 2 factor in front of the pressure term is treated implicitly. The implicitness of the density flux is treated by combining it with the momentum equation. For simplicity, in the present work we only consider the first order modified Lax-Friedrich scheme with local evaluation of the wave-speed in the current and neighboring cell.

Full time and space discretization: One dimensional case

For simplicity, we consider the domain Ω = [0, 1]. Using a uniform spatial mesh with ∆x = 1/M, M being an positive integer, the grid points are defined as

x j := j∆x, j = 0, 1, • • • , M.
The flux and Jacobian matrix of ( 19) become

F(U ) = ρ ε u ε ρ ε u 2 ε + α p(ρ ε ) , F′ (U ) = 0 1 -u 2 ε + α p ′ (ρ ε ) 2u ε , (21)
so, the wave speeds are

λ = u ε ± α p ′ (ρ ε ). (22)
Let U j be the approximation of U (x j ) and let

A j+ 1 2 (t) = max |λ j |, |λ j+1 | . ( 23
)
These are the local maximal wave-speeds in the current and neighboring cells. We discretize [START_REF] Kurganov | New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations[END_REF] in space as follows:

U n+1 j -U n j ∆t + F j+ 1 2 (U n+ 1 2 ) -F j-1 2 (U n+ 1 2 ) ∆x + Q j U n+1 = 0, (24)
where

F j±1/2 (U n+ 1 2 ) is the numerical flux F j+ 1 2 (U n+ 1 2 ) = 1 2 F + j+ 1 2 (U n+ 1 2 ) + F - j+ 1 2 (U n+ 1 2 ) (25)
and

F + j+ 1 2 (U n+ 1 2 ) = (ρ ε u ε ) n+1 j + A n j+ 1 2 ρ n ε j (ρ ε u ε ⊗ u ε ) n j + α p(ρ n ε j ) + A n j+ 1 2 (ρ ε u ε ) n j , F - j+ 1 2 (U n+ 1 2 ) = (ρ ε u ε ) n+1 j+1 -A n j+ 1 2 ρ n ε j+1 (ρ ε u ε ⊗ u ε ) n j+1 + α p(ρ n ε j+1 ) -A n j+ 1 2 (ρ ε u ε ) n j+1
and

QU n+1 j = 0 1-αε 2 ε 2 1 2∆x P(ρ n+1 ε j+1 ) -P(ρ n+1 ε j-1 ) . Let q = ρu, (26)
and F (1) , F (2) denote the first and second element of F respectively, we can rewrite the momentum discretization in [START_REF] Park | Multiple pressure variables methods for fluid flow at all Mach numbers[END_REF] as follows:

q n+1 ε j = q n ε j -∆tD x j F (2) ρ n ε , u n ε - 1 -αε 2 ε 2 ∆t 2∆x p(ρ n+1 ε j+1 ) -p(ρ n+1 ε j-1 ) . ( 27
)
Here

D x j u = u j+1/2 -u j-1/2 ∆x .
By substituting [START_REF] Van Der Heul | A conservative pressurecorrection method for flow at all speeds[END_REF] into the density equation in [START_REF] Park | Multiple pressure variables methods for fluid flow at all Mach numbers[END_REF], one gets

ρ n+1 ε j - (1 -αε 2 )∆t 2 4ε 2 ∆x 2 p(ρ n+1 ε j+2 ) -2p(ρ n+1 ε j ) + p(ρ n+1 ε j-2 ) = Dφ (ρ n ε , q n ε ), (28)
where [START_REF] Haack | Oscillatory Behavior of Asymptotic-Preserving Splitting Methods for a Linear Model of Diffusive Relaxation[END_REF]. We notice that (28) is a discretization of the elliptic equation [START_REF] Haack | All speed asymptotic preserving schemes for compressible flows[END_REF]. We can update q n+1 ε through [START_REF] Van Der Heul | A conservative pressurecorrection method for flow at all speeds[END_REF] afterwards. To obtain ρ n+1 ε in (28), a nonlinear system of equations needs to be solved. One possible way to simplify it is to replace ∇P(ρ n+1 ε ) by P ′ (ρ n ε )∇ρ n+1 ε , so that the following linear system is obtained:

Dφ (ρ n ε , q n ε ) = ρ n ε -∆tD x j F (1) (ρ n ε , u n ε ) + ∆t 2 2∆x D x j+1 -D x j-1 F (2) (ρ n ε , u n ε ) (29) is a discretization of φ (ρ n ε , u n ε ) in
ρ n+1 ε j - (1 -αε 2 )∆t 2 4ε 2 ∆x 2 p ′ (ρ n ε j+1 ) ρ n+1 ε j+2 -ρ n+1 ε j -p ′ (ρ n ε j-1 ) ρ n+1 ε j -ρ n+1 ε j-2 = Dφ (ρ n ε , q n ε ). ( 30
)
This is a five point scheme which is too much diffusive, especially near the shock. One possible improvement is that instead of (30), we use the following three points discretization

ρ n+1 ε j - (1 -αε 2 )∆t 2 ε 2 ∆x 2 p ′ (ρ n ε j+1 ) ρ n+1 ε j+1 -ρ n+1 ε j -p ′ (ρ n ε j ) ρ n+1 ε j -ρ n+1 ε j-1 = Dφ (ρ n ε , q n ε ). (31)
After obtaining ρ n+1 ε , we can substitute it into [START_REF] Van Der Heul | A conservative pressurecorrection method for flow at all speeds[END_REF] to get q n+1 ε j . To summarize, three schemes are proposed here: (28), ( 27); (30), ( 27) and (31), [START_REF] Van Der Heul | A conservative pressurecorrection method for flow at all speeds[END_REF]. To investigate the AP property, we take (30), ( 27) as an example. The proofs for the other two schemes are similar. By substituting the following expansion

ρ n ε j = ρ n 0c + ε 2 ρ n (2) j + • • • , q n ε j = q n 0c + εq n (2) j + • • • , (32)
into (30), the O( 1 ε 2 ) terms give that ρ n+1 (0) j = ρ n+1 0c is constant in space by using the periodic boundary condition, and thus:

ρ n+1 ε j = ρ n+1 (0)c + ε 2 ρ n+1 (2) j + • • • .
Summing (30) over all the grid points, one gets

ρ n+1 0c = ρ n 0c = ρ 0c , (33a)
which implies that ρ 0 is independent of time and space. Thus, the O(1) terms of (30) are

p ′ (ρ n 0c ) ρ n+1 (2) j+2 -ρ n+1 (2) j -p ′ (ρ n 0c ) ρ n+1 (2) j -ρ n+1 (2) j-1 = 0,
by recalling that the O(1) terms of both ρ n ε and q n ε are constant in space. Then the periodic boundary condition gives

ρ n+1 (2) j = ρ n+1 (2)c , (33b)
which gives that ρ n+1 (2) is also independent of space. Therefore from ( 2), [START_REF] Van Der Heul | A conservative pressurecorrection method for flow at all speeds[END_REF],

q n+1 0 j = q n 0 j = q n 0c . (33c)
In one dimension, (33) is the discretization of ( 15), ( 16), [START_REF] Klein | Asymptotic adaptive methods for multi-scale problems in fluid mechanics[END_REF] when periodic boundary conditions apply and thus is consistent with the incompressible limit. In fact all the three methods proposed here are AP.

Full time and space discretization: Two dimensional case

We consider the domain

Ω = [0, 1] × [0, 1].
For M 1 , M 2 two positive integers, we use a uniform spatial mesh ∆x = 1/M 1 , ∆y = 1/M 2 . The grid points are

(x i , y j ) := (i∆x, j∆y), i = 0, • • • , M 1 ; j = 0, • • • , M 2 Now U = (ρ ε , q (1)
ε , q

ε ) T and U i, j is the numerical approximation of U (x i , y j ). Let

G 1 (U ) =   ρ ε u ε1 ρ ε u 2 1 + α p(ρ ε ) ρ ε u 1 u 2   , G 2 (U ) =   ρ ε u ε2 ρ ε u 1 u 2 ρ ε u 2 2 + α p(ρ ε )   . ( 34
)
and

Q = 1 -αε 2 ε 2   0 0 0 ∂ x P 0 0 ∂ y P 0 0   .
Eq. ( 19) can be written as

∂ t U + ∂ x G 1 (U ) + ∂ y G 2 (U ) + QU = 0. Denote G 1 (U n+ 1 2 ) =   (ρ ε u ε1 ) n+1 ρ n ε (u n ε1 ) 2 + α p(ρ n ε ) ρ n ε u n ε1 u n 2   , G 2 (U n+ 1 2 ) =   (ρ ε u ε2 ) n+1 ρ n ε u n ε1 u n ε2 ρ n ε (u n ε2 ) 2 + α p(ρ n ε )   , Q =    0 0 0 1-αε 2 ε 2 D x P 0 0 1-αε 2 ε 2 D y P 0 0    , D x i j u = u i j+1 -u i j-1 2∆x , D y i j u = u i+1 j -u i-1 j 2∆y .
Now the eigenvalues of the two one-dimensional hyperbolic equations are

λ (1) = u 1 , u 1 ± α p ′ (ρ ε ), λ (2) = u 2 , u 2 ± α p ′ (ρ ε ).
The fully discrete scheme for the two dimensional problem is

U n+1 i j -U n i j ∆t + D x i j G 1 (U n+1/2 ) + 1 2 A i-1 2 , j D x i j--A i+ 1 2 , j D x i j+ U n +D y i j G 2 (U n+1/2 ) + 1 2 A i, j-1 2 D y i j--A i, j+ 1 2 D y i j+ U n + QU n+1 i j = 0, (35)
where

D x i j-u = u i j -u i-1 j ∆x , (D x i j+ u) = u i+1 j -u i j ∆x , D y i j-u = u i j -u i j-1 ∆y , D y i j+ u = u i j+1 -u i j ∆y ,
and

A n i+ 1 2 , j = max{|λ (1) i j |, |λ (1)
i+1, j |, |λ (2)
i j |, |λ (2)
i+1, j |}, A n i, j+ 1 2 = max{|λ (1) i j |, |λ (1) i, j+1 |, |λ (2) i j |, |λ (2)
i, j+1 |}.

(36) Let q ε be like in [START_REF] Van Der Heul | A conservative pressurecorrection method for flow at all speeds[END_REF]. Like in one dimension, we can substitute the expressions of q n+1 1i j , q n+1 2i j into the density equation and get the following discretized elliptic equation,

ρ n+1 εi j - (1 -αε 2 )∆t 2 4ε 2 1 ∆x 2 P(ρ n+1 εi+2, j ) -2P(ρ n+1 εi, j ) + P(ρ n+1 εi-2, j ) + (37) 1 ∆y 2 P(ρ n+1 εi, j+2 ) -2P(ρ n+1 εi, j ) + P(ρ n+1 εi, j-2 ) = Dφ i j (ρ n ε , q n 1ε , q n 2ε ),
where

Dφ i j (ρ n ε , q n 1ε , q n 2ε ) = ρ n ε -∆t D x i j q n ε1 + D y i j q n ε2 + 1 2 A i-1 2 D x i j--A i+ 1 2 , j D x i j+ + A i, j-1 2 D y i j--A i, j+ 1 2 D y i j+ ρ n ε +∆t 2 D x i j D x i j ρ n ε (u n ε1 ) 2 + α p(ρ n ε ) + D y i j D y i j ρ n ε (u n ε2 ) 2 + α p(ρ n ε ) +(D x i j D y i j + D y i j D x i j )ρ n ε u n ε1 u n ε2 + 1 2 D x i j (A i-1 2 , j D x i j--A i+ 1 2 , j D x i j+ )q n ε1 + 1 2 D y i j (A i-1 2 , j D x i j--A i+ 1 2 , j D x i j+ )q n ε2 + 1 2 D x i j (A i, j-1 2 D y i j--A i, j+ 1 2 D y i j+ )q n ε1 + 1 2 D y i j (A i, j-1 2 D y i j--A i, j+ 1 2 D y i j+ )q n ε2 . (38)
After obtaining ρ n+1 i j by (37), q n+1 1i j , q n+1 2i j can be updated by the momentum equation afterwards.

Similar to the one-dimensional case, the modified diffusion operator using a reduced stencil is as follows:

ρ n+1 εi j -∆t 2 1 -αε 2 ε 2 × × 1 ∆x 2 P ′ (ρ n εi, j+1 ) ρ n+1 εi, j+1 -ρ n+1 εi, j -P ′ (ρ n εi, j ) ρ n+1 εi, j -ρ n+1 εi, j-1 + 1 ∆y 2 P ′ (ρ n εi+1, j ) ρ n+1 εi+1, j -ρ n+1 εi, j -P ′ (ρ n εi, j ) ρ n+1 εi, j -ρ n+1 εi-1, j = φ (ρ n ε , q n 1ε , q n 2ε ). ( 39
)
Now we prove the AP property of our fully discrete scheme. Here only wellprepared initial conditions are considered, which means that there will be no shock forming in the solution. Then α can be chosen to be 0 to minimize the introduced numerical viscosity. Assuming that the expansions of ρ ε , u ε in (12) hold at time t n , when ε → 0, the O( 1ε 2 ) terms of (39) give 1

∆x 2 P ′ (ρ n 0i, j+1 ) ρ n+1 0i, j+1 -ρ n+1 0i, j -P ′ (ρ n 0i, j ) ρ n+1 0i, j -ρ n+1 0i, j-1 + 1 ∆y 2 P ′ (ρ n 0i+1, j ) ρ n+1 0i+1, j -ρ n+1 0i, j -P ′ (ρ n 0i, j ) ρ n+1 0i, j -ρ n+1 0i-1, j = 0.
When using periodic boundary conditions, one gets ρ n+1 0i j = ρ n+1 0c from [START_REF] Bonner | Compressible subsonic flow on a staggered grid[END_REF]. The time independence of ρ n+1 0 , similar to the one dimensional case, can be obtained by summing (39) over all the grid points. Accordingly we have

ρ n+1 εi j = ρ n 0c + ε 2 ρ n+1 (2)i j + • • • . ( 40
)
To prove the limiting behavior of u ε1 , u ε2 , we do not want to use the density equation because the diffusion operator with reduced stencil does not allow us to find the corresponding density equation. Therefore, we consider the O(1) term of (39),

ρ n+1 0i j -∆t 2 × × 1 ∆x 2 P ′ (ρ n 0c ) ρ n+1 (2)i, j+1 -ρ n+1 (2)i, j -P ′ (ρ n 0c ) ρ n+1 (2)i, j -ρ n+1 (2)i, j-1 + 1 ∆y 2 P ′ (ρ n 0c ) ρ n+1 (2)i+1, j -ρ n+1 (2)i, j -P ′ (ρ n 0c ) ρ n+1 (2)i, j -ρ n+1 (2)i-1, j = φ (ρ n 0 , q n 10 , q n 20 ). (41)
Moreover, noting the fact that

D x i j P(ρ n+1 ε ) = D x i j P(ρ n 0c + ε 2 ρ n+1 (2) + o(ε 2 )) = D x i j P(ρ n 0c ) + ε 2 D x i j ρ n+1 (2) P ′ (ρ n 0c ) + o(ε 2 ), = ε 2 D x i j ρ n+1 (2) P ′ (ρ n 0c ) + o(ε 2 ),
and similarly,

D y i j P(ρ n+1 ε ) = ε 2 D y i j ρ n+1 (2) P ′ (ρ n 0c ) + o(ε 2 ),
the O(1) terms of the momentum equations of (35) become

q n+1 01i j -q n 01i j ∆t + D x i j q 2 01 ρ 0 n + D y i j q 01 q 02 ρ 0 n + 1 2 A i-1 2 , j D x i j--A i+ 1 2 , j D x i j+ +A i, j-1 2 D y i j--A i, j+ 1 2 D y i j+ q 01 = D x i j P ′ (ρ n+1 0c )ρ n+1 (2) , (42a) q n+1 02i j -q n 02i j ∆t + D x i j q 01 q 02 ρ 0 n + D y i j q 2 02 ρ 0 n + 1 2 A i-1 2 , j D x i j--A i+ 1 2 , j D x i j+ +A i, j-1 2 D y i j--A i, j+ 1 2 D y i j+ q 02 = D y i j P ′ (ρ n+1 0c )ρ n+1 (2) . (42b)
Comparing (41) with ∆t * (D x i j (42a) + D y i j (42b)) , one gets

D x i j q n+1 01 + D y i j q n+1 02 = O(∆x∆t), ( 43
)
which is an approximation of [START_REF] Klainerman | Compressible and incompressible fluids[END_REF]. Moreover, it is obvious that (42) is a discretization of [START_REF] Klein | Asymptotic adaptive methods for multi-scale problems in fluid mechanics[END_REF]. Thus we obtain a full discretization of (4) in the limit ε → 0. Therefore, the two-dimensional scheme is also AP.

The ad-hoc parameter

In this section we illustrate how to choose ∆t and the parameter α by considering the simple state equation P(ρ ε ) = ρ ε . In this context, the fully discrete scheme (24) can be written as

ρ n+1 ε -ρ n ε ∆t + ∇ • q n+1 ε -∇ • q n ε + ∇ • q n ε = 0, q n+1 ε -q n ε ∆t + ∇ ρ n ε u n ε ⊗ u n ε + αρ n ε + 1-αε 2 ε 2 ∇ρ n+1 ε = 0. , ( 44
)
where ∇ is the centered difference while ∇ stands for the difference of fluxes. The latter is defined as follows (in one space-dimension for simplicity):

( ∇ • F(U n )) j = F j+ 1 2 (U n ) -F j-1 2 (U n ) ∆x ,
where the flux F j±1/2 is defined as in [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF]. By substituting ∇ • q n+1 ε q n from the momentum equation of (44) into its density equation, one gets

ρ n+1 ε -ρ n ε ∆t + ∇ • q n ε -∆t 1 -αε 2 ε 2 ∆ρ n+1 ε -∆t∇ • ∇ ρ n ε u n ε ⊗ u n ε + αρ n ε = 0. ( 45
)
The O(∆t) terms behave like a diffusion term which suppresses the oscillations at the discontinuity. Assuming that we use a first order explicit LLF scheme, the diffusions needed to damp out the oscillations in the mass and momentum equations are respectively:

1 2 (|u n ε | + 1 ε )∆x∆ρ n ε , 1 2 (|u n ε | + 1 ε )∆x∆q n ε . (46)
Here in (45), besides the O(∆t) terms, ∇ • q n ε also includes some numerical dissipation. By noting

ρ n+1 ε = ρ n ε -∆t∇ • q n+1 ε + O(∆t∆x),
the diffusion for ρ ε now is

1 2 |u ε | + √ α ∆x + ∆t ε 2 ∆ρ n ε + ∆t∆(ρ n ε u n ε ⊗ u n ε ) (47)
plus some higher order terms. Moreover, the diffusion for q ε is

1 2 |u ε | + √ α ∆x∆q n ε + ∆t 1 -αε 2 ε 2 ∆q n ε (48)
and some higher order terms. Comparing ( 46) and ( 47), (48), in order to suppress the oscillations at discontinuities we only need to have

1 2 |u n ε | + √ α ∆x + 1 -αε 2 ε 2 ∆t ≥ 1 2 |u n ε | + 1 ε ∆x, that is ∆t ≥ 1 2 ε∆x 1 + √ αε . ( 49
)
Moreover the CFL condition for the explicit part is

∆t ≤ σ ∆x max{|u ε | + √ α} , ( 50
)
where σ is the Courant number which is less than 1. We usually choose σ to be 0.5. Then the parameter α should satisfy

∆x 2∆t - 1 ε ≤ √ α ≤ σ ∆x ∆t -max{|u n | ε } (51)
according to (49), (50). Then the following constraint on ∆t should hold if we want the scheme to be stable and non-oscillatory

max{|u n ε |} + ∆x 2∆t ≤ σ ∆x ∆t + 1 ε . ( 52
)
The reason for the occurence of nonphysical oscillations when α = 0 lies in the fact that the diffusion is not large enough. In this case, with a simple reduction of ∆t, it is likely that the diffusion can no longer suppress the oscillations. This is why we need to introduce α to control the oscillations. But, from the analysis, no matter the value of α, as long as it is less than 1/ε 2 , the diffusion can never be sufficient when ∆t ≤ 1 4 ε∆x. In summary there is no specific way of choosing α < 1/ε that can guarantee that the nonphysical oscillations will disappear in any case. For well-prepared initial conditions in the low Mach number regime, because there is no shock formation in the solution, it is better to choose α as small as possible to get better accuracy, but if strong shocks exist in the solution, α should be big enough to suppress the oscillations. This is why the choice of α depends on the considered problem.

Numerical results

Three numerical examples will allow us to test the performances of the proposed schemes. In fact, three schemes are proposed in section 3 and 4, for example in one dimension: the scheme (24) without linearizing ∇P(ρ n+1 ε ) is denoted by "NL". We need to use Newton iterations to solve the nonlinear system. When ∇P(ρ n+1 ε ) is approximated by P ′ (ρ n ε )∇P(ρ n+1 ε ), the unknowns become a linear system. This scheme is represented by "L". "LD" denotes the scheme with the narrower stencil (31). Here we use well-prepared initial conditions of the form (3) and α = 1 for all the test cases.

In one dimension, let the computational domain be [a, b] and the mesh size be ∆x. The grid points are x j = a + ( j -1)∆x.

In the following tables, the L 2 norm of the relative error between the reference solutions u and the numerical ones

U e(U ) = U -u L 2 u L 2 = 1 M ∑ j |U j -u(x j )| 2 1 2 1 M e ∑ i |u(x i )| 2 1 2
are displayed.

Example 1 P(ρ ε ) and the initial conditions are chosen as

P(ρ ε ) = ρ 2 ε , ρ ε (x, 0) = 1, p ε (x, 0) = 1 -ε 2 /2 x ∈ [0, 0.2] ∪ [0.8, 1]; ρ ε (x, 0) = 1 + ε 2 , p ε (x, 0) = 1 x ∈ (0.2, 0.3]; ρ ε (x, 0) = 1, p ε (x, 0) = 1 + ε 2 /2 x ∈ (0.3, 0.7] ρ ε (x, 0) = 1 -ε 2 , p ε (x, 0) = 1 x ∈ (0.7, 0.8]
This example consists of several Riemann problems. Shocks and contact discontinuities are stronger when ε is bigger. We first check the difference of the three schemes (28), ( 30) and (31). The CFL condition for the linearized reduced stencil scheme (31) is discussed in (ii) and a fixed Courant number independent of ε is found numerically. Compared with the first order ICE method using local Lax-Friedrich discretization for (5), the improvement of removing nonphysical oscillation of our scheme is shown. We investigate the effect of α for different values of ε in (iii). In (iv), when α = 1, we numerically test the uniform convergence order. Finally, the AP property and its advantages are demonstrated in (v) by comparing with the fully explicit Lax-Fridrich scheme for the initial Isentropic Euler equation [START_REF] Bijl | A unified method for computing incompressible and compressible flows in boundary-fitted coordinates[END_REF]. When ε = 0.1, the initial density and momentum are displayed in Figure 1 and we can see the discontinuities clearly.

(i) In this example, we choose ε = 0.8, 0.3, 0.05 corresponding to the compressible, intermediate and incompressible regimes. The numerical results at T = 0.05 of "NL", "L" and "LD" are represented in Figure 2. Here ∆t is chosen to make all these three schemes stable and diminishing ∆t only will not improve much the numerical accuracy. The reference solution is calculated by an explicit Lax-Friedrich method [START_REF] Kurganov | New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations[END_REF][START_REF] Kurganov | Solution of two-dimensional Riemann problems for gas dynamics without Riemann problem solvers[END_REF] 1/20000. We can see that all these three methods can capture the right shock speed. The results of the three schemes are quite close, which implies that the linearization idea does simplify the scheme but the "LD" scheme does not really introduce less diffusion. When ε is small, though we can no longer capture all the details of the waves, the error is of the order ∆x which is the maximum information one can expect. Numerically, for different scales of ε, there is not much difference between these three methods. Thus in the following one dimensional examples, we only test the performance of the "LD" scheme.

(ii) Because of the explicit treatment of the flux terms in the momentum equation, the stability of the 'LD' scheme can be only guaranteed under the following CFL condition

∆t ≤ σ min i ∆x |u i | + αP ′ (ρ ε ) . ( 53
)
Here 0 < σ < 1 is the Courant number and is set up at initialization. Consistently with the fact that these three methods are AP, the Courant number does not depend on ε. Indeed, below, we numerically verify that σ is independent of ε. For ε = 0.8, 0.3, 0.05, the numerical Courant numbers are displayed in Table 1 and we can see numerically that the biggest allowed max{u} ∆t ∆x are close to 1 for all ε's. Therefore, σ = 0.9 is enough to guarantee stability and is numerically shown to be independent of ε. By contrast, the explicit local Lax-Friedrich scheme for the original Euler equation has a stability condition which becomes more and more restrictive as a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0. Here max{λ } denotes the maximum of max{λ j } defined in [START_REF] Munz | Linearized acoustic perturbation equations for low Mach number flow with variable density and temperature[END_REF] until T = 0.1 for all time steps.

(iii) The classical ICE method even in its conservative form introduces some nonphysical oscillations, no matter how small the time step is. These oscillations cannot be diminished by decreasing the time step. Their amplitude becomes smaller as the mesh is refined as long as the scheme is stable.

In this part we show that our method can suppress these oscillations numerically by choosing α = 1. When T = 0.01, for ε = 0.8, 0.3, 0.05, the numerical results of both our method with α = 1 and ICE calculated by ∆x = 1/200, ∆t = 1/20000 are displayed in Figure 3. The oscillations are more important for the ICE method and smooth away when α = 1. We can see that numerical nonphysical oscillations occur in the results of the ICE method when ε = 0.8, 0.3, but disappear when ε becomes small. This can also be seen from ( 47), (48). When ε is small the diffusion introduced by the implicitness is bigger. These oscillations also go away as time goes on due to dissipation.

(iv) When α = 1, the relative errors of the "LD" scheme for different ∆x, ∆t 2. Here ∆x, ∆t do not need to resolve ε and the reference solution is obtained by the explicit LLF scheme calculated with a very fine mesh ∆x = 1/1280, ∆t = 1/128000. We can see that good numerical approximations can be obtained without resolving the small ε. The convergence order is 1/2 when ∆t/∆x is fixed, uniformly with respect to ε. This convergence order when there are discontinuities is the same as the explicit LLF [START_REF] Leveque | Numerical methods for conservation laws[END_REF]. We can see from Table 2 that refinement in the time step does not improve the accuracy much (provided the Courant number is appropriately small, like σ = 0.7). Take ε = 0.8 as an example. When ∆x = 1/320, in order to obtain stability, ∆t should be less than 1/1920. It is demonstrated in Table 2 that the error calculated with ∆x = 1/320 does not decrease much when ∆t is changed from 1/2880 to 1/12800. Thus as long as the scheme is stable, we cannot use a smaller ∆t to obtain a better accuracy. This feature is the same as for standard hyperbolic solvers. resolved mesh ∆x = 1/2000, ∆t = 1/5000 are displayed in Figure 4, while the fully explicit Lax-Fridrich scheme is not stable with the same mesh size. We do capture the incompressible limit when ∆x, ∆t do not resolve ε.

Example 2: In this example we simulate the evolution of two collision acoustic waves by the "LD" scheme and test the convergence. We choose α = 1, ε = 0.1, ∆x = 1/100, ∆t = 1/1000. Here ∆t is chosen to stabilize the scheme and decreasing ∆t alone will not improve much the numerical accuracy. Similar to Klein's paper [START_REF] Klein | Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics I: one-dimensional flow[END_REF], P(ρ ε ) and the initial conditions are chosen as

P(ρ ε ) = ρ 1.4 ε , for x ∈ [-1, 1] ρ ε (x, 0) = 0.955 + ε 2 1 -cos(2πx) , u ε (x, 0) = -sign(x) √ 1.4 1 -cos(2πx) .
The initial density and momentum are displayed in Figure 5.

For ε = 0.1, the numerical results of the "LD" scheme at different times T are shown in Figure 6. The initial data approximate two acoustic pulses, one rightrunning and one left-running. They collide and their superposition gives rise to a maximum in the density. Then the pulses separate again. This procedure is demonstrated clearly in Figure 6.

Example 3

In this example, we show numerical results for the two dimensional case. Let P(ρ) = ρ 2 and the computational domain be (x, y) ∈ [0, 1] × [0, 1]. Be- cause no shock will form in this example, we choose α = 0 and the initial condition as follows:

ρ(x, y, 0) = 1 + ε 2 sin 2 (2π(x + y)), q 1 (x, y, 0) = sin(2π(x -y)) + ε 2 sin(2π(x + y)), q 2 (x, y, 0) = sin(2π(x -y)) + ε 2 cos(2π(x + y)).

The initial conditions for ε = 0.8 and numerical results at T = 1 calculated with ∆x = 1/20, ∆t = 1/80 are shown in Figure 7. Numerical tests show that a similar CFL condition is required as for the one-dimensional case. When ε = 0.05 at time T = 1, the numerical results with an unresolved mesh ∆x = 1/20, ∆t = 1/80 and a resolved mesh ∆x = 1/80, ∆t = 1/320 are displayed in Figure 8. We can see that the results using the coarse mesh are much 'smoother' than the one using the refined mesh. In this example the amplitude decay due to numerical diffusion cannot be ignored. When a coarse mesh is used, the first order method is known to have dissipation. This is mainly due to the numerical diffusion term, which smoothes out the solution. This phenomenon not only happens when ε is small but also when ε is O [START_REF] Bijl | A unified method for computing incompressible and compressible flows in boundary-fitted coordinates[END_REF]. We can also see from Figure 8 that when ε is small, D x p 1ε + D y p 2ε is close to 0. As a comparison, the numerical solutions of the incompressible limit (4) with and without numerical viscosity are shown in Figure 9. The latter is obtained by a difference method based on a staggered grid configuration [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF]. This staggered difference method is attractive for incompressible flows, since no artificial terms a) -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0. are needed to obtain stability and suppress the oscillations. Because of the stable pressure-velocity coupling, solutions with almost no viscosity can be obtained. The viscosity introduced here is of the form A 2 ∆x where A is given by (36). We can see that the amplitude of the wave decays as time evolves even though the viscosity is only O(∆x). In the limit of ε → 0, (42) generates a discretization of the incompressible limit with O(∆x) numerical diffusion terms. This is why the results for ∆x = 1/20, ∆t = 1/80 in Figure 8 are close to those with viscosity in Figure 9. When the meshes are refined, less diffusion is introduced and the solution becomes closer to the solution with no viscosity. The scheme indeed catches the incompressible Euler limit and good numerical approximations can be obtained without resolving ε, which confirms the AP property that is proved in section 4. However we need to take care of the numerical diffusion when coarse meshes are used. One possible way to improve this is to use less diffusive shock capturing schemes at first order or higher order schemes using the MUSCL strategy for instance [START_REF] Degond | Polynomial upwind schemes for hyperbolic systems[END_REF][START_REF] Kurganov | New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations[END_REF][START_REF] Kurganov | Solution of two-dimensional Riemann problems for gas dynamics without Riemann problem solvers[END_REF], or to use staggered grid discretizations.

Conclusion

We propose an all speed scheme for the Isentropic Euler equation. The key idea is the semi-implicit time discretization, in which the low Mach number stiff pressure term is divided into two parts, one being treated explicitly and the other one implicitly. Moreover, the flux of the density equation is also treated implicitly. The parameter which tunes the explicit-implicit decomposition of the pressure term allows to suppress the nonphysical oscillations. The numerical results show that the oscillations around shocks of O(ε 2 ) strength can be suppressed by choosing α = 1. The low Mach number limit of the time semi-discrete scheme becomes an elliptic equation for the pressure term, so that the density becomes a constant when ε → 0. In this way, the incompressible property is recovered in the limit ε → 0. Implemented with proper space discretizations, we can propose an AP scheme which can capture the incompressible limit without the need for ∆t, ∆x to resolve ε.

In this paper we demonstrate the potential of this idea by using the first order Lax-Friedrich scheme with local evaluation of the wave speeds. Though this first order method is quite dissipative, we can observe that the scheme is stable independently of ε and that the CFL condition is ∆t = O(∆x) uniformly in ε. It can also capture the right incompressible limit without resolving the mach number. Higher order space discretizations like the MUSCL method [START_REF] Degond | Polynomial upwind schemes for hyperbolic systems[END_REF][START_REF] Kurganov | New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations[END_REF][START_REF] Kurganov | Solution of two-dimensional Riemann problems for gas dynamics without Riemann problem solvers[END_REF]

Figure 1 :

 1 Figure 1: Example 1. When ε = 0.1, the initial density and momentum are displayed.

Figure 2 :

 2 Figure 2: Example 1. When T = 0.05, ∆x = 1/200, ∆t = 1/2000, the density and momentum of the "NL", "L" and "LD" schemes for isentropic Euler equation are represented respectively by dashed, dash dotted, and dotted lines. The solid line is the reference solution calculated by an explicit Lax-Friedrich method [19, 20] with ∆x = 1/500, ∆t = 1/20000. a): ε = 0.8; b): ε = 0.3; c): ε = 0.05. Left: density; Right: momentum. For all ε's, these three lines are so close to each other that '-.-.' and '...' are not visible in the figure.

ε

  goes to zero. Thus the CFL condition of the standard hyperbolic solver ∆t = O(ε∆x) is considerably improved.

  ε

Figure 3 :

 3 Figure 3: Example 1. When T = 0.01, the density and momentum for different ε are presented. The solid and dashed lines are the numerical results of our scheme and ICE with ∆x = 1/200, ∆t = 1/20000 respectively. a) ε = 0.8; b) ε = 0.3; c)ε = 0.05.

Figure 4 :

 4 Figure 4: Example 1. By using the "LD" scheme, the density (left) and momentum (right) for ε = 0.005 at T = 0.01 are represented. The circles are the results for ∆x = 1/20, ∆t = 1/500 and the solid line is calculated with ∆x = 1/2000, ∆t = 1/5000.

Figure 5 :

 5 Figure 5: Example 2. When ε = 0.1, the initial density and momentum are displayed.

Figure 6 :

 6 Figure 6: Example 2. When ε = 0.1, the density and momentum of the "LD" scheme at different times are represented: a) T = 0.01; b) T = 0.02; c) T = 0.04; d) T = 0.06. e) T = 0.08. All these pictures correspond to ∆x = 1/50, ∆t = 1/1000.

Figure 7 :Figure 8 :Figure 9 :

 789 scheme which can capture the incompressible limit without the need for ∆t, ∆x to resolve ε.In this paper we demonstrate the potential of this idea by using the first order Lax-Friedrich scheme with local evaluation of the wave speeds. Though this first order method is quite dissipative, we can observe that the scheme is stable independently of ε and that the CFL condition is ∆t = O(∆x) uniformly in ε. It can also capture the right incompressible limit without resolving the mach number. Higher order space discretizations like the MUSCL method[START_REF] Degond | Polynomial upwind schemes for hyperbolic systems[END_REF][START_REF] Kurganov | New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations[END_REF][START_REF] Kurganov | Solution of two-dimensional Riemann problems for gas dynamics without Riemann problem solvers[END_REF] can be built

  with ∆x = 1/500, ∆t =

		1.02												1.01
		1.015												1.008
														1.006
		1.01
														1.004
		1.005												1.002
	density	1											momentum	1
		0.995												0.998
														0.996
		0.99
														0.994
		0.985												0.992
		0.98	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	0.99	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
								x												x

Table 1 :

 1 Example 1. The numerical Courant numbers for different ε.

	0.8	max λ 4.24 1/100 1/340 3.40 1.25 ∆x stable∆t ∆x ∆t u ∆t ∆x
	0.8	6.35 1/200 1/970 4.85 1.31
	0.8	6.58 1/400 1/2420 6.05 1.09
	0.8	6.70 1/800 1/5460 6.82 0.982
	0.3	2.64 1/100 1/260 2.60 1.02
	0.3	2.70 1/200 1/510 2.55 1.06
	0.3	2.76 1/400 1/1000 2.50 1.10
	0.3	2.81 1/800 1/2050 2.56 1.10
	0.05 2.43 1/100 1/260 2.60 0.93
	0.05 2.44 1/200 1/490 2.45 1.00
	0.05 2.45 1/400 1/960 2.40 1.02
	0.05 2.46 1/800 1/1920 2.40 1.03

Table 2 :

 2 Example 1. T = 0.1, the L 2 norm of the relative error between the reference solution which is calculated with a very fine mesh ∆x = 1/1280, ∆t = 1/128000 and the numerical results for different ε with different ∆x, ∆t are displayed.

	ε 0.8 1/20 ∆x 0.8 1/40 0.8 1/80 0.8 1/160 1/1440 1.985 * 10 -1 1.75 2.048 * 10 -1 1.36 ∆t e(ρ ε ) ratio ratio e(p ε ) -1.197 -1/180 9.739 * 10 -1 1/360 5.959 * 10 -1 1.63 7.484 * 10 -1 1.16 1/720 3.467 * 10 -1 1.72 4.180 * 10 -1 1.31 0.8 1/320 1/2880 1.126 * 10 -1 1.76 8.477 * 10 -2 1.79 0.8 1/320 1/12800 1.126 * 10 -1 --8.539 * 10 -2 0.05 1/20 1/70 4.679 * 10 -3 --1.355 * 10 -1 0.05 1/40 1/140 3.305 * 10 -3 1.42 9.574 * 10 -2 1.42 0.05 1/80 1/280 2.353 * 10 -3 1.40 6.758 * 10 -2 1.42 0.05 1/160 1/560 1.655 * 10 -3 1.42 4.430 * 10 -2 1.53 0.05 1/320 1/1120 1.094 * 10 -3 1.51 2.538 * 10 -2 1.75 0.05 1/320 1/12800 6.012 * 10 -4 -9.303 * 10 -3 -
	(v) We emphasize the AP property in this final part. For ε = 0.005, the numer-ical results at T = 0.01 with unresolved mesh ∆x = 1/20, ∆t = 1/500 and

Acknowledgments

This work was supported by the french 'Commissariat à l'Energie Atomique (CEA)' (Centre de Saclay) in the frame of the contract 'ASTRE', # SAV 34 160.

into this framework. This is the subject of current work.

This paper provides a framework for the design of a class of all speed schemes. Compared with the ICE method [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF][START_REF] Harlow | A numerical fluid dynamics calculation method for all flow speeds[END_REF] and some recent work by Jin, Liu and Hauck [START_REF] Haack | All speed asymptotic preserving schemes for compressible flows[END_REF], the idea is simpler and more natural. This framework can also be easily extended to the full Euler equation and flows with variable densities and temperatures. These extensions and applications [START_REF] Munz | Linearized acoustic perturbation equations for low Mach number flow with variable density and temperature[END_REF] will be the subject of future work.