Many physical phenomena in hydrodynamics and magnetohydrodynamics are governed by the compressible Euler equations which represent the fundamental conservation principles of mass, momentum and energy. In many problems from meteorology, geophysics, combustion etc., one often encounters a situation where a characteristic fluid velocity is much lesser than a corresponding sound velocity in the medium. In such cases, the ratio of these two velocities typically gives rise to a singular perturbation parameter, such as the Mach number or the Froude number. When the motion of the fluid is humble in comparison to the motion of sound waves in the medium, i.e. when the Mach number tends to zero, the fluid can be treated as almost incompressible. From a mathematical point of view, one can say that solutions of the compressible Euler equations are close to those of their incompressible counter parts as the Mach number approaches zero. In many seminal works, e.g. [18, 19, 28], a rigorous convergence analysis of the solutions of the purely hyperbolic compressible Euler system to those of the mixed hyperbolic-elliptic incompressible Euler system in the limit of zero Mach number can be found.
It is widely known in the literature that standard explicit Gudunov-type finite volume compressible flow solvers suffer from a lot pathologies at low Mach numbers. The stiffness arising from stringent CFL stability restrictions, loss of accuracy due to the creation of spurious waves, lack of stability due to the dependence of numerical viscosity on the Mach number, and inability to respect the transitional behaviour of the continuous system of governing equations are some of the commonly observed ailments, to name but a few. The stiffness due to CFL restrictions severely imposes the timesteps to be as small as the order of the Mach number, resulting in a drastic slowdown of the solver in a low Mach number regime. The inaccuracy of the numerical solution, and inability of the numerical method to respect the transitional behaviour of the equations are mainly due to the loss of information in the passage from continuous to discrete level. We refer the reader to [10, 13, 20], and the references therein for a detailed account of the above-mentioned anomalies.
Design of stable and accurate numerical schemes for weakly compressible or low Mach number flows is a challenging task due to the aforementioned difficulties and more. In the literature, several approaches to develop numerical schemes which work for small as well as large Mach numbers, addressing one or more of the above issues, can be found. In , Bijl and Wesseling used the standard MAC-type finite difference method for the incompressible Euler equations to derive a scheme which can simulate flows at a wide range of Mach numbers. Another approach is due to Munz et al. 
, in which the wellknown SIMPLE method for incompressible flows is extended to weakly compressible flows using the insight gained from an asymptotic analysis. Recently, Feistauer and Kučera developed an all Mach number flow solver in the context of high order discontinuous Galerkin methods. We refer the interested reader also to, e.g. [3, 21, 29], for a different approach using the so-called sound-proof models which eliminate sound waves completely and provide good approximations for low speed atmospheric flows.
The convergence of solutions of the compressible Euler equations to those of the incompressible Euler equations in the limit of zero Mach number, and the associated challenges in numerical approximation is an active area of research. The consistency of the limit of a compressible flow solver, if it exists, with the incompressible limit system, and the stiffness arising from stringent stability requirements are usually addressed in the framework of the famous asymptotic preserving (AP) methodology. The notion of AP schemes was initially introduced by Jin  for kinetic transport equations; see also [8, 16] for a comprehensive review of this subject. The AP framework provides a systematic and robust tool for developing numerical discretisation techniques for singular perturbation problems. The procedure takes into account the transitional behaviour of the governing equations of the problem, and its stability requirements are independent of the singular perturbation parameter. Hence, an AP scheme for the compressible Euler equations automatically transforms to an incompressible solver when the Mach number goes to zero. As discussed in , semi-implicit time-stepping techniques provide a systematic approach to derive AP schemes; see, e.g. [5, 7, 9, 24, 30], for some semi-implicit AP schemes for the Euler or shallow water equations.
The presence of singular perturbation parameters gives rise to multiple scales in time as well as space. As mentioned above, a widely used strategy to resolve multiple timescales, and to derive an AP scheme is the use of semi-implicit time discretisations. The property of being AP of a numerical approximation is primarily that of the particular time discretisation used. Implicit-explicit Runge-Kutta (IMEX-RK) schemes 
offer a precise and robust approach to define high order semi-implicit AP schemes for singularly perturbed ordinary differential equations. These schemes have been extensively used and analysed, e.g. in[1, 25, 26]
for stiff systems of ordinary differential equations, and also time-dependent partial differential equations. The decisive step in the implementation of an IMEX-RK method to a system of conservation laws, such as the Euler or shallow water equations, is a splitting of the fluxes into the so-called stiff and non-stiff components in order that the resulting scheme is AP.
Recently, in [10, 11], the authors have presented the results of a detailed study on the accuracy of Godunov-type schemes at low Mach numbers. It has been shown in these papers that the origin of inaccuracies can be attributed to several factors: the number of space dimensions, particular discretisation strategies and numerical diffusion in the scheme used, cell geometry chosen, and so on. Based on the work of Schochet , Dellacherie in 
has presented the results of an extensive study on the accuracy of Gudunov-type schemes at low Mach numbers on Cartesian meshes. At the core of this analysis is an estimate from which states that when the Mach number is small, a solution to the compressible Euler equations is close to that of the incompressible equations whenever the corresponding initial data are close, cf. also [18, 19]. Further, an analysis of the nature of solutions of the linear wave equation system at low Mach numbers using Hodge decompositions showed that a sufficient condition to ensure the above-mentioned estimate in the linear case is the invariance of a so-called ‘well-prepared space’ of constant densities and divergence-free velocities. In the light of this invariance property, Dellacherie established that the origin of inaccuracies of standard Godunov-type schemes in two and three-dimensional geometries using the notion of first-order modified equations. Finally, a sufficient condition to avoid the generation of spurious waves, and to ensure accuracy at low Mach numbers for nonlinear schemes for the Euler equations is proposed by modifying the numerical viscosity, and changing the discretisation of the pressure gradient term.
In this paper, we study a class of second order accurate finite volume schemes for the isentropic Euler equations, employing the IMEX-RK time-stepping procedure. We prove that both the time semi-discrete and space-time fully discrete schemes are asymptotically consistent with the limiting incompressible system in the zero Mach number limit, and that their stability constraints are independent of the Mach number which implies the AP property. We also define a notion of accuracy at low Mach numbers, hereafter designated as the asymptotic accuracy, in terms of the invariance of the well-prepared space as in . It is shown theoretically as well as numerically that both the time semi-discrete and space-time fully discrete schemes are asymptotically accurate. The rest of this paper is organised as follows. In Section 2 we mention the zero Mach number limits of the isentropic Euler and the linear wave equation systems, and recall the main convergence results from . Section 3 is dedicated to the defintion of an AP scheme which is also accurate low Mach numbers. In Section 4 we introduce the IMEX-RK time semi-discrete schemes, and in Section 5 we show their AP property and asymptotic accuracy. Space-time fully-discrete schemes, and their analysis are presented in Section 6. In Section 7 numerical evidences are provided to support the theoretical claims made. Finally, we close the paper with some concluding remarks in Section 8.
2. Zero Mach Number Limits of the Euler and the Wave Equation Systems
2.1. Zero Mach Number Limit of the Euler System
The scaled, compressible Euler equations are given by
where the independent variables are time and space , and the dependent variables are the density and velocity of the fluid . In order to close the system (2.1), an isentropic equation of state , where is a positive constant, is assumed. Here, the parameter is the reference Mach number defined by
with and being a reference fluid speed, and a reference sound speed, respectively.
for all dependent variables. Doing a scale analysis, and using appropriate boundary conditions, cf. e.g. , yields the multiscale representation
for the density and the incompressible Euler system:
for the velocity . Here, the second order pressure survives as the incompressible pressure.
Loosely speaking, the results of [18, 19] shows that when the initial data are almost incompressible, the solutions of the compressible Euler equations (2.1) are good approximations of those of the incompressible Euler equations (2.5). However, in order to make a precise statement, we recall the following convergence result due to Schochet . Note that the non-dimensionalised, isentropic Euler system (2.1) can be recast in a non-conservative, evolution form as
Here, is the convective operator with a time scale of order and is the acoustic operator with a time scale of the order of . The system (2.6) is supplied with periodic boundary conditions
where denotes the -dimensional torus.
Note that is the subspace of spatially constant densities and divergence-free velocities. In other words, it is the ‘incompressible’ subspace and hence, hereafter, is referred to as the ‘well-prepared’ space. Having defined and , the following Helmholtz-Leray decomposition can be immediately written.
Thus, for any , there exists a unique and , such that . Let us also define the Helmholtz-Leray projection of any onto the space as
Let be a solution of the initial value problem for the system (2.6), i.e.
and let be a solution of the projected subsystem:
where . Then, for , there holds the estimate:
It is easy to see that the projected subsystem (2.15)-(2.16) is nothing but the incompressible Euler system (2.5). The essence of the above theorem is that when the initial data are well prepared, i.e. they are taken in , solutions of the compressible Euler system (2.1) approximate those of the incompressible Euler system (2.5) when . The theorem, clearly, is in agreement with the results obtained in [18, 19].
2.2. Zero Mach Number Limit of the Linear Wave Equation System
The rest of this section is devoted to the analysis of the linear wave equation system as a prototype model of the compressible Euler equations (2.6). Linearising the system of equations (2.6) about a constant state , we get the following wave equation system with advection:
In order to present an analysis of the solutions of (2.18), we first introduce a few definitions. For any and in , we consider the innerproduct
where the innerproducts appearing on the right hand side are the usual innerproducts in . The norm generated by the above inner product is given by
Clearly, the norm is equivalent to the usual norm in . We also consider the energies
which are, respectively, the total energy, the incompressible energy, and the acoustic or compressible energy. Since , we have . In the following, we state the corresponding linearised version of Schochet’s result, i.e. Theorem 2.1.
Theorem 2.3 ().
Let be a solution of the IVP for the wave equation system (2.18), i.e.
and let be a solution of the IVP:
where . Let be the Helmholtz-Leray decomposition of . Then the following holds.
is the solution of (2.23) with initial condition .
Moreover, there holds the energy conservation:
and as a consequence, the following estimate holds:
for the unknowns . Since the boundary conditions are periodic, and (2.29) is linear, the divergence-free condition (2.30) on now forces to be a constant. Hence, , and (2.29) then reduces to a linear transport equation for , and the divergence-free condition (2.30) is now required only at time . Thus, the projected subsystem (2.25)-(2.26) and the system (2.29)-(2.30) are equivalent, and are incompressible.
The above theorem lies at the centre of the analysis of numerical schemes presented in . However, numerical discretisation procedures introduce numerical diffusion, dispersion or higher order derivative terms, depending on the order of the scheme under consideration. Hence, the Theorem 2.3 has to be generalised to accomodate more general spatial differential operators arising from the discretisation. Such a generalisation is presented in  which also gives a sufficient condition to ensure the estimate (2.28).
Theorem 2.5 ().
Let be a solution of the IVP:
which is assumed to be well-posed in , with a linear spatial differential operator. Then the following conclusions hold.
3. Asymptotic Preserving and Asymptotically Accurate Schemes
As mentioned in the introduction, the AP property must be necessary for any numerical scheme to survive in the passage of the limit . However, the accuracy of the scheme also deteriorates in the low Mach number limit, and hence it is essential for the scheme to maintain its accuracy. The two primary aims of the present work are to design and test an IMEX-RK finite volume scheme which is AP, and to maintain the accuracy in the limit . In what follows, we formally state the definitions of the AP property and asymptotic accuracy.
3.1. AP Property
Numerical solution of singular perturbation problems is, in general, a difficult task mainly because of the existence of multiple scales, the dependence of the stability characteristics on the perturbation parameter, and the transitional nature of the governing equations. A numerical approximation scheme may not resolve the scales consistently, its stability requirements might deteriorate, and in the singular limit the scheme might approximate a completely different set of equations than the actual limiting system. AP schemes provide a robust approach as a remedy against these problems. AP schemes were originally developed in the context of kinetic transport equations; see [8, 15, 16], and the references therein for more details.
Let denote a singularly perturbed problem with the perturbation parameter . Let denote the limiting system of for . A discretisation of , with being the discretisation parameter, is called AP, if
is a consistent discretisation of the problem , called the asymptotic consistency, and
the stability constraints on are independent of , called the asymptotic stability.
In other words, the following diagram commutes:
3.2. Asymptotic Accuracy
It has to be noted that the Theorem 2.5, and its conclusions are valid for a linear scheme for a linear system of equations. Nonetheless, motivated by this result, and analogous considerations from , we propose the asymptotic accuracy of a scheme to be its ability to preserve the estimate (i) in Theorem 2.5. In other words, a numerical scheme for the compressible Euler system (2.1) should possess the ability to maintain solutions close to the incompressible solutions in for all times, whenever the initial data are well-prepared. In , such a notion was designated to be a scheme being accurate at low Mach numbers. Since the invariance of the well prepared space is a sufficient condition to maintain the estimate in the linear case, we make the following definition.
A numerical approximation for the compressible Euler system (2.1) is said to be asymptotically accurate (AA), if it leaves the incompressible subspace invariant.
It has to be noted both the AP and AA properties are to be satisfied by the time semi-discrete as well as space-time fully-discrete schemes. In the following sections we establish that the IMEX-RK schemes considered in this paper possess both the AP and AA properties.
4. Time Semi-discrete Scheme
4.1. Implicit-Explicit (IMEX) Runge-Kutta (RK) Time Discretisation
The IMEX-RK schemes are designed for the numerical integration of stiff systems of ordinary differential equations (ODEs) of the form
where , and is usually known as the stiffness parameter. The functions and are called, respectively, the non-stiff part and the stiff part of the system (4.1). Stiff systems of ODEs of the type (4.1) typically arise in the modelling of semiconductor devices, study of kinetic equations, theory of hyperbolic systems with relaxation etc.; see, e.g. , for a comprehensive treatment of such systems.
The matrices and are matrices such that resulting scheme is explicit in and implicit in . However, in order to reduce the number of implicit evaluations in the intermediate stages (4.2), we consider only the so-called diagonally implicit Runge-Kutta (DIRK) schemes in which for , and for . The coefficients and and the weights and are fixed by the order conditions; see [25, 26] for details. For the sake of completion, and to draw reference for the analysis carried out later, we write down the order conditions for a first order and second order IMEX-RK scheme as follows:
The conditions (4.4), (4.5) and (4.6) are, respectively, the consistency, the first order, and the second order order conditions. To further simplify the analysis of the schemes presented in this paper, we restrict ourselves only to two types of DIRK schemes, namely the type-A and type-CK schemes which are defined below; see [17, 26] for details.
An IMEX-RK scheme with the Butcher tableau given in Figure 1 is said to be of
type-A, if the matrix is invertible;
type-CK, if the matrix , can be written as
where and is invertible.
In the results presented in the later sections, we have used both the first order Euler(1,1,1), and second order ARS(2,2,2) schemes for time discretisations; see Figure 2 for their Butcher tableaus. Here, in the triplet , is the number of stages of the implicit part, the number gives the number of stages for the explicit part and gives the overall order of the scheme. For a detailed study of IMEX-RK schemes we refer the reader to [14, 17, 25, 26] and the references therein.
4.2. IMEX-RK Time Discretisation of the Euler System
where is the momentum. Applying IMEX-RK time discretisation, i.e. treating explicitly and implicitly, results in the following semi-discrete scheme.
The stage of an -stage IMEX-RK scheme for the Euler system (2.1) is defined as
The numerical solution at time is given by
In the above, and throughout the rest of this paper, we follow the convention that a repeated index always denotes the summation with respect to that index. The indices and are used to denote, respectively, the summation in the explicit and implicit terms, i.e. they assume values in the sets and .
Though the intermediate stages (4.8)-(4.9) consist of two implicit steps, the resolution of the scheme (4.8)-(4.9) is quite easy. As in , we eliminate between (4.8) and (4.9), and obtain the nonlinear elliptic equation:
for . Here, we have denoted by
those terms that can be explicitly evaluated. Once is known after solving the elliptic equation (4.12), can be explicitly evaluated from (4.9). Finally, the updates and are calculated explicitly using (4.10)-(4.11) with the values obtained from the intermediate stages.
In order to further reduce the nonlinear nature of the elliptic equation (4.12), and to make the numerical implementation simpler, in our computations we transform the nonlinear elliptic equation for into a semilinear elliptic equation for the pressure .
5. Analysis of Time Semi-discrete Scheme
5.1. AP Property
Establishing the AP property consists of showing the consistency of the scheme with that of the incompressible system as with its stability requirements independent of . In order to show the former, we perform an asymptotic analysis of the time discrete scheme (4.8)-(4.11), and to show the latter, we use a linear stability analysis.
5.1.1. Asymptotic Consistency
Suppose that the data at time are well-prepared, i.e. and have the form:
where and . Then at the intermediate stages, and defined by (4.8) and (4.9) satisfy and , i.e. they are well-prepared as well, and as a consequence of this the numerical solutions and are well-prepared. In other words, the semi-discrete scheme (4.8)-(4.11) is asymptotically consistent with the incompressible limit system in the sense of Definition 3.1.
We use induction on the number of stages to prove the theorem. To this end, we prove the consistency for the first stage, i.e. , which corresponds to a fully implicit step. We plugin the ansatz (2.3) for each of the dependent variables in (4.8)-(4.11). Equating to zero the terms in system (4.8)-(4.9) yields
Hence, the zeroth order density is spatially constant. In an analogous way we can show that the first order density is also spatially constant.
Next, we consider the terms in the mass update (4.8) to obtain
We Integrate the above equation (5.4) over a spatial domain and use the Gauss’ divergence theorem to get
If the boundary conditions of the problem are periodic or wall, then the right most integral in the above equation (5.5) vanishes; yielding . As the zeroth order density is a constant, as well as are constants. Using this in equation (5.4), we obtain the divergence condition:
for the leading order velocity . This completes the proof for .
Now proceeding as in the case of , we can obtain . Hence, for the stage the zeroth order density is same as that of , and in turn we also obtain the divergence condition:
Finally, to prove that the numerical solution at time is also well prepared, we consider the terms in the elliptic equation (4.12) and get
As , and for all , we have
Considering the terms in the final mass update (4.10) gives
Since and , we immediately get from (5.12) that
and similarly for . To get the divergence-free condition for the zeroth order velocity of the numerical solution, we equate to zero the terms in the final velocity update (4.11), and take its divergence to get
Note that the above proof is valid for any -stage IMEX-RK scheme. However, it uses the fact that for , which is the property of an RK scheme of type-A. The result also holds for a type-CK scheme in which the first step is trivial, and then on the proof follows similar lines as that of a type-A scheme.
5.1.2. Linearised -stability and -invariance
In this section, we analyse the correction terms arising from the time discretisations, and their effect on the asymptotic accuracy and stability of the resulting scheme. In order to analyse these correction terms, we follow the standard modified partial differential equation approach; see also  for a related analysis on stability. It has been shown in  that explicit Godunov-type schemes in a low Mach number regime are accurate only in one dimension, and that they are inaccurate in two and three dimensions. Specifically, the analysis in  shows that there is an acoustic time scale so that a Godunov-type scheme suffers from inaccuracies in a time interval of . A cure proposed in  is to remove the numerical diffusion terms from the momentum updates, and to use central differencing for the pressure gradient term. However, diffusion terms are required due to stability constraints and deleting them may not be a feasible option always; see also  for a related discussion. In the following, we show that using IMEX-RK time discretisation has the following advantages: no inaccuracy arises in two and three dimensions, it is not necessary to delete the diffusion coefficients to maintain stability, and the time-steps are independent of .
In order to establish the asymptotic stability of the IMEX-RK time discretisation for the isentropic Euler system (2.1), we carry out a thorough linear stability analysis for the linear wave equation system (2.18) as a linearised model. We believe that the results of linear stability analysis holds good also for the nonlinear Euler system (2.1). Analogous to Section 4, an IMEX-RK time discrete scheme for the linear wave equation system (2.18) is defined as follows.
The stage of an -stage IMEX-RK scheme for the wave equation system (2.18) is given by
The numerical solutions and at time are defined as
In the following, we derive the modified partial differential equation (MPDE) corresponding to a general second order accurate time discrete scheme of the form (5.17)-(5.20). We show that the solution of the MPDE is energy-dissipative, and the time-steps are independent of . Hence, the time discrete scheme (4.8)-(4.11) for the Euler system is linearly asymptotically -stable. In addition, we also prove that the MPDE keeps the well-prepared space invariant.
We consider a general second order accurate IMEX-RK time discretisation in (5.17)-(5.20). Expanding the unknown functions in a Taylor series, and making use of the conditions (4.4)-(4.6) for the IMEX-RK coefficients, results in the following MPDE:
Here, the operators and contain the third and fourth derivatives of the unknown functions and , respectively. Since these expressions are quite long, we present them only in the Appendix. After using the second-order order conditions (4.4)-(4.6), the MPDE (5.21) is free from first and second order derivative terms.
The rate of change of the energy , defined in (2.22), is given by
We use the MPDE (5.21) in (5.22), and integrate the resulting terms by parts. It is easy to see that the third order derivatives do not contribute as they all vanish in view of periodic boundary conditions, and hence we get