Artificial diffusion for convective and acoustic low Mach number flows I: Analysis of the modified equations, and application to Roe-type schemes

Three asymptotic limits exist for the Euler equations at low Mach number - purely convective, purely acoustic, and mixed convective-acoustic. Standard collocated density-based numerical schemes for compressible flow are known to fail at low Mach number due to the incorrect asymptotic scaling of the artificial diffusion. Previous studies of this class of schemes have shown a variety of behaviours across the different limits and proposed guidelines for the design of low-Mach schemes. However, these studies have primarily focused on specific discretisations and/or only the convective limit. In this paper, we review the low-Mach behaviour using the modified equations - the continuous Euler equations augmented with artificial diffusion terms - which are representative of a wide range of schemes in this class. By considering both convective and acoustic effects, we show that three diffusion scalings naturally arise. Single- and multiple-scale asymptotic analysis of these scalings shows that many of the important low-Mach features of this class of schemes can be reproduced in a straightforward manner in the continuous setting. As an example, we show that many existing low-Mach Roe-type finite-volume schemes match one of these three scalings. Our analysis corroborates previous analysis of these schemes, and we are able to refine previous guidelines on the design of low-Mach schemes by including both convective and acoustic effects. Discrete analysis and numerical examples demonstrate the behaviour of minimal Roe-type schemes with each of the three scalings for convective, acoustic, and mixed flows.


Introduction
Low Mach number flows with non-zero divergence -distinct from properly incompressible divergence-free flows -are of interest in a wide range of applications. Flows with significant heat transfer or external body forces experience compressibility and buoyancy effects, for example environmental flows [7] and cooling systems for electronic devices and turbomachinery [63,11]. Low Mach number acoustic phenomena are also of interest [28], for example noise prediction for wind turbines or aeroplanes during take-off and landing. Some flows experience both of these effects, such as low Mach number combustion, where heatrelease induced compressibility and acoustic interactions both play important roles in governing the flow evolution [62,13]. Some flows contain regions of both low and high Mach number in the same domain, for example aerofoils at high angle of attack where the background flow has low Mach number but regions of higher Mach numbers occur in the vicinity of the aerofoil [50]. To properly simulate any of these cases requires a numerical scheme which is accurate at low Mach number.
There are three common approaches for designing numerical methods for low Mach number flow. The first uses asymptotic analysis and other assumptions to simplify the full governing equations (either Euler or Navier-Stokes) into a reduced set of equations which can then be solved using specialised solution strategies [48], such as the Boussinesq buoyancy equations [51,43] or multiple pressure variable (MPV) methods [31,57]. This approach leads to schemes which are highly specialised and often very efficient, however they are restricted to the regimes where the simplifying assumptions hold. The second approach extends pressure based methods for incompressible flows -such as SIMPLE or fractional step -to the low Mach number regime. These methods have proved particularly useful for low Mach number aeroacoustics [77,79]. These methods can leverage the large existing research and code bases for incompressible methods, however they may only handle mixed Mach number flow with moderate, not high Mach number. The third approach is the one considered in this paper: extending density based methods for compressible flows into the low Mach number regime. These methods face efficiency and accuracy problems in this regime, but Turkel 1993 [70] states three reasons justifying the effort placed into overcoming these issues which are still valid today. Firstly, these methods can handle mixed Mach number flows, a major advantage over the other approaches. Secondly, these methods handle very naturally flows with large thermal effects which cause large density variations and strong coupling of the energy equation. Lastly, there is a large body of research and software built upon these methods which can be taken advantage of, including effective acceleration techniques and handling of complex geometries [3].
Collocated density based compressible numerical methods face two substantial problems in the low Mach number regime: efficiency and accuracy. The efficiency problem is easily understood by considering the wavespeeds of the Euler equations: as M = u/a → 0, the condition number (u + a)/u → 1/M , and the convective wavespeed u evolves very slowly relative to the acoustic wavespeed a. If the step size of an iterative scheme (for time-accurate evolution or for steady-state convergence) is chosen according to the acoustic wavespeed, then a large number of steps are required to resolve the convective features. The common solutions to this problem are semi-implicit methods which remove the acoustic CFL limit [77], or low Mach preconditioners which improve the condition number and restore efficiency [71,69]. The design and analysis of low Mach number preconditioning is a substantial research topic in itself, but is not the focus of the current work -although the impact of the artificial diffusion on efficiency will be discussed. See Turkel 1999 [69] for a review of low Mach number preconditioning.
The accuracy problem is less easily demonstrated, and is the main focus of this paper. Using asymptotic analysis, three limits for the Euler equations can be found at low Mach number. The limit most often of interest is the convective single-scale limit, where only convective features exist, and solutions approach the incompressible limit as M → 0 [30]. At the acoustic single-scale limit, only acoustic features exist, and the solutions approach those of the wave equations for linear acoustics. Finally, a single-space-scale multipletime-scale (or multiple-space-scale single-time-scale) analysis produces a mixed convective-acoustic limit with convective features on a slow timescale (short space-scale) co-existing with acoustic features on a fast timescale (long space-scale) [48]. The accuracy problem describes the inability of some schemes to produce solutions which match those of the desired limit (usually the convective limit).
Volpe 1993 [76] and Godfrey et al. 1993 [21] showed that a naïve application of transonic numerical schemes to flows at the convective low Mach number limit can produce inaccurate, physically inconsistent results. Godfrey et al. showed that calculating the artificial diffusion based on the preconditioned system instead of the original Euler system produces vastly improved results. Since then, many papers have been published presenting different methods for modifying the artificial diffusion to achieve accurate low Mach number results. These use a variety of strategies, building on central schemes [70,74,73,75], flux-difference [21,78,26,25,23,44,65,15,53,49,50,5,32,35,34,33,20], and flux-vector splitting [18,39,37,61,58,36] methods. The majority of studies consider only the convective limit, although a number also address flows with acoustic features [74,73,75,50,5,61,58,46]. There have been a number of excellent review and analysis articles published on artificial diffusion at low Mach number, covering the discrete equations of Roe-type schemes at the different limits [26,25,23,33,5], the modified equations [15], and the relationship to preconditioning [46,74,73,75]. We cover the literature on low Mach number artificial diffusion schemes and their analysis in more detail in section 2.
In this paper we review the role of the artificial diffusion in the class of collocated, density based schemes at low Mach number. We have endeavoured to find a formulation of the problem which is as simple as possible whilst still demonstrating the most important findings across this research area. We hope that this approach will help highlight the root causes of these findings, and enable it to act as an introduction for those with little experience in this field.
To achieve this, we consider all three low Mach number limits -the purely convective, purely acoustic, and mixed convective-acoustic -which are introduced in section 1. After this, we can review the previous literature in section 2. Previous studies usually begin by showing that a classical transonic scheme is inaccurate for the convective limit, then proceed to 'fix' this scheme. Instead, in section 3 we use the continuous modified equations in the entropy variables to design artificial diffusion schemes which naturally match each of the three limits independently of the specific discretisation. In section 4 we apply single-and multiple-scale asymptotic expansions to the modified equations of each scheme. In section 5 we transform the artificial diffusion to the conservative variables to find the equivalent finite-volume flux functions. We compare to previous works, and proceed to repeat the asymptotic expansions on the discrete equations. Finally, in section 6 we show a series of numerical examples which verify the results of the previous sections. In a sequel paper † we extend this to three families of convection-pressure flux splittings: AUSM [40], Zha-Bilgen [80] and Toro-Vasquez [68].

Asymptotic expansion of the Euler equations at low Mach number
We review the theoretical results for the asymptotic solutions of the Euler equations with an ideal gas at low Mach number using single-and multiple-timescale analysis. We present only the most relevant results, and work in the entropy variables w = (p, u, v, s) T . For a more extensive analysis in the conserved variables, see Müller 1999 [48], and for a multiple space scale analysis see [31]. First, all variables are non-dimensionalised: where tildes indicate local dimensional quantities, ∞ indicates reference dimensional quantities, and underlining indicates vector-valued quantities. E and H are the specific total energy and enthalpy respectively. The Euler equations in entropy variables become: ∂ t p + u · ∇p + γp∇ · u = 0 (2a) ρ∂ t u + M −2 ∇p + ρu · ∇u = 0 (2b) ∂ t s + u · ∇s = 0 (2c) Where M = √ γu∞ a∞ is the reference Mach number. The relation γp = ρa 2 has been used in the 3rd term in (2a), which is only valid for a perfect gas or barotropic gases which are polytropic (including isothermal or isentropic gases). It is not valid for a general barotropic gas, although the same results can still be obtained in this case. The entropy equation (2c) is not required for barotropic gases.

Single timescale convective limit
Using M as a small parameter, we treat the system (2) as a perturbation problem and expand all variables as power series of M : ψ(x, t, M ) = ψ (0) (x, t) + M ψ (1)   Relations (4a) and (4b) mean that the zeroth and first order pressure terms vary only in time. If p (1) is zero at all points on the boundary at all times then the pressure expansion (3) can be replaced with: The spatial variations of the non-dimensional convective pressure are O(M 2 ), so the dimensional pressure variations are O(ρu 2 ), i.e. independent of the Mach number. For a barotropic gas, this implies that the density also has O(M 2 ) spatial variations. The last relation (4c) is the evolution equation for the convective velocity. Expanding the pressure equation (2a) and making use of relations (4a) and (4b), we find the following O(M 0 ) and O(M 1 ) relations: The relations (5a) and (5b) imply that the divergence of the zeroth and first order velocities are spatially uniform and react everywhere instantaneously to temporal variations in the background pressure. Lastly, we expand the entropy equation. Every order n of the expansion is identical to the original equation, showing that the entropy is simply convected.
Note that there are no acoustic effects in this single-scale expansion, which should not be surprising given our choice of the convective timescale L ∞ /u ∞ to non-dimensionalise the time derivatives in equation (2).

Two timescale convective-acoustic limit
A two-timescale, one space scale asymptotic analysis can be used to include acoustic effects. Defining an additional non-dimensional time τ using the acoustic speed: leads to the power series expansion: The time derivatives at constant x and M are now: As for the convective limit, we start with the three lowest order relations from the velocity equation: The zeroth order pressure is again constant in space, however the same is no longer true for the first order pressure term, which now varies with zeroth order velocity fluctuations on the acoustic timescale. The second order pressure is still the relevant pressure for the zeroth order velocity fluctuations on the convective timescale (10c). Next, the pressure equation is expanded using (8) and making use of the relation (10a): Relations (10a) and (11a) imply that the background pressure varies only on the convective timescale. Relation (11b) shows that the leading order velocity divergence now has spatial variations on the acoustic timescale related to first order pressure variations. In fact, relations (10b) and (11b) are the equations for linear acoustics at low Mach number, with a source term d t p (0) (see [48] for more details). The relation (10c) is again the equation for the convective velocity variations. The correct expansion for the pressure at the two-timescale limit is (8), except the zeroth order term which becomes p (0) (x, t, τ ) = p (0) (t). Note that a purely acoustic limit of (2) can be found with a single-timescale expansion with the acoustic timescale τ , which results in the relations (10a,10b,11a), and (11b) without the d t p (0) term. We now expand the entropy equation using the two timescale expansion. This time we find: From which we see that the zeroth order entropy is constant on the acoustic timescale, which is consistent with a leading order approximation of isentropic sound waves.
To summarise the important points of this section, we have briefly covered two asymptotic limits of the inviscid homogeneous Euler equations at low Mach number. Firstly, a convective limit with spatially uniform velocity divergence and O(M 2 ) spatial pressure variations associated with convective features on the slow timescale t. Secondly, a mixed convective-acoustic limit which has spatially varying divergence and O(M ) spatial pressure variations associated with acoustic features on the fast timescale τ in addition to the convective features on the slow timescale. A single-scale acoustic limit can be constructed using only the fast timescale τ , which contains only the acoustic variations on the fast timescale.

Previous work
As highlighted in the introduction, Godfrey et al. showed in 1993 that calculating the artificial diffusion for a flux-difference-splitting scheme based on the preconditioned system produced accurate results for steady convective flows, whereas calculating the diffusion from the unmodified system produced physically inconsistent results [21]. This approach was adopted by a number of other researchers for both fluxdifference-splitting methods [78,26] and central difference methods [69,9,74]. Turkel et al 1994 andTurkel 1999 [71, 69] used the modified equations and the known low Mach number scalings for the convective limit to show that the artificial diffusion of a standard transonic scheme is poorly balanced for this limit, having velocity diffusion which is too high and pressure diffusion which is too low. They also showed that the preconditioned artificial diffusion rectifies both of these issues, having pressure and velocity diffusion terms which are larger and smaller by a factor of M respectively compared to the unpreconditioned scheme. We note that artificial diffusion is not the only approach to achieving stability at low Mach number. For example the kinetic energy consistent scheme of Subbareddy & Candler 2009 [64] is accurate and stable at low Mach number without artificial diffusion. However, artificial diffusion is the predominant method for stabilising collocated schemes in many applications.
Guillard & Viozat 1999 [26] used asymptotic expansions of the discrete equations to show that the unmodified Roe scheme [56] produces discrete solutions which do not match the convective asymptotic solutions of the exact Euler Equations, in that they contain spurious spatial oscillations of the first order pressure p (1) . In a later paper [25] they showed that this is because the diffusion of the unmodified scheme is derived from the acoustic component of the Euler equations, so produces a residual in the acoustic pressure p (1) even for initial conditions containing purely convective variations. Asymptotic expansion of the discrete equations of the preconditioned Roe scheme shows that this scheme produces discrete solutions which do match the exact asymptotic convective limit -the smaller velocity diffusion preventing the creation of spurious p (1) modes, and the larger pressure diffusion preventing chequerboard modes on p (1) .
The preconditioned diffusion proved to be popular and successful, especially because it can be used in combination with preconditioned timestepping, resulting in a scheme which is stable under a convective CFL condition, eliminating the efficiency problem. However, Birken & Meister 2005 [2] showed that for time accurate timestepping the preconditioned diffusion has a timestep limit of ∆t ∼ O(M 2 ) for stability of an explicit scheme, compared to ∆t ∼ O(M ) for the unmodified scheme, so requires implicit timestepping to be practical for time-accurate simulations.
Following the finding of Birken & Meister, a number of papers proposed a different approach to creating a diffusion scheme which is accurate at low Mach number. Thornber et al. 2008 [66] showed that the inaccuracies of standard Godunov schemes at low Mach number are linked to spurious entropy generation by the velocity diffusion, and subsequently proposed a modification to Roe's flux which reduced the velocity diffusion by a factor of M , resulting in accurate solutions at low Mach number [67]. Other methods of achieving this reduction in the velocity diffusion have been proposed, including those by Thornber et al. [65], Dellacherie [15], Rieper (LMRoe) [53], and Miczek et al. [47]. All of these methods retain the ∆t ∼ O(M ) stability limit of the original scheme, but with improved low-Mach accuracy.
Dellacherie 2010 and Dellacherie et al 2016 [15,14] used the modified equations for the linear acoustic waves with artificial diffusion to rigorously show that reducing the velocity diffusion of a standard Godunov scheme by at least a factor of M will prevent the production of spurious acoustic p (1) modes and ensure accuracy for the convective low Mach number limit, provided the initial conditions are wellprepared † . This shows that the increased pressure diffusion in the preconditioned artificial diffusion is in fact not necessary for low-Mach accuracy. However, it has been shown [15,23] that the increased pressure diffusion introduces a Brezzi-Pitkäranta type stabilisation [4], which prevents pressure chequerboard instabilities on collocated schemes, and has previously been applied to finite-element [27] and finite-volume † Well-prepared initial conditions are solutions which are close to the solution space of the convective limit i.e. contain vanishingly small acoustic components. See Schochet [59] and Dellacherie et al [15,14] for a formal definition.
A third approach, the all-speed Roe scheme, introduced by Li & Gu 2008 [32,35] reduces not only the velocity diffusion by a factor of M , but also the pressure diffusion compared to the standard Roe scheme, so that the diffusion is essentially a diagonal upwinding on the convective velocity scale. The reduction in velocity diffusion means that this scheme is accurate for convective low Mach number flow, but is susceptible to severe chequerboard instabilities. Li & Gu resolve this by reintroducing a pressure diffusion term into the physical mass flux using a momentum interpolation method similar to that proposed by Rhie & Chow [52], and used by Mary and Sagaut [44] and in AUSM schemes [40,38,37]. The scheme was later extended to unsteady flows by adding a timestep dependence to the pressure diffusion [34].
The discrete forms of these three approaches are reviewed and compared by Li & Gu 2013 [33] and Guillard & Nkonga 2017 [23]. Li & Gu compare the scaling of the coefficients of the pressure and velocity diffusion terms, and propose a set of guidelines on these scalings for accuracy and stability. Guillard & Nkonga review the origin of the accuracy problem with respect to the multiple possible low Mach number limits, and also highlight the dependence on grid type of the various low-Mach schemes. For a first order Roe scheme on a simplex mesh (triangles in 2D or tetrahedrons in 3D), the velocity degrees of freedom are reduced such that the jumps in the normal velocity over cell interfaces vanish. This eliminates the velocity diffusion, meaning that the unmodified Roe scheme actually provides accurate results for this cell geometry, despite its inaccuracy on other meshes. See [54,55,24,17] and section 3.3 of [23] for details of this particular behaviour.
The literature covered so far has focused on the convective limit. A number of papers have also covered schemes for flows with acoustic effects. In a series of papers, Venkateswaran and Merkle showed that preconditioned dual-time schemes are ill conditioned at the acoustic limit and display poor convergence [72], and that the preconditioned diffusion is inaccurate for flow with acoustics, having excessive pressure dissipation, whereas the unmodified scheme is both efficient and accurate for purely acoustic flow [74,46]. They developed an 'enhanced' diffusion scheme having the same pressure diffusion as the unmodified scheme but reduced velocity diffusion, as in the preconditioned scheme [73,75]. This diffusion scheme was shown to produce accurate results with reasonable convergence for both convective and acoustic flows. The later low-Mach fix of Thornber/Dellacherie/Rieper [67,15,53] is in a sense equivalent to the earlier 'enhanced' scheme, although appears to have been developed independently. Thornber and Rieper present numerical evidence that their schemes are accurate for 1D acoustics, although their analyses do not cover acoustic effects.
The enhanced/Low-Mach fix scheme is generally suitable for both convective and acoustic flows. However, Potsdam et al. 2007 [50] and Sachdev et al. 2012 [58] showed numerically that under certain conditions, it is susceptible to slight instabilities on the convective pressure field and acoustic velocity field (the instability of the acoustic velocity grid-mode was also demonstrated and analysed in [16]). By introducing timestep dependence into the diffusion formulation, they developed adaptive schemes which vary between different diffusion schemes depending on whether acoustic waves are resolved. This adaptive methodology was shown to maintain accuracy over a range of Strouhal numbers, and is applicable to Roe-type [50,6] and AUSM type [58] schemes.
Whereas the adaptive methodology aims to return to the preconditioned or unmodified scheme in situations where the enhanced/LMRoe scheme may produce oscillatory solutions, Bruel et al. 2019 [5] tried to modify the LMRoe scheme to eliminate the acoustic velocity instability and the related degradation in the CFL condition (σ < 1 2 ) of this scheme. Using discrete multiple-scale asymptotic analysis, they confirmed the suitability of the Roe, preconditioned Roe and LMRoe schemes for purely acoustic, purely convective, and either acoustic or convective flows respectively. By reintroducing off-diagonal diffusion terms into LMRoe the acoustic instabilities and CFL degradation were eliminated, although the final scheme produces symmetry-breaking solutions for some convective flows which the authors attribute to the scheme not maintaining Galilean invariance.

Design of artificial diffusion at low Mach number
In this section, we find the appropriate asymptotic scaling for the artificial diffusion terms of a numerical scheme, such that the scheme provides accurate results at low Mach number. We consider each of the three low Mach number regimes in turn. The diffusion scaling for purely convective flow is well-known, and its form was derived in Turkel et al 1994 andTurkel 1999 [71, 69]. We repeat this derivation here for completeness, before reapplying Turkel's method to derive the diffusion scaling for purely acoustic flow. By examining the two schemes, we can then design a diffusion scaling suitable for mixed convective-acoustic flow.
In equation (13) we augment the continuous x-split 2-dimensional Euler equations with artificial diffusion terms with a form that mirrors that of the flux Jacobian, which maximises the degrees of freedom without increasing the coupling between equations. The accuracy problem for the convective limit occurs only for multi-dimensional flow [15]. However, the vast majority of schemes in the class under consideration use some form of dimension-splitting (e.g. finite-difference derivatives along grid lines, or finite-volume interface fluxes over cell faces), and the dimension-split form is simpler than the full multi-dimensional form whilst still displaying the most important low-Mach behaviours, which justifies its use over a truly multi-dimensional analysis here. The dimension-split results can be readily extended to the y-split part of the 2D equations, to the 3D equations, and to the non-dimension-split equations. For analysis of the multi-dimensional modified equations see [15], and for a recent example of a truly multi-dimensional low Mach number scheme see [1].
We have written the artificial diffusion as linear second order terms for simplicity, but any consistent diffusive term is acceptable † . The results of the analysis still hold if, for example, we assume that x ψ ∼ O(M n ) for non-linear or fourth order diffusion terms respectively. This underlines the fact that we are interested in the limit of vanishing Mach number at a fixed, finite mesh spacing. If we also took the limit of vanishing mesh spacing, then the diffusion terms on the right hand side of the system (13) would disappear, as is necessary for a consistent numerical scheme.
The aim of this section is now to derive the appropriate asymptotic scaling for the coefficients A ij , such that the solutions to the modified system (13) have the same scaling behaviour as solutions of the original Euler equations at each low Mach number regime. Turkel's method for deriving these scalings for purely convective flow (which we will also use for purely acoustic flow) is as follows: 1. Find the limiting form of the left-hand side (L) of equations (13) by retaining only the largest term(s) as M → 0. 2. Force the artificial diffusion terms to be retained in the limit by choosing the order of the coefficients A ij such that all terms on the right-hand side (R) are the same order as the largest term(s) in L, i.e. R ∼ O(L) ‡ .
Because the mesh spacing is fixed, we can assume that ∂ x ψ ∼ O(M n ) =⇒ ∂ xx ψ ∼ O(M n ) and the scalings found in section 1 can be used for both L and R in step 1. These scalings are summarised in table 1. We can immediately find the artificial diffusion scaling for the linear vorticity and entropy fields (13c,13d). For these equations, L ∼ O(M 0 ) in both the single-and multiple-timescale limits, so for R ∼ O(L) as M → 0, we require: This is consistent with a standard upwind scheme with a diffusion coefficient on the order of the advecting velocity u, which is appropriate given that (13c,13d) have the form of the linear scalar advection equation. The pressure and velocity equations (13a,13b) form a coupled subsystem which is the distinguishing feature between the different low Mach number regimes and, by extension, between the different schemes. The rest of this section is concerned with finding the correct scaling for the artificial diffusion coefficients for this subsystem only, assuming the scaling (14) for the linear fields.

Artificial diffusion for purely convective flows
First, we find the required scaling of the coefficients A c ij , i, j = 1, 2 for purely convective flow. Using the scalings in table 1, the largest terms in L for (13a, 13b) are O(M 0 ). The diagonal coefficients A 11 and A 22 should be O(M −2 ) and O(M 0 ) respectively for the diagonal artificial diffusion terms to remain in the limit equations. By the same argument, the off-diagonal coefficients A 12 and A 21 should be O(M 0 ) and O(M −2 ) respectively § . Collecting these results, the convective scaling of the diffusion matrix is: 15) † A consistent diffusion term is proportional to ∆x n for some positive n, where ∆x is the grid spacing, so vanishes as ∆x → 0. ‡ Strictly, this should be R ∼ Θ(L), because if R ∼ o(L) then R ∼ O(L) still holds, but the diffusion R would vanish asymptotically. We will continue to use big-O instead of big-theta notation to remain consistent with the rest of the literature, with the understanding that a stricter interpretation is required for (at least the diagonal) diffusion coefficients.
§ It is acceptable for the off-diagonal terms to be o(L) and vanish in the limit, however the diagonal terms should be Θ(L) for stability.

Convective
Acoustic Table 1: Low Mach number scaling of the various terms in equations (13a, 13b) for purely convective or acoustic variation. Note that the time derivatives scaling for the acoustic variations of p and u make use of equation (9).
Retaining only the leading order O(M 0 ) terms from the modified equations (13a,13b), the limit equations for convective flow variations and convective diffusion scaling are † : The limit equations (16) are identical to the O(M 0 ) pressure and velocity relations (5a,4c) from the single-scale asymptotic expansion, with the addition of artificial diffusion terms. They also closely resemble the governing equations of Chorin's Artificial Compressibility Method [10], which was influential in much of the early work on low Mach number preconditioning [70,71,78]. As is appropriate for close-to-incompressible flow, the magnitude of both limit equations becomes independent of the Mach number, and at steady-state the velocity divergence will approach zero (up to the value of the diffusion terms) according to the pressure equation (16a).

Artificial diffusion for purely acoustic flows
Next we consider flow with purely acoustic variation, repeating the process above to obtain the correct scaling of the artificial diffusion A a in this regime. Using the acoustic flow variations from table 1, the largest terms in L for (13a,13b) are O(M 0 ) and O(M −1 ) respectively. For all artificial diffusion terms to be retained in the limit, we require: The limit equations using this diffusion scaling and acoustic flow variations are: which are the equations governing linear acoustics at low Mach number, as we found in section 1 with the relations (11b) and (10b), plus the artificial diffusion terms. Most schemes designed for transonic flow approach a diagonal approximation to this scheme, which resembles upwinding of the two acoustic waves (see appendix A). Note that if the off-diagonal terms vanish (A 12 ∼ o(M 0 ) and A 21 ∼ o(M −2 )), then the limit equations (18) are equivalent to Dellacherie's first order modified equations of a standard Godunov scheme for the linear acoustic equations (equations (61) and (62) in [15] with ν r = ν u k = a * ∆x/2M ), which are shown to destroy convective low Mach number accuracy. In practice, this scaling is not used for low Mach number simulations because it will create spurious acoustic waves from any non-trivial convective variations (as we will see later), and if the flow has zero (or enforced) convective variations then specialist acoustics solvers are much more suitable.

Artificial diffusion for mixed convective-acoustic flows
Acoustic variations are often the result of, or coexist with, convective phenomena, making it desirable to have a scheme with acceptable limit equations for both convective and acoustic variations. By examining the limit equations (16) and (18), we can identify two requirements for an acceptable form of the limit equations. The first is that all terms in L from the relevant physical relations are retained, i.e. the momentum and divergence relations (5a,4c) for convective flow, and the linear acoustics relations (11b,10b) for acoustic flow, which is necessary for the solutions to the modified equations to have the same scaling as solutions to the original equations. The second property is that some diffusion terms are retained in R, which will be necessary once we attempt to construct stable discrete schemes. Next we will see if the limit equations for convective diffusion with acoustic variations, and for acoustic diffusion with convective variations satisfy these requirements. The limit equations for convective diffusion scaling A c , equation (15), and acoustic flow variations are: resulting in a parabolic equation for the acoustic pressure p (1) . This means that the convective scheme A c effectively filters acoustic variations from the solution, just as the single-timescale expansion of the original Euler equations does, making this scheme unsuitable for purely acoustic or mixed convectiveacoustic flow. Conversely, the velocity equation (19b) retains all the required terms in L from the acoustic velocity relation (10b), although the diagonal diffusion term vanishes asymptotically. The limit equations for the acoustic diffusion scaling A a , equation (17), and convective flow variations are: resulting in a parabolic equation for the convective velocity u (0) . This will smooth out any convective variations of the velocity, making this scheme unsuitable for purely convective or mixed convective-acoustic flow. On the other hand, the pressure equation (20a) retains all the required terms in L from the convective pressure relation (5a), but the diagonal diffusion term is asymptotically vanishing. This over-damping of the velocity field and under-damping of the pressure field is a well-known problem with conventional (transonic) compressible schemes at the single-scale low Mach number limit [26,25,23].
We have shown that, unsurprisingly, neither the convective scheme A c nor the acoustic scheme A a completely fulfil the requirements for a mixed convective-acoustic scheme. However, the pressure limit equation (20a) with A a and the velocity limit equation (19b) with A c both retain the relevant terms in L, even if some diffusion terms vanish. This suggests we can form a scheme for mixed flow by combining the acoustic pressure diffusion A a 11 , A a 12 and the convective velocity diffusion A c 21 , A c 22 : which we call the mixed scheme. We shall see during the discrete analysis in section 5 that most modern schemes use this scaling, although to the authors' knowledge the only studies which explicitly identify the use of acoustic scaling on the pressure equation and convective scaling on the velocity equation are those of Venkateswaran and co-workers [73,75,50,58]. The limit equations for the mixed diffusion scaling with convective flow variations are: and with acoustic flow variations are: from which we see that for both convective and acoustic flow variations, all terms in L are retained compared to the desired form of the limit equations (16,18). The limit equations for the convective velocity and the acoustic pressure (22b,23a) retain both artificial diffusion terms, however the convective pressure and acoustic velocity limit equations (22a, 23b) retain only the off-diagonal diffusion terms. This may degrade the stability of the corresponding discrete scheme, especially if the off-diagonal diffusion terms vanish. The issues of degraded stability on the convective pressure and acoustic velocity were numerically demonstrated in [58], and the acoustic velocity instability was analysed in [16,5]; here we show that the cause of both can be identified in the continuous setting. Note that for vanishing offdiagonal terms, the limit equations (23) are equivalent to Dellacherie's first order modified equations for a low Mach Godunov scheme applied to the linear acoustic equations (equations (61) and (62), or (75), in [15] with ν r = a * ∆x/2M and ν u k = 0), which are shown to be accurate for convective low Mach number flow.

Adaptive schemes
We have identified three scalings of the diffusion matrix -convective, acoustic, and mixed -each suited to different flow regimes. However, it would be useful to have a single numerical method that could select the most appropriate scaling, either so the same method could be used for multiple problems, or because different flow regimes exist in a single problem [50]. The major difference between the schemes is whether they allow acoustic features, so a sensible metric for which scheme to use would be whether acoustic waves can be resolved. This can be estimated by the ratio of the acoustic timescale τ and the simulation timestep ∆t: When ∆t > τ acoustic waves cannot be temporally resolved and M u is small. When ∆t < τ acoustic waves can be resolved and M u is large. The parameter τ /∆t is normally presented as the unsteady Mach number M u . This parameter was first introduced by Venkateswaran & Merkle 1995 [72] to improve the conditioning of a low Mach preconditioned dual-time scheme, and has since been used to control the diffusion scaling [75,50,58,6]. In practice, M u is bounded between M and 1 so M u → 1 as ∆t → 0 and M u → M as ∆t → ∞. The coefficients A 11 and A 22 can be made to vary between the convective scaling when M u ≈ M and the acoustic scaling when M u ≈ 1 using: The effect of adaptive diffusion on accuracy and efficiency is discussed more in [50] and [58] but will not be covered in any more detail here.
In this section we have shown how artificial diffusion schemes can be derived for purely convective or purely acoustic low Mach number flow using Turkel's method [69] of balancing the artificial diffusion terms R with the physical terms L in the limit as M → 0. These schemes were then combined to create a mixed scheme suitable for both convective and acoustic flow, although this scheme has asymptotically vanishing artificial diffusion on the convective pressure and acoustic velocity equations.
We emphasise that these three diffusion scalings are not novel to the present study -all are in use today for low Mach number simulations, with the exception of the purely acoustic scheme, for the reasons discussed above. However, we have shown that many of the low-Mach behaviours of this class of schemes can be demonstrated independently of any specific discretisation.
In deriving and analysing these schemes so far, we have enforced the variations to obey either the convective or acoustic scalings from table 1. While this is useful for designing the diffusion schemes, we do not yet know how the schemes will perform if the variations are not enforced.
2. Transport equation relations. The convective momentum relation (4c), and the acoustic relations (10b,11b) should retain some artificial diffusion terms for the scheme to remain stable. 3. Divergence constraint relations. It can be argued that relations (5a,5b) should be matched exactly, so that the O(M 0 ) velocity field becomes divergence-free [33]. However, the inf-sup condition means that some discrete schemes are susceptible to pressure-velocity decoupling (chequer-board modes) in the incompressible limit. Certain collocated schemes [27,19] can avoid this issue by using Brezzi-Pitkäranta stabilisation [4], which introduces a pressure diffusion term into the continuity equation.
As noted in [23,15] this technique can be applied to density-based schemes, so it can also be argued that the velocity divergence should be zero only up to the value of some pressure diffusion term. Both choices are valid, and the impact of this choice will be discussed later.

Expansion of the Euler equations with convective diffusion
First, we carry out a single timescale expansion of (13a,13b) with convective diffusion scaling (15).
The relation (26c) is exactly the physical momentum relation (4c) plus all the artificial diffusion terms, as we required, and as found in the limit equation (16b Now we see that the spatial variations of p (0,1) are constrained by the two systems (26a,27a) and (26b,27b) respectively. Both systems admit the desired constant solutions, but we cannot yet tell whether they also admit non-constant solutions † . We will assume p (0,1) have constant initial conditions, and boundary conditions which are either periodic, or constant in both time and space. As for the physical relations (4), p (1) is assumed to have initial and boundary conditions equal to zero. These assumptions lead to constant solutions, satisfying the pressure variation relations (4a,4b). See Guillard & Viozat [26] for a more in-depth discussion of the conditions under which the discrete equivalents of these relations (which we will see later) lead to constant p (0,1) solutions. Note that it is useful to have restraints on both the first and second derivatives, despite the more involved reasoning about the constancy of p (0,1) . Assuming that central gradient approximations ) will be used in the final discrete schemes, the discrete ∂ x p = 0 constraint alone could lead to odd-even decoupling and allow non-constant grid modes. These modes are suppressed by the additional constraint ∂ xx p = 0.
Assuming constant p (0,1) then ∂ t p (0) → d t p (0) and relation (27c) is exactly the physical divergence relation (5a), plus all the artificial diffusion terms. The retention of the artificial pressure diffusion term allows for Brezzi-Pitkäranta stabilisation as discussed above. Therefore, all relevant relations from the single-timescale expansion are reproduced by the convective diffusion scheme.
The O(M −2 ), O(M −1 ) and O(M 0 ) pressure relations are: The p (0) relations (28a,29a) are identical to the single-timescale relations (26a,27a), so p (0) is constant under the same assumptions as before. If we assume no forcing of p (0) on the acoustic timescale -which seems reasonable -then ∂ τ p (0) = 0 and ∂ t p (0) → d t p (0) , and relation (29b) becomes a parabolic equation for p (1) equivalent to relation (19a). This will quickly damp any acoustic variations in the pressure, confirming that the convective diffusion scaling A c is unsuitable for simulations with acoustic variations.

Expansion of the Euler equations with acoustic diffusion
The single-timescale expansion of the modified equations (13a,13b) with A a leads to the following O(M −2 ) and O(M −1 ) velocity relations: and the O(M −1 ) and O(M 0 ) pressure relations: The p (0) relations (30a,31a) lead to constant p (0) under the assumptions described above. On the other hand, the relation (30b) means that ∂ x p (1) = 0 unless the velocity field is trivial and the term A 22 ∂ xx u (0) = 0, making the acoustic diffusion scheme unsuitable for convective low Mach number flows. For a diagonal scheme at steady-state, relations (30b,31b) resemble the Stokes flow equations, indicating that the solution will be dominated by the diffusive effects, balanced by variations in p (1) . The production of unphysical acoustic modes by A (−1) 22 , even from a well-prepared initial field, has been shown numerous times for specific discrete schemes [26,67,33,53], and in the continuous setting in [25,15].
A two timescale expansion of equations (13a,13b) with acoustic diffusion gives the O(M −2 ) and O(M −1 ) velocity relations: and the O(M −1 ) and O(M 0 ) pressure relations: The p (0) relations (32a,33a) lead to constant p (0) under the assumptions described above. The relations (32b,33b) are exactly the physical acoustic relations (10b,11b) plus all artificial diffusion terms, so we expect this scheme to be able to simulate at least purely acoustic flow. However, this scheme is still unsuitable for mixed convective-acoustic flow, as it will be unable to properly resolve the convective component as shown in the single timescale expansion.

Expansion of the Euler equations with mixed diffusion
The p (0) relations (34a,35a) lead to constant p (0) under the assumptions described above. On the other hand p (1) is constrained by relation (34b) and the continuity relation (35b). While the constant vector is clearly a solution to (34b) under the previous assumptions, it is now less clear whether it is the only solution. As we shall see later, almost all schemes which use the mixed diffusion scaling have asymptotically diagonal diffusion. In this case, (34b) becomes ∂ x p (1) = 0, so the physical relation (4b) is matched exactly, although as discussed earlier, the discrete scheme will be susceptible to odd-even decoupling on p (1) without the constraint on ∂ xx p (1) . At steady-state, the continuity relation (35b) for a diagonal scheme becomes (1) . This means that if p (1) is constant then u (0) is divergence-free, just as the incompressibility condition requires, but cannot be stabilised against chequerboard modes on p (≥2) in the manner of the Brezzi-Pitkäranta.
The last relation, velocity relation (34c), exactly matches (26c), which matches the physical relation (4c) with both diffusion terms retained. As every relevant physical relation is matched by the mixed diffusion scheme, we can expect that discrete schemes with matching modified equations will be suitable for simulations of purely convective low Mach number flow, although potentially at the risk of chequerboard modes on p (≥2) .
Two timescale expansion of (13a,13b) with mixed diffusion scaling (21) The p (0) relations (36a,37a) again lead to constant p (0) under the assumptions described above. Relations (36b,37b) match the physical acoustic relations (10b,11b), plus some artificial diffusion terms. The presence of ∂ xx p (1) in relation (37b) now provides the appropriate diffusion on the acoustic pressure variations.
On the other hand, the relation (36b) lacks a diagonal diffusion term, just as in the limit equation (23b), so a diagonal scheme will have no diffusion on the acoustic velocity variations. Relations (36b,37b) are almost equivalent to the modified equations of the scheme proposed by Bruel et al. (section 3.1.5.1 in [5]). This scheme differs from the current scheme in the precise form of the diffusion in the off-diagonal terms, but importantly the Mach number scaling of these terms matches ours. These differences will be discussed further in section 5. The mixed diffusion scheme A m is suitable for both convective and mixed convective-acoustic low Mach number flows. however, the scheme is not without disadvantages, having potentially degraded stability on the convective pressure and acoustic velocity variations.

Spectral radius estimates
We now estimate the scaling of the spectral radii Λ(A) of A c , A a and A m as M → 0. The spectral radius of the combined flux components determines the CFL bound for stability of explicit timestepping schemes and affects the convergence of implicit schemes. The flux Jacobian of the Euler equations has spectral radius u + a ∼ O(M −1 ), so the artificial diffusion will induce a more stringent stability bound if Λ(A) ∼ ω(M −1 ). Because diffusion matrices are positive (semi-)definite, we can estimate the scaling of the spectral radius with the scaling of the trace † . By inspection, the traces of the coefficient matrices (15), (17) and (21) tr(A a ) and tr(A m ) are the same order as the physical spectral radius, so should not change the stability requirements -up to a constant factor. On the other hand, tr(A c ) is larger than the physical spectral radius by an order of M −1 , so will require a CFL bound which decreases one order faster than the physical CFL bound in order to remain stable as M → 0. This scaling was first proven by Birken & Meister 2005 [2] for the specific case of preconditioned matrix-Rusanov artificial diffusion. The estimate (38) also agrees with [15] (equation (80) with κ r = M −1 ) that this restriction will hold for any scheme whose pressure diffusion converges to the convective limit as M → 0. For example, in the companion paper we will see † The trace cannot grow faster than the spectral radius for any fixed rank matrix. The trace of a matrix can grow more slowly than the spectral radius due to cancellation. For example, the trace of the flux Jacobian of the Euler equations in N dimensions is (N + 2)u, which is O(M 0 ), even though the spectral radius u + a is O(M −1 ). However, cancellation is impossible for a positive (semi-)definite matrix by definition. As the trace of a positive (semi-)definite matrix cannot grow asymptotically faster or slower than the spectral radius, it is an appropriate estimate for the scaling of the spectral radius.
Interface jump in normal velocity (u l − u i ) · n il Table 2: Nomenclature used for the discrete asymptotic expansions.
that this is true for existing AUSM schemes that converge to this limit, which explains the problems encountered in [45,29,8].
In this section, we have used single and multiple-scale asymptotic analysis to investigate each artificial diffusion scheme at the convective and mixed convective-acoustic limits. We have demonstrated several known properties of discrete low Mach number schemes. Including: the spurious forcing of the acoustic pressure p (1) by the acoustic scheme A a ; the over-damping of the acoustic pressure by the convective scheme A c ; the susceptibility to chequer-board modes on the convective pressure p (2) and instabilities on the acoustic velocity variations of the mixed scheme A m ; and the CFL limit of the convective scheme. So far we have worked in the entropy variables and the continuous setting to simplify the presentation, but in the next section we transfer the artificial diffusion to the conservative variables.

Analysis of the discrete Euler equations with artificial diffusion
The continuous analysis of the previous section is general to most schemes with the modified equations (13). In this section we demonstrate how the findings of the continuous analysis transfer to the discrete setting in the particular case of a first order cell-centred finite volume Roe-type scheme. First we identify the form of the interface flux in the conserved variables that matches the diffusion in the entropy variables (13), which we can then compare against previous Roe-type schemes from the literature. We then carry out the same six asymptotic expansions as in the previous section, but this time on the discrete equations. We finish the section by comparing the von Neumann symbols for each of the three schemes.

A general form for low Mach number finite-volume schemes
We consider the semi-discrete equations for a first order cell-centred finite-volume scheme in conservative variables q = (ρ, ρu, ρE) T , where the time derivatives are left continuous using a method of lines approach. Using the nomenclature of Guillard & Viozat [26] (table 2), the evolution in time of the solution at each cell i is described by the ODE: Where is the central approximation of the exact physical flux f (q) and f d il is the artificial diffusion flux between cells i and l. This nomenclature is general to any cell-centred scheme, but we will consider only quadrilaterals in 2D and hexahedra in 3D because, as mentioned in section 2, the accuracy problem for the convective limit disappears on simplex meshes [54,55,24,17,23]. In the non-dimensional entropy variables and a 2D face-aligned coordinate frame, we choose the elements A ij of the Jacobian of the diffusive flux to have the form: The first term is the convective upwinding. The second term is the diffusion on the pressure (first row) and normal velocity (second row) which was discussed in the previous sections. |v| is some O(M 0 ) velocity scale which is included to ensure correct dimensionality. µ α (α ∈ {u, 11, 12, 21, 22}) are positive expressions whose form and Mach number scaling are specific to each discrete scheme. The precise forms of the elements of this Jacobian are somewhat arbitrary, so long as their scaling can be chosen to match one of the three diffusion scalings A c,a,m -achieved here by varying the scaling of the coefficients µ α . This particular form is chosen because it simplifies the diffusion coefficients in the conserved variables. The off-diagonal terms may be negative so long as positive (semi-)definiteness of the Jacobian is maintained -this will be discussed in more detail below. Transforming (40) to the dimensional conserved variables, the diffusive flux between cells i and l is: where the interface velocity and pressure perturbations δU and δp are defined as: This is precisely the Liu & Vinokur form [42], which was also used by Weiss & Smith [78] and Li & Gu [33] for low Mach number Roe-type fluxes. The first term is the natural upwinding for the convective system. From the definitions of δU and δp we can see that the second and third terms are, respectively, the diffusion on the pressure and velocity equations in the entropy variables. We expand the coefficients µ α as: α is a constant positive real valued diffusion coefficient, for example = 1 gives the standard first order upwind diffusion, or 1/64 < < 1/32 can be used for fourth derivative diffusion [3]. Usually α are all equal for a particular scheme, although they do not have to be. The most common example of this is 11 for the convective scheme, which should instead be chosen in the recommended range for the Brezzi-Pitkäranta stabilisation coefficient. For the rest of this section we will use α = 1 for simplicity. ν α ∼ O(M 0 ) is some (non-dimensionalised) scheme specific expression and is the distinguishing feature between different schemes with the same diffusion scaling. The exponent n determines whether the diffusion has the convective, acoustic or mixed scaling, with the required exponents listed in table 3. The M −2 coefficients on µ 11 and µ 21 in (40) mean that the convective scaling exponents are all zero (i.e. M independent). Clearly, the convective upwinding coefficient should be µ u = u for all regimes. † Although expanding out µ α in this manner adds some complexity, it allows the separation of: the magnitude of the diffusion ( α ); the type of diffusion scaling (M n ); and the specific discrete form (ν α ) of each scheme. The full expressions for the diffusive fluxes can be found in appendix B.

Previous guidelines
We can compare the scalings in table 3 with guidelines put forward by previous authors. Dellacherie proved in [15] that for Godunov type upwind schemes, the velocity diffusion must be µ 22 ∼ O(M 0 ) for accuracy at the convective limit, which was also shown in [26,69]. It is clear that our findings for the convective and mixed scalings agree with this scaling. Dellacherie primarily studies the mixed scaling, but notes that the convective diffusion scaling (specifically the preconditioned Roe-Turkel scheme) also satisfies In their survey of low Mach number Roe-type schemes for the convective limit, Gu & Li [33] offered three guidelines for the design of low Mach number schemes. The first guideline states that µ 22 ∼ O(M 0 ) or smaller is necessary for accuracy at the convective limit, as previously discussed for [15]. The second guideline states that µ 11 should be between O(M ) and O(M 0 ) inclusive, where µ 11 ∼ O(M ) may allow some small pressure chequerboards, but µ 11 ∼ O(M 0 ) will suppress all pressure chequerboards. This † µu = uνu could be used where νu is some non-linear function used to control the dissipation level, such as in [29]. guideline matches the analysis for the convective and mixed scalings in the previous section, but again by including acoustic effects in the analysis we can see why this choice exists, and which flows each scaling is suitable for. Gu & Li also state that if µ 11 ∼ O(M 0 ) only in the continuity equation, very little improvement in the control of pressure chequerboards is seen compared to µ 11 ∼ O(M ). Equations (40) and (41,42) show that if µ 11 in δU is different in each equation then the equivalent diffusion in the entropy variables will not exactly match (40), which could explain the degraded performance. The third guideline states that a cut-off Mach number should only be used in denominators, i.e. when it decreases the diffusion. The issue of cut-off Mach numbers has not been covered here, but broadly speaking this guideline means, firstly, that µ 22 will not be increased beyond O(M 0 ) by a cut-off Mach number, which would compromise the accuracy for convective flow features, and secondly, µ 11 can still be decreased in the convective scaling, which alleviates the stability issue related to the O(M −2 ) spectral radius. See [33] for details of this implementation.

Classification of existing Roe-type schemes
Now that we have an expression for the artificial diffusion in the conserved variables (41,42), we can classify existing schemes as having either convective, mixed or acoustic diffusion scaling at low Mach number. Table 4 details many low-Mach Roe-type schemes from the literature. Some are adaptive schemes, and have different scalings for large timesteps using the convective CFL, σ u = u∆t/∆x ≈ 1, and for small timesteps using the acoustic CFL, σ a = a∆t/∆x ≈ 1. We have also noted whether each diffusion scheme is asymptotically diagonal (µ 12 , µ 21 ∼ o(M 0 )) or upper/lower triangular in the entropy variables. The discrete form of µ ij for a number of schemes in table 4 can be found in [33]. Several trends are apparent in this table. Almost all of the earlier schemes use the convective scaling [21,70,78,26,25,44]. These early schemes use a diffusive Jacobian derived from the preconditioned physical Jacobian (except [44]), which leads naturally to the convective diffusion scaling (see [69] and [26] for more detail). On the other hand, many of the more recent schemes use the mixed diffusion [73,75,50,67,65,15,53,58,49,5] even when the target regime is the convective limit. These schemes are derived primarily by identifying that for the original Roe scheme, which has the acoustic scaling, the accuracy problem is caused by A (−1) 22 . These schemes selectively reduce A 22 by a factor of M , thus reaching the mixed scheme. Most of these studies do not discuss the acoustic low Mach capability of the scheme -exceptions being Venkateswaran and coauthors [73,75,50,58] and Bruel et al. [5] † . All non-adaptive mixed schemes are diagonal except for the scheme of Bruel et al. [5].
The main development of adaptive schemes is due to Venkateswaran and co-authors [72,73,75,50,58], to allow accurate simulation of acoustic phenomena with the mixed scheme, but to return to the more favourable convergence and stability properties of the convective scaling when ∆t is too large to resolve the acoustic waves. To the authors' knowledge, [73] is the earliest use of the mixed diffusion scaling in the literature. The scalar diffusion schemes in [73,75,50] are diagonal for both large and small timesteps. The matrix diffusion schemes in [58] return to the preconditioned Roe scheme for large timesteps, where all four diffusion coefficients are well balanced, but for small timesteps they approach mixed schemes which are upper/lower triangular in the entropy variables (i.e. either µ 12 or µ 21 is o(M 0 )). ‡ The first scheme (equation (14) in [58]) has adaptive scaling on the δU diffusion terms µ 11 and µ 12 , meaning they both return to the values in the original Roe scheme for small timestep -µ 11 becomes O(M ) and µ 12 becomes o(M 0 ), so A is lower triangular. Conversely, the second scheme (equation (15) in [58]) has adaptive scaling on the pressure diffusion terms µ 11 and µ 21 , so becomes upper triangular. The authors state that the matrix diffusion scheme of Potsdam et al. [50] is almost equivalent to the second scheme in [58]. Despite these differences, Sachdev et al. report no visible differences between results with these three schemes. † Dellacherie [15] analyses the acoustic equations in depth, but makes only a minor distinction between the convective and mixed diffusion scalings, instead focusing on the necessity of reducing µ22 compared to the acoustic scaling. ‡ This can be seen by identifying that Cφ T and C p in equations (12,13) in [58] correspond to the δU and δp terms respectively in (41,42).
Interestingly, the unsteady momentum interpolation method of Li & Gu [34] also appears to be an adaptive scheme. The earlier All-speed scheme [32,35] approaches the convective scheme of Mary & Sagaut [44] at low Mach number. However, in [34] the pressure diffusion µ 11 is scaled by the timestep so will be a factor of M smaller with σ a ≈ 1 than with σ u ≈ 1, hence approaching the mixed scaling for small timesteps. This can be seen from equations (19,20) and figures 3-4 of [34] where slight chequerboards are seen for small timesteps but are controlled for larger timesteps.
The final point to note from table 4 is that a few schemes reduce at least some components of the convective upwinding term by a factor of M [65,15,49]. We cannot assess the effects of this change using our analysis so far, but we will return to consider it at the end of the section when we examine the von Neumann symbols of the first order scheme (39).

Behaviour of the off-diagonal diffusion
The off-diagonal diffusion terms in (42) contain the switch U il /|U il | = ±1. This ensures that these terms have a diffusive character but also gives the possibility of negative diffusion coefficients. † Previous schemes often use just U il , but we prefer a non-dimensional O(1) expression. In 2D the off-diagonal terms result in anisotropic diffusion of the form (±∂ xx u ± ∂ yy v) on all equations through δU , and ±∂ xx p on the xmomentum equation and ±∂ yy p on the y-momentum equation through δp. The sign on ∂ xx u and ∂ xx p will be the same, as will the sign on ∂ yy v and ∂ yy p, but the sign may differ between the x and y derivatives. As stated earlier, off-diagonal anti-diffusion is tolerable so long as positive (semi-)definiteness is maintained. This requires µ 11 µ 22 − µ 12 µ 21 ≥ 0, where equality implies semi-definiteness. For asymptotically diagonal or triangular schemes this holds trivially. This covers all acoustic schemes, almost all mixed schemes and some convective schemes in table 4. The full convective schemes are due to the preconditioned diffusion, and the inequality holds for these schemes. ‡ For the mixed scheme µ 11 µ 22 vanishes asymptotically for both convective and acoustic variations, therefore µ 12 µ 21 must also vanish asymptotically else the diffusion will be negative definite and the scheme unstable. § The one exception to this is the full scheme of Bruel et al [5], where the off-diagonal terms are designed to remove the acoustic instability and to obtain a CFL limit σ a < 1 compared to σ a < 0.5 for the diagonal mixed scheme. This scheme modifies the off-diagonal diffusion to be isotropic with the form ∇ 2 u + ∇ 2 v in δU and ∇ 2 p on all momentum equations through δp. Positive semi-definiteness is ensured by enforcing that the off-diagonal terms have opposite sign, but which term is positive and which negative is indeterminate. However, the scheme produces asymmetric solutions for flow around a circular cylinder, which Bruel et al attribute to this free choice of sign and the resulting loss of Galilean invariance.

Expansion of the discrete Euler equations with artificial diffusion
In this section we carry out single-and multiple-scale asymptotic expansions of the discrete equations (39) with diffusive flux (41,42) for each of the three scalings in table (3). For simplicity we have taken α = 1 throughout. We will keep the discussion of these expansions brief, because most points have already been discussed in the previous section, and many of the expansions have been presented in the literature (although usually for one specific flux scheme). The discrete single-timescale expansion with the convective diffusion is investigated in detail in the seminal papers from Guillard & coauthors [26,25], and the most important features are also reported in the reviews by Guillard & Nkonga [23] and Li & Gu [33]. The discrete single-timescale expansion with the acoustic diffusion can be found in a several papers, as it shows why classical transonic schemes fail at low Mach number, see [26,25,23,53]. Several of these papers also present the discrete single-timescale expansion with the mixed diffusion [15,53,32,36,33,23]. To the authors' knowledge, the only previous study in the literature to carry out multiple-timescale expansions of the discrete equations is Bruel et al. [5], who present multiple-timescale expansions of the baratropic Euler equations with all three schemes.

Expansion with convective diffusion
First, we carry out the single timescale expansion of the discrete equations (39) with convective diffusion scaling. The relations are exactly equivalent to those found by Guillard & Viozat [26], except for the more general notation for the diffusion coefficients, so will only be discussed briefly -consult [26] for more in-depth discussion on these discrete equations, they are included here mainly for completeness. The † The authors would like to thank the anonymous reviewer who demonstrated that these terms did not have a diffusive form in an earlier draft. Correcting this led to the discussion in section 5.2.3.
The right-hand sides of (44a,44c) are in the form of the discrete Laplacian. In N dimensions, equations (44) are two elliptic systems each of N + 2 equations for p (0,1) , which are the discrete equivalents of the continuous relations (26a,27a) and (26b,27b). Under a small number of reasonable assumptions, they enforce constant p (0,1) over the entire domain [26]. Using the non-dimensional O(M 0 ) equation of state [48], the O(M 0 ) relations are: Where we have used constant p (0) to imply constant ρE (0) and ρH (0) to simplify relation (45c). If the zeroth order entropy s (0) is constant, then the zeroth order temperature and density T (0) and ρ (0) also become constant. The second term in (45a) is the discrete divergence, and relations (45a,45c) are the discrete equivalents of the continuity relation (27c). At steady state on a regular grid, the pressure diffusion takes exactly the form of the Brezzi-Pitkäranta stabilisation [4] used with a finite-volume scheme by Eymard et al. [19], so this scheme is expected to be free of chequerboard instabilities. The momentum relation (45b) is the discrete equivalent of relation (26c) † , having all the necessary terms in L, and properly balanced diffusion in R. The discrete single scale expansion with the convective diffusion scaling therefore reproduces all the necessary relations for the convective low Mach number limit, so is a consistent scheme for this limit. Next, we carry out the two timescale expansion with the convective diffusion scaling. The O(M −2 ) relations are identical to those from the single-timescale expansion, so we assume p (0) and related zeroth order thermodynamic quantities are constant again. The O(M −1 ) relations are: Using the O(M 0 ) equation of state, we can see that relations (46a,46c) are the discrete equivalents to (29b), which is is a parabolic equation for p (1) , assuming ∂ τ p (0) → 0. This means p (1) will become constant on the acoustic timescale τ unless p (0) is forced on the acoustic timescale (which would violate the physical relation (11a)), so the convective diffusion scaling is indeed unsuitable for acoustic or mixed convective-acoustic simulations.

Expansion with acoustic diffusion
The single timescale expansion of the discrete equations (39) with acoustic diffusion scaling leads to the O(M −2 ) momentum relation: and the O(M −1 ) relations: The relations (47,48a,48c) are the same elliptic system for p (0) as we have seen previously, so we assume p (0) is constant. We have used a constant p (0) to simplify relation (48b), which is a discretisation of the continuous relation (30b). The first order pressure p (1) cannot be constant unless the jumps in normal velocity ∆ il U (0) are zero at every cell face -this is even clearer if the scheme is diagonal, which the acoustic scheme usually is. As previously stated, this condition is actually enforced for simplex meshes, leading to a divergence free solution with the correct pressure scaling [54,55,24,17]. This is why we restricted our analysis to quadrilaterals and hexahedra at the beginning of this section. For these cell types, (48b) forces either trivial velocity fields or non-constant p (1) , so this scheme is unsuitable for convective low Mach number flows. The two timescale expansion with acoustic diffusion scaling leads to the same O(M −2 ) momentum relation (47), and the following O(M −1 ) relations: The O(M 0 ) density relation is: Relations (47,49a,49c) are now an unsteady elliptic system on the acoustic timescale for p (0) , leading variations in p (0) to be dissipated on the acoustic timescale. Constant p (0) has been used to simplify relations (49b,50), which for constant entropy (i.e. constant ρ (0) ) are the discrete equivalents of the linear acoustics relations (32b,33b). This shows that the scheme is suitable for simulating purely acoustic low Mach number flows (for example 1D shocktube flows), but not mixed convective-acoustic flows due to the issues highlighted in the single-scale expansion.

Expansion with mixed diffusion
Lastly, we carry out the expansions with the mixed diffusion scaling. For much of the discussion we assume that the scheme is diagonal for simplicity, but we have retained both off-diagonal diffusion terms in the expansions on the understanding that at least one of µ 12 and µ 21 must vanish according to the discussion in section 5.2.3. The single timescale expansion leads to the O(M −2 ) momentum relation: and the O(M −1 ) relations: Relations (51,52a,52c) are the elliptic system for p (0) . Assuming constant p (0) , (52b) becomes an elliptic equation for p (1) , for which constant p (1) may not be the only solution. For diagonal or upper triangular schemes the right hand side of relation (52b) vanishes. On a regular 2D Cartesian grid, the momentum relations then become p i,j−1 = 0, which will, on their own, admit chequerboard modes. The O(M 0 ) relations are: Relations (53a,53c) are discrete equivalents of the continuity relation (35b). Assuming µ 12 ∼ o(M 0 ), at steady state these relations mean that p (1) is free of chequer-boards iff the discrete divergence of the zeroth order velocity is zero. In practice, we know of no case -either in the literature or in our own studies -where this scheme has produced steady solutions with non-constant p (1) , and Rieper [53] provides an argument that this scheme will damp the divergence on the convective timescale on Cartesian grids. Although approaching a divergence free solution as M → 0 is appealing, it means that this scheme will be susceptible to chequerboard-modes on p (2) , as we shall see in numerical examples later. Lastly, the momentum relation (53b) is consistent with the continuous relation (34c). The scheme therefore reproduces all the relevant asymptotic relations, and is suitable for simulating convective low Mach number flow, but with the caveat of potential chequer-board modes on p (2) . The two timescale expansion with mixed diffusion scaling gives the same O(M −2 ) momentum relation and O(M −1 ) density and energy relations as for the two timescale expansion with acoustic diffusion (47,49a,49c), which lead to constant p (0) . Subsequently, the O(M −1 ) momentum relation is: and the O(M 0 ) relations are: The relations (54) and (55a,55c) are the discrete equivalents of the linear acoustics relations (36b,37b). It can be seen that, as in the continuous analysis, the velocity diffusion vanishes from the acoustic momentum relation (54), while the acoustic relations (55a,55c) retain the pressure diffusion term. Therefore this scheme is consistent for acoustic flows, however potentially has instabilities in the acoustic velocity variations. The momentum relation (55b) is the discrete equivalent of (36c), confirming that this scheme is suitable for mixed convective-acoustic flows.
The discrete expansions of the first order finite-volume scheme have confirmed that all findings of the continuous analysis in sections 3 and 4 carry over to the discrete setting for this particular discrete form, with each diffusion scaling behaving as expected for both the single-and multiple-scale limits.

Stability of the discrete scheme
We now consider the stability of the discrete schemes by finding expressions for the eigenvalues of the artificial diffusion Jacobian and the von Neumann symbols of the fully discrete first order scheme.

von Neumann symbols of the first order scheme
We now find the stability limits for the linearised first order scheme with each diffusion scaling using von Neumann analysis. Discretising equation (39) using central differences in space and forward Euler in time, we have the linearised update for the solution vector q n i in cell i at time level n: where J is the physical flux Jacobian and A is the diffusive flux Jacobian. In this section we assume all diffusion matrices are diagonal in the entropy variables, as it simplifies the presentation. For stability analysis with the off-diagonal terms included, see [2] for the convective scheme and [5] for the mixed scheme. Substituting q n j =q n e ikj∆x , where k = nπ∆x/l is the discrete wavenumber, we find the amplification matrix in Fourier space is:Ĝ where I is the identity matrix, S k/2 = sin(k/2) and S k = sin(k). For stability, the eigenvalues λ i of (Ĝ(k)) must be |λ i | ≤ 1 ∀ k. We consider the 1D equations, as stability of the 1D scheme is a necessary condition for stability of the multidimensional schemes. In the entropy variables the entropy equation decouples and the associated eigenvalue is identical for all three diffusion schemes, and is exactly that of the first order upwind scheme for scalar advection: where σ a = a∆t/∆x is the acoustic CFL number. This eigenvalue depends on the convective CFL number σ a M , as we would expect, and is stable for 0 ≤ (σ a M ) ≤ 1. The eigenvalues for the pressure-velocity subsystem using the acoustic diffusion scheme are: Where we can identify the upwind discretisations of the convective system in the second term, and the forward/backward travelling acoustic waves in the third term. The stability is dominated by the acoustic terms so, to leading order, the scheme is stable for 0 ≤ σ a ≤ 1. The full expression for the eigenvalues for the convective scheme is: The second term is the convective subsystem, the third term is due to the large (σ a M −1 ) pressure diffusion, and the radicand has contributions from both the pressure diffusion and from the acoustic waves. The leading order term is the (σ a M −1 ) pressure diffusion, so the stability limit will scale a factor of M worse than the usual stability limit (which depends on σ a ), as originally shown for preconditioned diffusion by Birken & Meister [2] and as predicted for all convective schemes by the spectral radius estimates (38) in section 4. Contours of the real and imaginary parts of λ c 1 and λ c 2 for varying k and M are plotted in figure 1. The eigenvalues have two distinct regions, either side of the locus M = tan(k/2)/2 where the radicand goes to zero (labelled R = 0 in figure 1). In the low wavenumber region to the left of this locus the radicand is positive so this term is imaginary, and to the right in the moderate/high wavenumber region the radicand is negative and this term is real. We can simplify the square root in each region using binomial expansions for k << M and k >> M . In the low wavenumber region this leads to the following expression: and in the high wavenumber region: These expressions can be interpreted from the perspective of a multiple-space scale approach, where acoustic and convective phenomena occur on large and small spatial scales respectively, instead of the multiple-time scale approach used in the rest of the paper. In the low wavenumber 'acoustic' region, the leading order imaginary term for both eigenvalues is due to the σ a forward/backward travelling acoustic waves. However, these terms are dominated by the (σ a M −1 ) diffusion, so to leading order both eigenvalues resemble those of a diffusion equation, corresponding to the rapid damping of acoustic waves by this scheme. Long wavelength variations in both pressure and (normal) velocity are due to acoustic effects and must be damped.
On the other hand, in the high wavenumber 'convective' region the leading order imaginary terms in both eigenvalues are on the (σ a M ) convective scale, and only λ c 2 is dominated by the (σ a M −1 ) diffusion. At short wavelengths, (normal) velocity is no longer an acoustic quantity, but joins entropy (and vorticity in 2/3D) as a convected quantity. The pressure is still dominated by the large diffusion, which means that the Brezzi-Pitkäranta type stabilisation is retained across all wavenumbers.
The boundary between the two regions reduces approximately linearly with M because of the scale separation between acoustic and convective phenomena as M → 0. In real applications implicit timestepping is used for this scheme to overcome the stringent (σ a M −1 ) stability limit, but the von Neumann analysis of the explicit first order scheme provides an interesting insight into the behaviour of the scheme.
The first two eigenvalues of the mixed scheme are: We can identify the convective second term proportional to (σ a M ). The acoustic third term is proportional to (σ a ) but the diffusive term S 2 k/2 is half that of the standard upwind scheme, and the advective term has an additional contribution from the second term in the radicand. The magnitudes and imaginary component of the eigenvalues λ m 1,2 , ignoring the convective term, are plotted in figure 2 alongside the equivalent values for the standard upwind scheme. For the low wavenumber range 0 ≤ k < 2 arctan(2) ≈ 2.2, the radicand is real and the scheme behaves similarly to the upwind scheme except with lower diffusion, evident in figure 2a, and an additional dispersion error, evident in figure 2b. For the high wavenumber range 2 arctan(2) ≤ k ≤ π, the radicand is imaginary and the scheme is entirely diffusive. Figure 2a shows that |λ i | bifurcates in the high wavenumber range. As π − k = k → 0, the eigenvalues λ m 1,2 can be expressed as: where C k/2 = cos(k/2). λ m 1 approaches the pure diffusion scheme as k → π, as does the upwind scheme. To leading order, λ m 2 approaches 1 independently of σ a , indicating that there is a grid mode which is undamped by the mixed scheme. This was shown by Dellacherie [16] by considering the time evolution of the energy of the grid mode using the mixed scheme. The low wavenumber range is stable under the reduced CFL condition 0 < σ a < 0.5, which is the same as found by Dellacherie [15], and Bruel et al. [5]. The high wavenumber range is (linearly) stable under the CFL condition 0 < σ a < 1, although nonlinear interactions in the full scheme will excite the undamped grid mode in many circumstances.
Lastly, we return to the question of reducing the convective upwinding by a factor of M, as done by Thornber [65] and Oßwald et al. [49] in the velocity equations. The eigenvalue of the upwind scheme for the convective terms (61) using a diffusion coefficient of α|U |/2 instead of |U |/2 is: which is linearly stable under the CFL condition 0 ≤ (σ a M ) ≤ α. If α = M and σ a ≈ 1, then the upwind scheme with reduced diffusion is linearly stable. This argument also justifies reducing the convective upwinding for every component, not just the velocity. However, because we rely on σ a ≈ 1 this argument does not hold for preconditioned schemes where σ u ≈ 1 and σ a ≈ 1/M .

Numerical examples
We now show numerical examples which demonstrate the behaviours of the three low Mach number schemes. The discrete equations for a cell-centred finite volume scheme (39) are solved using a first order scheme in both time and space. The first order explicit Euler scheme is used for time integration. Diagonal diffusion is used for all examples. Following the discussion of section 5.2.3, the mixed scheme with full diffusion is unstable with the diffusion form (41,42), and upper/lower triangular mixed diffusion was found by Potsdam et al [50] and Sachdev et al [58] to have the same behaviour as the diagonal scheme. For the acoustic scheme diagonal diffusion is equivalent to standard upwinding (see appendix A). For the convective scheme almost no difference was seen between full or diagonal diffusion, so the diagonal diffusion results are shown for consistency with the other schemes. If higher order reconstructions and time integration are used, higher resolution results are obtained, but the qualitative behaviour remains the same, although smaller differences between the fluxes are observed due to the reduction in the diffusion. All simulations are carried out with an ideal gas with ratio of specific heats γ = 1.4 and specific gas constant R = 287.058.

One dimensional examples
In one dimension, the only solutions which are compatible with the convective limit (4,5) are trivial, with both constant pressure and velocity [15]. As such, all flow variations are either acoustic or entropy waves, and one dimensional examples can be used to isolate the behaviour of the numerical fluxes for purely acoustic flow.

Isolated soundwave
The first test case is an isolated soundwave in one dimension. This test case will show the differences between the fluxes for a smooth purely acoustic flow. The pressure p(x, t) is initialised with a sinusoidal profile, the density ρ(x, t) is initialised assuming isentropic flow, and the velocity u(x, t) is initialised using the Riemann invariant for a forward travelling sound wave:   (62) and (66).
The acoustic scheme has diffused the pressure peaks roughly twice as much as the mixed scheme. The convective scheme however has almost entirely smoothed out the wave, which has travelled only a small portion of a wavelength. These results are entirely consistent with our earlier analysis which found that the acoustic and mixed schemes are both suitable for smooth acoustic flow, whereas the convective scheme is overly diffusive of acoustic variations.

Low Mach number shocktube
The second one dimensional case is a low Mach number shock tube. This case demonstrates the performance of the different fluxes for discontinuous purely acoustic flow, and was first used by Sachdev et al. [58] for testing adaptive schemes. The left and right initial states are: The very small pressure difference produces a contact wave moving to the right at M ∞ = 0.0001 between two receding weak shockwaves. Figure 3b shows a close-up of the Mach number distributions between the two shock waves obtained with the acoustic and the mixed flux schemes after 96 timesteps with a CFL of σ = 0.4. The acoustic flux gives a monotone solution, as expected for a first order, upwind scheme. However, the solution found with the mixed flux has significant oscillations originating at the shockwaves -an undamped grid mode as predicted by the von Neumann analysis.
From the one-dimensional examples we can verify that: the convective scheme is completely unsuitable for flow with acoustic variations; the acoustic scheme is suitable for flows with both smooth and discontinuous acoustic variations; and the mixed scheme is suitable for flows with smooth acoustic variations, but has too little diffusion to properly handle acoustic variations close to the grid scale.

Two dimensional examples
Non-trivial solutions to the convective limit exist in two and higher dimensions. We present two numerical examples in 2D, one which tests the schemes' capabilities for steady purely convective flow, and one which tests the capabilities for unsteady mixed convective-acoustic flow.

Circular cylinder
This classic test case will demonstrate the performance of each flux scheme for steady purely convective flow. The correct solutions should tend to the incompressible solution as M → 0. The farfield state is: The background pressure is calculated from p ∞ = ρ ∞ /(γM 2 ∞ ). The cylinder is centred at the origin with radius r = 1, and is meshed using an O-mesh which extends out to 50 radii from the origin, with 64 and 48 cells in the circumferential and radial directions. The first row of cells next to the cylinder has a radial height of 0.036r. The cell height grows at a rate of 1.11, with the final row having a radial height of 4.9r. Curvature-corrected boundary conditions are used at the inviscid cylinder wall [12], and ghost cells at freestream boundaries are set using the upstream velocity and entropy and the downstream static pressure. Convergence is reached by pseudo-timestepping at a CFL of σ = 0.4, which is accelerated by local timestepping and preconditioning using Weiss & Smith's preconditioning matrix [78] (preconditioning is only used for the convective and mixed schemes, as preconditioning does not reduce the spectral radius of the acoustic scheme [21]).
The results for the convective flux is shown in figure 4b, and closely resemble the exact pressure distribution from the classical incompressible potential solution in figure 4a. The solution found with the acoustic flux is shown in figure 4c, and is visibly very different from the incompressible solution. The artificial diffusion dominates over the other terms, and the flow more closely resembles a Stokes flow than the inviscid solution, with significantly larger pressure variations. Lastly, the solution found with the mixed flux is shown in figure 4d and is very similar to the convective scheme and the incompressible solution. However, on close inspection there is a small chequer-board mode in the radial direction. This is more easily seen in figure 5 which shows the difference between the exact pressure and the pressure found with the mixed and convective flux schemes. The error for the convective flux is smooth with no chequer-board modes, however the error for the mixed flux clearly shows chequer-board modes in the radial direction.

Soundwave-vortex interaction
The last numerical example is the interaction of a one-dimensional soundwave passing through a low Mach number Gresho vortex [22,41] with stationary background conditions. This provides a very simple test of a numerical scheme's ability to not only simulate both convective and acoustic features, but also their interaction, and has been used previously by [5,47].
The domain is (x, y) ∈ [−2h, 2h) × [−h, h) with periodic boundary conditions in both dimensions, and is discretised with 256 × 128 square cells. The initial conditions are shown in figure 6a, composed of the superposition of a Gresho vortex centred at (0, 0) with a diameter 0.4h, and a right-travelling soundwave with a Gaussian profile centred along the line (−h, y). The soundwave profile is calculated according to equations (69), with: Where G(x) is the Gaussian function with standard deviation of 0.1h, c 1 is chosen such that the maximum acoustic velocity is 1, and p ∞ = ρ ∞ /(γM 2 ∞ ) again. The Gresho vortex has circumferential velocity u θ (r) and pressure p(r) according to [47]: The results for each flux after a single acoustic period (4h/a ∞ ) and after ten acoustic periods (equivalent to 1/10th of a convective period 4h/(a ∞ M ∞ )) are shown in figures 6 and 7 respectively, using a CFL of σ = 0.4.  Figure 6b shows the solution found with the convective scheme after one acoustic time. The vortex shape is very well preserved, with barely any reduction in peak velocity. The soundwave on the other hand has been completely dissipated and is no longer visible. The only indication of the soundwave is the slight asymmetry in the vortex velocity over the y axis, which happens as the soundwave is smeared out over the domain.
The solution from the acoustic scheme is shown in figure 6c. The vortex has a drastically reduced peak velocity, and is misshapen, having become aligned with the grid due to the anisotropic artificial diffusion on the face-normal velocity. The soundwave is still visible, having travelled once around the domain. It is also damped, although not to an unsurprising degree for a first order upwind scheme, given there were fewer than 20 cells across the initial width of the soundwave.
As predicted by the asymptotic analysis, the results for the mixed scheme (figure 6d) combine the favourable behaviour of both the convective and acoustic schemes. Both the vortex and the soundwave are well preserved. The vortex shows comparable diffusion to the convective scheme (although without the asymmetry), and the soundwave is less damped than that found with the acoustic scheme, in agreement with the von Neumann analysis in the previous section.
The results for the convective scheme after 10 acoustic periods (0.1 convective periods) are shown in 7b. The vortex shape is mostly preserved with only minimal distortion around the corners, where the flow is most misaligned with the grid. The peak vortex velocity is around 75% of the initial peak value. The acoustic scheme results are shown in figure 7c. The vortex is even more distorted than in figure 6c, and has polluted a large part of the domain. After integration over 10 periods, the soundwave is barely visible on the left side of the figure behind this distortion. The result for the mixed scheme again retains the favourable properties of the other two schemes (figure 7d). The velocity magnitude of the vortex is similar to that found with the convective scheme, but with marginally less distortion to the shape. The soundwave is almost completely diffused, although less so than with the acoustic scheme.
We now investigate the effect of reduced convective upwinding, as proposed in [65,53,49]. Figure 8b shows the solution after 10 acoustic times found with the mixed scheme where the convective upwinding -the first term in (40) -has been reduced by a factor of M . The peak vortex velocity is 87% of the initial value compared to 75% with the original mixed scheme (figure 8a), a moderate improvement. Inspection of the diffusion matrix (40) shows that for the mixed scheme the µ 22 term is of the same magnitude as the convective upwind term. Figure 8c shows the solution found by reducing only the µ 22 term by a factor of M . The peak velocity is 82% of the initial value, a comparable improvement to reducing the convective   upwind term. Figure 8d shows the results found by reducing both the convective upwinding and µ 22 terms by a factor of M compared to the original mixed scheme. The improvement is striking, with a peak vortex velocity of 99% of the initial value. Some small instabilities can be seen around the corners of the vortex, where the flow is most misaligned with the grid and the 1D von Neumann analysis of section 5 is least valid. Although figure 8d shows a significant reduction in the diffusion without significant loss of stability, recall from the von Neumann analysis that the stability of this scheme relies on σ a ∼ O(1), making it unsuitable for preconditioned schemes where σ u ∼ O(1).
The convective upwinding reductions in [53,65,49] do not exactly correspond to those discussed here. Rieper [53] reduces the jump in the face-normal velocity component by a factor of M , corresponding to reducing the µ 22 term from the acoustic scheme scaling to the mixed scheme scaling, and reducing by a factor of M the convective upwinding on the face-normal component of velocity only. Thornber et al. [65] and Oßwald et al. [49] reduce the jumps in all velocity components by a factor of M , so for vortical flows are equivalent to the reduced convective upwinding scheme in figure 8b. However, reducing the entire convective upwinding term also reduces the diffusion on the entropy convection, so can be expected to give better results for cases with significant entropy variations such as heat transfer simulations.
In this section we have used numerical examples to verify the findings of sections 3, 4 and 5. The one dimensional examples demonstrated that the acoustic scheme is well suited to purely acoustic flow, while the convective scheme is unable to resolve any acoustic phenomena. The mixed scheme is suitable for smooth acoustic flow -resolving these phenomena with less diffusion than the acoustic scheme -although is unstable for discontinuous acoustic flow, possessing an undamped grid mode. Using steady flow around a two dimensional cylinder, we verified that the acoustic scheme is unsuitable for purely convective flow, whilst the convective scheme approaches the incompressible limit as M → 0. The mixed scheme also approaches the incompressible limit, avoiding the catastrophic failure of the acoustic scheme, however we saw it is susceptible to chequerboard modes on the second order convective pressure. Finally we saw that neither the acoustic nor convective scheme is suitable for mixed convective-acoustic flow. The mixed scheme on the other hand can resolve both convective and acoustic features in the same flow, and can also achieve this with dramatically reduced convective upwind diffusion.

Conclusions
Low Mach number flows include a range of applications of scientific and engineering interest. Collocated, density based solvers for compressible flow are one method of simulating low Mach number flow, but care must be taken to ensure that the artificial diffusion does not compromise their accuracy in this regime. In this paper we have reviewed the behaviour of the artificial diffusion in this class of schemes at low Mach number using the modified equations. By considering both the convective and acoustic low Mach number limits, we have shown how three artificial diffusion scalings naturally arise in the entropy variables -one suitable for purely convective flow, one for purely acoustic flow, and one for flow with both convective and acoustic features. Single-and multiple-scale asymptotic expansions of the modified equations established the behaviours of each scheme for different flows. The convective scaling is compatible with convective flows, but damps out acoustic variations on a very fast timescale. The acoustic scaling is compatible with acoustic flows, but will produce spurious pressure waves for convective flow. The mixed scaling is compatible with both convective, acoustic, and mixed flow, but has vanishing diffusion on the convective pressure and the acoustic velocity, which may lead to pressure chequerboard modes and acoustic grid-mode instabilities respectively.
Transforming the artificial diffusion to a Roe-type finite volume scheme in the conserved variables enabled us to compare to a number of existing low Mach number methods. Each of these methods matched one of the three scalings, and our analysis agrees with previous theoretical analyses and wellknown behaviours of these existing methods. The convective and mixed scalings conform with previous guidelines for accurate schemes for convective low Mach number flow [15,33], but by considering accuracy requirements for both convective and acoustic effects, we were able to explain why there are two possible scalings suitable for the convective limit, and what their relative advantages and disadvantages are. The price of the mixed scheme's flexibility is its compromised stability and the additional constraints on the off-diagonal diffusion. It is the authors' belief that remedying the stability of this scheme, and in particular whether the off-diagonal terms can be leveraged to achieve this, is the most pertinent open question in this research area. Asymptotic expansion of the discrete equations showed that all findings of the continuous analysis apply for the first order Roe-type finite volume scheme, and von Neumann analysis confirmed the stability estimates obtained from the continuous analysis. Finally, four numerical examples demonstrated the performance of each diffusion scaling for acoustic, convective and mixed convective-acoustic flow.
There is a significant body of literature investigating this class of schemes in the limit of vanishing Mach number. We have shown that the most important behaviours of this class of schemes at low Mach number can be found and explained in a simple manner using the continuous modified equations in the entropy variables. This form can be used to compare schemes and predict their capabilities independently of the specific discretisation, as well as to provide guidelines for the development of novel low-Mach schemes.

Competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.