A quantitative comparison of high-order
asymptotic-preserving and asymptotically-accurate
IMEX methods for the Euler equations with non-ideal
gases

Giuseppe Orlando, Sebastiano Boscarino, Giovanni Russo

► To cite this version:

Giuseppe Orlando, Sebastiano Boscarino, Giovanni Russo. A quantitative comparison of high-order asymptotic-preserving and asymptotically-accurate IMEX methods for the Euler equations with non-ideal gases. Computer Methods in Applied Mechanics and Engineering, 2025, 442, pp.118037. 10.1016/j.cma.2025.118037. hal-04907464v4

HAL Id: hal-04907464

https://hal.science/hal-04907464v4

Submitted on 22 Oct 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A quantitative comparison of high-order asymptotic-preserving and asymptotically-accurate IMEX methods for the Euler equations with non-ideal gases

Giuseppe Orlando$^{(1)}$, Sebastiano Boscarino$^{(2)}$, Giovanni Russo$^{(2)}$

(1) CMAP, CNRS, École polytechnique, Institut Polytechnique de Paris
Route de Saclay, 91120 Palaiseau, France
giuseppe.orlando@polytechnique.edu

(2) Department of Mathematics and Computer Science, University of Catania
95125 Catania, Italy
sebastiano.boscarino@unict.it, giovanni.russo1@unict.it

**Keywords:** Asymptotic-preserving, Asymptotically-accurate, Euler equations, IMEX, Semi-implicit, Discontinuous Galerkin, Non-ideal gas.
Abstract

We present a quantitative comparison between two different Implicit-Explicit Runge-Kutta
(IMEX-RK) approaches for the Euler equations of gas dynamics, specifically tailored for the low
Mach limit. In this regime, a classical IMEX-RK approach involves an implicit coupling between
the momentum and energy balance so as to avoid the acoustic CFL restriction, while the density
can be treated in a fully explicit fashion. This approach leads to a mildly nonlinear equation
for the pressure, which can be solved according to a fixed point procedure. An alternative
strategy consists of employing a semi-implicit temporal integrator based on IMEX-RK methods
(SI-IMEX-RK). The stiff dependence is carefully analyzed, so as to avoid the solution of a
nonlinear equation for the pressure also for equations of state (EOS) of non-ideal gases. The
spatial discretization is based on a Discontinuous Galerkin (DG) method, which naturally allows
high-order accuracy. The asymptotic-preserving (AP) and the asymptotically-accurate (AA)
properties of the two approaches are assessed on a number of classical benchmarks for ideal
gases and on their extension to non-ideal gases.
# 1 Introduction

The Euler equations of gas dynamics represent the standard mathematical model for a wide range of applications in fluid mechanics, mechanical engineering, and environmental engineering [76]. Several numerical methods have been developed over the years, which can generally be divided into two categories, depending on a dimensionless parameter called Mach number. The Mach number $M$ represents the ratio between the local fluid velocity and the speed of sound in the medium [44]. For moderate to high Mach numbers, compressible effects have to be taken into account and numerical discretization strategies typically rely on Godunov-type shock-capturing schemes [42, 49, 63, 65, 82, 86]. On the other hand, in the low Mach number regime, the flow can be considered weakly compressible or even incompressible.

Explicit time discretization methods are very popular for high Mach number flows [65, 68, 82]. The time step has to satisfy a Courant-Friedrichs-Lewy (CFL) condition, given by the mesh size divided by the fastest wave speed [17, 33]. For moderate to high Mach numbers, this restriction is not a problem, since one is interested in resolving all the waves. However, for flows characterized by low values of the Mach number, severe time step restrictions may be required by these schemes. In this regime, acoustic waves usually carry a negligible amount of energy and, therefore, one may not be interested in resolving them. Hence, the system becomes stiff and stability limitations on the time step are much stricter than the restrictions imposed by accuracy.

The use of implicit and semi-implicit time discretization methods, so as to avoid the acoustic CFL restriction, has a long tradition in low Mach number flows [7, 17, 25, 58, 71, 72, 74]. Several numerical methods for weakly compressible flows have been proposed in the literature, see, e.g., [27, 32, 35, 61, 67, 85] and the references therein. Since the seminal paper [25], an effective approach to deal with low Mach flows is given by pressure-based algorithms. Indeed, an implicit treatment of the pressure gradient term within the momentum equation and of the pressure work term in the energy equation is sufficient to remove the acoustic CFL restriction and to decouple acoustic and transport effects [25].

Another aspect to consider is that weakly compressible flows are characterized by multiple length and time scales. Under specific assumptions, the compressible Euler equations converge to the incompressible ones as the Mach number goes to zero [38, 57]. Hence, the density remains constant along the fluid particle trajectories and the pressure acts as a Lagrange multiplier to enforce incompressibility of the flow [17, 58]. Robust numerical methods for the Euler equations should recover the incompressible limit for vanishing Mach number. For this purpose, the concept of asymptotic-preserving (AP) schemes has been introduced [47]. We also refer to [16] for a review of asymptotic-preserving methods for quasilinear hyperbolic systems with stiff relaxation. A numerical method for the compressible Euler equations is said to be asymptotic-preserving if its stability condition does not depend on the Mach number $M$ and if it provides a consistent discretization of the incompressible Euler equations as $M \to 0$. The AP property of the aforementioned approach was proven in [74].

Preserving incompressibility and resolving vortex dynamics are among the main purposes of numerical discretizations for weakly-compressible flows and high-order methods can help to reach this goal. The aim of the present work is to provide a quantitative comparison between two different Implicit-Explicit Runge-Kutta (IMEX-RK) approaches for the Euler equations, specifically tailored for the low Mach limit. The first approach uses an IMEX-RK solver, as proposed in [71] and validated in [72] for atmospheric applications. A key feature of this approach is the implicit treatment of acoustic waves, while material waves are handled explicitly. This method combines IMEX-RK schemes, carefully designed for stability and accuracy, with a time-stepping size that is independent of the Mach number $M$. The spatial discretization is based on the Discontinuous Galerkin (DG) method [40], which naturally allows for high-order accuracy and has proven highly effective for a wide range of computational fluid dynamics problems, across various flow regimes (see, e.g., [11, 30, 31, 54]). Additionally, the IMEX-DG method can handle a general equation of state (EOS) [71, 74], for which only a few studies have
been devoted [1, 32], and that, as we will discuss later, poses additional numerical challenges. However, this approach leads to a mildly nonlinear equation for the pressure also for ideal gases. In order to avoid this nonlinearity, an alternative strategy consists of employing a semi-implicit temporal integrator based on IMEX-RK schemes (SI-IMEX-RK), similar to the one adopted in [13]. This method leads to a linearized equation for both the pressure and the EOS. We refer to this scheme as the SI-IMEX-DG.

Furthermore, in this work, we show that the time discretization methods satisfy the AP property and the asymptotically accuracy (AA) property, i.e. they maintain their high-order accuracy also in the case of $M \to 0$. Moreover, if specific boundary conditions are considered, the limit model can differ from the incompressible Euler equations and depends on the employed EOS [74]. Specific numerical treatments for non-ideal EOS in the framework of the IMEX-DG scheme were presented in [71]. In [19], for non-ideal EOS, a nonlinear equation is solved through a Newton method. We propose here a novel different strategy, so as to avoid any nonlinear equation for the semi-implicit approach.

The paper is structured as follows. In Section 2, we briefly recall the mathematical model and its limit as $M \to 0$. In Section 3, we present the numerical method. More specifically, we outline the IMEX and the SI-IMEX time discretization methods. Moreover, we provide suitable strategies to deal with a general class of EOS. Some details of the DG formulation will also be discussed, specifying some advantages and disadvantages of this method for low Mach number flows. In Section 4, some numerical results to assess the properties of the two methods and to compare them are presented. Finally, some conclusions and perspectives for future work are discussed in Section 5.

## 2 The mathematical model

Let $\Omega \subset \mathbb{R}^d$, $1 \le d \le 3$ be a connected open bounded set with a sufficiently smooth boundary $\partial\Omega$ and denote by $\mathbf{x}$ the spatial coordinates and by $t$ the temporal coordinate. The mathematical model consists of the fully compressible Euler equations of gas dynamics, written in non-dimensional form as follows [17, 19, 58, 74]:

$$ 
\begin{aligned}
& \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0 \\
& \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \frac{1}{M^2} \nabla p = \mathbf{0} \\
& \frac{\partial \rho E}{\partial t} + \nabla \cdot [(\rho E + p) \mathbf{u}] = 0.
\end{aligned}
\quad (1) $$

Here, $\rho$ is the density, $\mathbf{u}$ is the fluid velocity, $p$ is the pressure, and $E$ is the total energy per unit of mass. Moreover, $M \equiv u_0/\sqrt{p_0/\rho_0}$, where $u_0$, $p_0$ and $\rho_0$ are reference fluid speed, pressure and density, respectively. It is related to the Mach number $M_0$, i.e. the ratio between a typical fluid velocity and a typical speed of sound. For a $\gamma$-law gas, for example, it is $M = M_0\sqrt{\gamma}$. We are mainly interested in the so-called low Mach regime, i.e. $M \ll 1$, for which material waves are much slower than acoustic waves. The previous set of equations has to be completed by an equation of state (EOS). Further details on the EOS will be discussed in the upcoming Section 2.1. The total energy $\rho E$ can be rewritten as $\rho E = \rho e + M^2\rho k$, where $e$ denotes the internal energy and $k = |\mathbf{u}|^2/2$ the kinetic energy. For the sake of convenience, we also introduce the specific enthalpy $h = e + p/\rho$, and we notice that the energy flux can be rewritten as

$$ (\rho E + p) \mathbf{u} = \left( e + M^2 k + \frac{p}{\rho} \right) \rho \mathbf{u} = (h + M^2 k) \rho \mathbf{u}. \quad (2) $$
Hence, system (1) reads equivalently as follows:

$$
\begin{align}
& \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0 \nonumber \\
& \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \frac{1}{M^2} \nabla p = \mathbf{0} \tag{3} \\
& \frac{\partial \rho E}{\partial t} + \nabla \cdot [(h + M^2 k) \rho \mathbf{u}] = 0. \nonumber
\end{align}
$$

2.1 The equation of state

System (3) has to be completed with an equation of state (EOS). In this work, we focus on the
classical ideal gas law, the stiffened gas EOS (SG-EOS) [62], and the general cubic EOS [87].
The equation that links together pressure, density, and internal energy for an ideal gas is given
by [87]

$$
p = (\gamma - 1) \rho e = (\gamma - 1) \left( \rho E - \frac{1}{2} M^2 \rho \mathbf{u} \cdot \mathbf{u} \right). \quad (4)
$$

Notice that (4) is valid only for a constant value of the specific heats ratio γ [87]. The analogous relation for the SG-EOS reads as follows:

$$
p = (\gamma - 1)(\rho e - \rho q_{\infty}) - \gamma \pi_{\infty} = (\gamma - 1) \left( \rho E - \frac{1}{2} M^2 \rho \mathbf{u} \cdot \mathbf{u} - \rho q_{\infty} \right) - \gamma \pi_{\infty}, \quad (5)
$$

with $q_\infty$ and $\pi_\infty$ representing constant parameters that determine the characteristics of the
fluid. Notice that, for $q_\infty = \pi_\infty = 0$ in (5), we recover (4). The last relation that we consider is
the so-called general cubic EOS, for which the link between pressure, density, and temperature
can be expressed as follows [81, p. 221], [87, p. 119]:

$$
p = \frac{\rho R_g T}{1 - \rho b} - \frac{a(T)\rho^2}{(1 - \rho b r_1)(1 - \rho b r_2)}. \quad (6)
$$

After some algebraic manipulations, (6) can be expressed as [81, p. 222]

$$
z^3 + c_0 z^2 + c_1 z + c_2 = 0, \tag{7}
$$

which is a cubic polynomial for the compressibility factor $z = p/(\rho R_g T)$. The parameters $c_0, c_1$, and $c_2$ depend on the thermodynamic state. Moreover, $R_g = R/m^*$ denotes the specific gas constant, with $R$ being the gas constant and $m^*$ the molar mass of the gas. Notice that for $a=b=0$, the expression of the pressure of an ideal gas is recovered. As discussed in [81, 87], it is convenient to express thermodynamic functions such as the internal energy or the enthalpy as the sum of a contribution due to the ideal gas and a *residual* contribution due to non-ideality. Hence, after some manipulations, the equation linking together internal energy, density, and temperature, is given by [81, p. 231], [87, p. 116]

$$
e = e^{\#}(T(p, \rho)) + \frac{a(T(p, \rho)) - T(p, \rho) \frac{da(T(p, \rho))}{dT}}{b} \frac{1}{r_1 - r_2} \log \left( \frac{1 - \rho br_1}{1 - \rho br_2} \right). \quad (8)
$$

Here, $e^{\#}(T(p, \rho))$ denotes the internal energy of an ideal gas, which is function solely of the
temperature $T$, $r_1$ and $r_2$ are suitable constants, whereas the parameters $a(T)$ and $b$ determine
fluid characteristics [87]. More specifically, $a(T)$ is related to intermolecular forces, whereas $b$,
the so called co-volume, takes into account the volume occupied by the molecules. Notice that,
for $r_1 \to 0$ and $r_2 \to 0$, then $\frac{1}{r_1-r_2} \log(\frac{1-\rho br_1}{1-\rho br_2}) \to -\rho b$, which corresponds to the van der Waals
EOS. For $r_1 = -1 - \sqrt{2}, r_2 = -1 + \sqrt{2}$, we get the Peng-Robinson EOS [77], [87, p. 118]. If
$c_v = \frac{de^{\#}}{dT}$ is constant, as in the case of calorically perfect gas, relation (8) can be simplified

$$
e = c_v T(p, \rho) + \frac{a(T(p, \rho)) - T(p, \rho) \frac{da(T(p, \rho))}{dT}}{b} \frac{1}{r_1 - r_2} \log \left( \frac{1 - \rho b r_1}{1 - \rho b r_2} \right). \quad (9)
$$
## 2.2 Asymptotic expansion

In this Section, we analyze the asymptotic limit of (3) as $M \to 0$. Consider the following expansion for density, velocity, and pressure, respectively:

$$ \rho(\mathbf{x}, t) = \bar{\rho}(\mathbf{x}, t) + M\rho'(\mathbf{x}, t) + M^2\rho''(\mathbf{x}, t) + o(M^2) \quad (10) $$

$$ \mathbf{u}(\mathbf{x}, t) = \bar{\mathbf{u}}(\mathbf{x}, t) + M\mathbf{u}'(\mathbf{x}, t) + M^2\mathbf{u}''(\mathbf{x}, t) + o(M^2) \quad (11) $$

$$ p(\mathbf{x}, t) = \bar{p}(\mathbf{x}, t) + Mp'(\mathbf{x}, t) + M^2p''(\mathbf{x}, t) + o(M^2) \quad (12) $$

From now on, to simplify the notation, we omit the explicit dependence on space and time for all the variables. Substituting (10) and (11) into the continuity equation in (3), the leading order term relation is

$$ \frac{\partial \bar{\rho}}{\partial t} + \nabla \cdot (\bar{\rho} \bar{\mathbf{u}}) = 0. \quad (13) $$

The leading order term relation for the momentum balance reduces to

$$ \nabla \bar{p} = \mathbf{0}. \quad (14) $$

Hence, $\bar{p}$ is a function solely of time. Analogously, for the first order term, we obtain

$$ \nabla p' = \mathbf{0}. \quad (15) $$

In addition, the second order term relation reads as follows:

$$ \frac{\partial \bar{\rho} \bar{\mathbf{u}}}{\partial t} + \nabla \cdot (\bar{\rho} \bar{\mathbf{u}} \otimes \bar{\mathbf{u}}) + \nabla p'' = \mathbf{0}, \quad (16) $$

where $p''$ represents a dynamical pressure [32, 58, 85]. Finally, the leading order term relation for the energy equation reads as follows:

$$ \frac{\partial \bar{\rho} \bar{e}}{\partial t} + \nabla \cdot (\bar{\rho} \bar{h} \bar{\mathbf{u}}) = 0, \quad (17) $$

where $\bar{e} = e(\bar{\rho}, \bar{p})$ and $\bar{h} = h(\bar{\rho}, \bar{p})$. Since $\bar{\rho}\bar{e} = \bar{\rho}\bar{h} - \bar{p}$, we obtain

$$ \frac{\partial \bar{\rho} \bar{h}}{\partial t} - \frac{\partial \bar{p}}{\partial t} + \bar{\mathbf{u}} \cdot \nabla (\bar{\rho} \bar{h}) + \bar{\rho} \bar{h} (\nabla \cdot \mathbf{u}) = 0, \quad (18) $$

or equivalently, considering $\bar{\rho}\bar{h} = (\bar{\rho}\bar{h})(\bar{\rho}, \bar{p})$,

$$ \frac{\partial \bar{\rho} \bar{h}}{\partial \bar{\rho}} \left( \frac{\partial \bar{\rho}}{\partial t} + \bar{\mathbf{u}} \cdot \nabla \bar{\rho} \right) + \frac{\partial \bar{\rho} \bar{h}}{\partial \bar{p}} \left( \frac{\partial \bar{p}}{\partial t} + \bar{\mathbf{u}} \cdot \nabla \bar{p} \right) - \frac{\partial \bar{p}}{\partial t} + \bar{\rho} \bar{h} (\nabla \cdot \bar{\mathbf{u}}) = 0. \quad (19) $$

Thanks to (13) and (14), we obtain

$$ \left( \bar{\rho} \bar{h} - \frac{\partial \bar{\rho} \bar{h}}{\partial \bar{\rho}} \bar{\rho} \right) (\nabla \cdot \bar{\mathbf{u}}) + \left( \frac{\partial \bar{\rho} \bar{h}}{\partial \bar{p}} - 1 \right) \frac{d\bar{p}}{dt} = 0, \quad (20) $$

or, since

$$ \bar{\rho} \bar{h} - \frac{\partial \bar{\rho} \bar{h}}{\partial \bar{\rho}} \bar{\rho} = - \bar{\rho}^2 \frac{\partial \bar{h}}{\partial \bar{\rho}} $$

and

$$ \frac{\partial \bar{\rho} \bar{h}}{\partial \bar{p}} - 1 = \frac{\partial \bar{\rho} \bar{e}}{\partial \bar{p}}, $$

equivalently,

$$ -\bar{\rho}^2 \frac{\partial \bar{h}}{\partial \bar{\rho}} (\nabla \cdot \bar{\mathbf{u}}) + \frac{\partial \bar{\rho} \bar{e}}{\partial \bar{p}} \frac{d\bar{p}}{dt} = 0. \quad (21) $$
We assume $\frac{\partial \bar{h}}{\partial \bar{\rho}} \neq 0$, as it holds away from vacuum. If $\frac{d\bar{p}}{dt} = 0$, we recover the incompressibility constraint

$$ \nabla \cdot \bar{\mathbf{u}} = 0. \qquad (22) $$

Summing up, the asymptotic limit of (3) is

$$ \begin{aligned} \frac{\partial \bar{\rho}}{\partial t} + \nabla \cdot (\bar{\rho} \bar{\mathbf{u}}) &= 0 \\ \nabla \bar{p} &= \mathbf{0} \\ \nabla p' &= \mathbf{0} \end{aligned} \qquad (23) $$

$$ \begin{aligned} \frac{\partial \bar{\rho} \bar{\mathbf{u}}}{\partial t} + \nabla \cdot (\bar{\rho} \bar{\mathbf{u}} \otimes \bar{\mathbf{u}}) + \nabla p'' &= \mathbf{0} \\ -\bar{\rho}^2 \frac{\partial \bar{h}}{\partial \bar{\rho}} (\nabla \cdot \bar{\mathbf{u}}) + \frac{\partial \bar{\rho} \bar{e}}{\partial \bar{p}} \frac{\partial \bar{p}}{\partial t} &= 0. \end{aligned} $$

System (23) represents a more general asymptotic limit with respect to the EOS and the boundary conditions of the compressible Euler equations for vanishing Mach number [34, 74]. We notice that we can rewrite the last equation of system (23) as follows:

$$ \nabla \cdot \bar{\mathbf{u}} = \left( \frac{\partial \bar{\rho} \bar{e}}{\partial \bar{p}} \right) \left( \bar{\rho}^2 \frac{\partial \bar{h}}{\partial \bar{\rho}} \right)^{-1} \frac{d\bar{p}}{dt} = - \frac{1}{\bar{\rho} c^2 (\bar{\rho}, \bar{p})} \frac{d\bar{p}}{dt}, \qquad (24) $$

where c denotes the speed of sound. Indeed, from the definition of $\bar{h}$, one has

$$ \left(\bar{\rho}^2 \frac{\partial \bar{h}}{\partial \bar{\rho}}\right) \left(\frac{\partial \bar{\rho} \bar{e}}{\partial \bar{p}}\right)^{-1} = \rho \left(\frac{\partial \bar{e}}{\partial \bar{\rho}} - \frac{\bar{p}}{\bar{\rho}^2}\right) \left(\frac{\partial \bar{e}}{\partial \bar{p}}\right)^{-1}. $$

From the first principle of thermodynamics, denoting by s the specific entropy, one has

$$ T ds = de + pd \frac{1}{\rho} = \frac{\partial e}{\partial \rho} d\rho + \frac{\partial e}{\partial p} dp - \frac{p}{\rho^2} d\rho = \left( \frac{\partial e}{\partial \rho} - \frac{p}{\rho^2} \right) d\rho + \frac{\partial e}{\partial p} dp. $$

As $c^2 = \frac{\partial p}{\partial \rho}|_s$, one has

$$ c^2 = \left. \frac{\partial p}{\partial \rho} \right|_s = \left( \frac{\partial e}{\partial p} \right)^{-1} \left( \frac{p}{\rho^2} - \frac{\partial e}{\partial \rho} \right) = - \left( \frac{\partial e}{\partial p} \right)^{-1} \frac{\partial h}{\partial \rho}, \qquad (25) $$

and therefore

$$ \frac{1}{\rho c^2} = - \frac{1}{\rho} \left( \frac{\partial e}{\partial p} \right)^{-1} \frac{\partial h}{\partial \rho}, $$

which proves (24) when applied to the lowest order terms in the asymptotic expansions (10), (12). For more details consult [74, 87, 89]. Under periodic or free-slip boundary conditions, thanks to the divergence theorem, we have

$$ \int_{\Omega} \nabla \cdot \bar{\mathbf{u}} d\Omega = 0, $$

so that, by integrating (24) on $\Omega$, we find $\frac{d\bar{p}}{dt} = 0$. On the other hand, as one can easily notice from (24), a time-dependent pressure with large amplitude variations imposed by a Dirichlet outflow boundary condition leads to a non-incompressible flow, i.e. $\nabla \cdot \bar{\mathbf{u}} \neq 0$ and depending on the specific EOS. Consider, e.g., the ideal gas law (4). We get

$$ \frac{\partial \bar{\rho} \bar{e}}{\partial \bar{p}} = \frac{1}{\gamma - 1} \qquad \bar{\rho}^2 \frac{\partial \bar{h}}{\partial \bar{\rho}} = - \frac{\gamma}{\gamma - 1} \bar{p}, \qquad (26) $$
so that (24) reduces to

$$ \nabla \cdot \bar{\mathbf{u}} = - \frac{1}{\gamma} \frac{d \log \bar{\rho}}{dt}. \qquad (27) $$

Hence, the compressibility of a fluid described by the ideal gas law (4) is uniform in space and changes only in time. This is no longer valid for a general EOS [74].

# 3 The numerical framework

In the low Mach number limit, pressure gradients terms, which are proportional to $1/M^2$, yield stiff components for the resulting semi-discretized ODE system [25, 66, 71]. Implicit-Explicit Runge-Kutta (IMEX-RK) methods [13, 55] are widely employed for ODE systems that include both stiff and non-stiff components, to which the implicit and explicit schemes are applied, respectively. Therefore, an implicit coupling between the energy equation and the momentum one is appropriate, while the continuity equation can be treated in a fully explicit fashion. The spatial discretization is based on the Discontinuous Galerkin (DG) method, which easily allows for high-order accuracy. We refer to [40] for a general introduction to the method. In this Section, we review some well-known concepts of IMEX-RK schemes. Then, we present the IMEX and the SI-IMEX time discretization for system (3), respectively. Finally, some details concerning the spatial discretization will be also discussed.

## 3.1 IMEX Runge-Kutta schemes

Implicit-Explicit Runge-Kutta (IMEX-RK) methods find extensive application in the numerical solution of PDEs, such as hyperbolic systems with relaxations [14, 75], convection-diffusion equations [8] and convection-diffusion-reaction equations [55, 56]. Let us start considering the following initial value problem for a system of ODE's

$$ \frac{dy}{dt} = \mathbf{f}_E(y, t) + \mathbf{f}_I(y, t), \quad y(0) = y_0, \qquad (28) $$

where $y(t) \in \mathbb{R}^M, M \ge 1$ and we assume that $\mathbf{f}_E$ and $\mathbf{f}_I: \mathbb{R}^M \times \mathbb{R} \to \mathbb{R}^M$ are Lipschitz functions of $y(t)$. We assume that the term $\mathbf{f}_I$ is stiff and the term $\mathbf{f}_E$ is non-stiff.

An s-stage IMEX-RK scheme applied to system (28) takes the form:

$$ \mathbf{v}^{(l)} = \mathbf{v}^n + \Delta t \sum_{m=1}^{l-1} \tilde{a}_{lm} \mathbf{f}_E \left( t^n + \tilde{c}_m \Delta t, \mathbf{v}^{(m)} \right) + \Delta t \sum_{m=1}^{s} a_{lm} \mathbf{f}_I \left( t^n + c_m \Delta t, \mathbf{v}^{(m)} \right), (29a) $$

$$ \mathbf{v}^{n+1} = \mathbf{v}^n + \Delta t \sum_{l=1}^{s} \tilde{b}_l \mathbf{f}_E \left( t^n + \tilde{c}_l \Delta t, \mathbf{v}^{(l)} \right) + \Delta t \sum_{l=1}^{s} b_l \mathbf{f}_I \left( t^n + c_l \Delta t, \mathbf{v}^{(l)} \right). \qquad (29b) $$

where the quantities $\mathbf{v}^{(l)}$ for $l = 1, \dots, s$, are called internal stages and approximate the exact solution $y(t)$, at time $t = t^n + c_l \Delta t$, whereas $\mathbf{v}^{n+1}$ is the numerical solution that approximates the exact solution $y(t)$ at time $t = t^n + \Delta t$. An s-stage IMEX-RK method is defined by two $s \times s$ real matrices $\tilde{\mathbf{A}} = \{\tilde{a}_{lm}\}$ and $\mathbf{A} = \{a_{lm}\}$, where the matrix $\tilde{\mathbf{A}}$ corresponds to the explicit method and is a lower triangular matrix with zero diagonal, i.e., $\tilde{a}_{lm} = 0$ for $l \le m$, while $\mathbf{A}$ is the one corresponding to the implicit scheme. We consider Diagonally Implicit Runge-Kutta (DIRK) methods for the implicit scheme so that $a_{lm} = 0$ for $l < m$. The use of a DIRK method for the treatment of $\mathbf{f}_I$ provides a sufficient condition to guarantee that the function $\mathbf{f}_E$ is always evaluated explicitly. The method is also characterized by the quadrature nodes $\tilde{c} = (0, \tilde{c}_2, \dots, \tilde{c}_s)^\top$, $\mathbf{c} = (c_1, c_2, \dots, c_s)^\top$, given by the usual relation

$$ \sum_{m=1}^{s} a_{lm} = c_l \qquad \sum_{m=1}^{s} \tilde{a}_{lm} = \tilde{c}_l, \quad l = 1, \dots, s, \qquad (30) $$
and by the weights $\tilde{\mathbf{b}}^{\top} = (\tilde{b}_{1}, \tilde{b}_{2}, \dots, \tilde{b}_{s})$ and $\mathbf{b}^{\top} = (b_{1}, b_{2}, \dots, b_{s})$ in $\mathbb{R}^{s}$. IMEX-RK methods can be represented in the usual Butcher notation [23]

$$ \begin{array}{c|c} \tilde{\mathbf{c}} & \tilde{\mathbf{A}} \\ \hline \tilde{\mathbf{b}}^{\top} & \mathbf{c} \\ \end{array} \quad \begin{array}{c|c} \mathbf{A} & \mathbf{b}^{\top} \\ \hline \end{array} $$

Notice that the relation (30) is a usual assumption for Runge-Kutta methods [48].

It is useful to characterize the different IMEX-RK methods presented in the literature in two main types according to the structure of the matrix of the DIRK method. Following [14], we have

**Definition 3.1.** An IMEX-RK method is said to be of **type I** [13, 75] if the matrix **A** is invertible. It is said to be of **type II** [13, 55] if the matrix **A** can be written in the form

$$ \mathbf{A} = \begin{pmatrix} 0 & 0 \\ \mathbf{a} & \mathbf{A} \end{pmatrix}, $$

with $\mathbf{a} = (a_{21}, \dots, a_{s1})^{\top} \in \mathbb{R}^{s-1}$ and the matrix $\mathbf{A} \in \mathbb{R}^{(s-1) \times (s-1)}$ is invertible. In the special case $\mathbf{a} = 0, b_1 = 0$, the method is said of **type ARS** (see [8]) and the DIRK method is reducible to a method using $s-1$ stages.

Schemes of type II allow some simplifying assumptions, that make order conditions easier to treat for the construction of higher order schemes [56]. On the other hand, schemes of type I are more suited to a theoretical analysis [12, 15] because of the invertibility of **A**.

**Definition 3.2.** We call an IMEX-RK method stiffly accurate (SA) if the corresponding DIRK method is stiffly accurate, namely [88]

$$ a_{sl} = b_l, \quad i = 1, \dots, s. \tag{31} $$

All the IMEX-RK schemes employed in this work for the numerical simulations are stiffly accurate (see Appendix A).

The asymptotic properties of IMEX-RK methods are strongly related to the L-stability of the implicit part of the scheme. An implicit Runge-Kutta scheme is said to be L-stable [88] if it is **A**-stable and $R(z) \to 0$ as $z \to \infty$, where $R(z)$ is the stability function of the DIRK scheme. Following the result in [88], L-stability is typically obtained combining the **A**-stability property with the SA property. However, for methods of type II, this combination does not necessarily lead to a L-stable scheme for the implicit part, because the matrix **A** is not invertible [15]. For SA schemes of type II, a supplementary condition is required to obtain L-stability, i.e. [15]

$$ \mathbf{e}_s^{\top} \mathbf{A}^{-1} \mathbf{a} = \sum_{m=2}^{s} \hat{w}_{sm} a_{m1} = 0, \tag{32} $$

where $\mathbf{e}_s^{\top} = (0, \dots, 0, 1)^{\top}$ and $\hat{w}_{lm}$ denotes the elements of the inverse of $\mathbf{A}$. One can easily verify that all the implicit companion methods reported in Appendix A are L-stable.

In the sequel, to identify the different IMEX-RK schemes, we shall use the notation $(s, \sigma, p)$, where $s$ is the number of function evaluations of the implicit companion method, $\sigma$ is the number of function evaluations of the explicit companion method, and $p$ is the order of the IMEX scheme. In this work, we employ second, third, and fourth order time discretization schemes (see Appendix A).

## 3.2 IMEX time discretization for the Euler equations

In this Section, we outline the IMEX time discretization for the Euler equations (3). Following [25, 36], we consider an implicit treatment of the pressure gradient term within the momentum
equation and of the pressure work term in the energy equation, while the continuity equation is discretized in a fully explicit fashion. A generic stage reads therefore as follows [71, 74]:

$$ \rho^{(l)} = \rho^n - \Delta t \sum_{m=1}^{l-1} \tilde{a}_{lm} \nabla \cdot (\mathbf{q})^{(m)} $$

$$ \mathbf{q}^{(l)} + \frac{1}{M^2} a_{ll} \Delta t \nabla p^{(l)} = \mathbf{q}^n - \frac{\Delta t}{M^2} \sum_{m=1}^{l-1} a_{lm} \nabla p^{(m)} - \Delta t \sum_{m=1}^{l-1} \tilde{a}_{lm} \nabla \cdot (\mathbf{q} \otimes \mathbf{u})^{(m)} \quad (33) $$

$$ \mathcal{E}^{(l)} + a_{ll} \Delta t \nabla \cdot (h\mathbf{q})^{(l)} = \mathcal{E}^{n} - \Delta t \sum_{m=1}^{l-1} a_{lm} \nabla \cdot (h\mathbf{q})^{(m)} - \Delta t M^2 \sum_{m=1}^{l-1} \tilde{a}_{lm} \nabla \cdot (k\mathbf{q})^{(m)}, $$

where

$$ \mathbf{q}^{(l)} = (\rho \mathbf{u})^{(l)} \qquad \mathbf{u}^{(l)} = \frac{\mathbf{q}^{(l)}}{\rho^{(l)}} \qquad \mathcal{E}^{(l)} = (\rho E)^{(l)}. \quad (34) $$

Notice that, substituting formally $\mathbf{q}^{(l)}$ into the energy equation and taking into account the definitions $\rho E = \rho e + M^2 \rho k$ and $h = e + p/\rho$, the following nonlinear Helmholtz-type equation for the pressure is obtained:

$$ \rho^{(l)} \left[ e(p^{(l)}, \rho^{(l)}) + M^2 k^{(l)} \right] - a_{ll}^2 \frac{\Delta t^2}{M^2} \nabla \cdot \left[ \left( e(p^{(l)}, \rho^{(l)}) + \frac{p^{(l)}}{\rho^{(l)}} \right) \nabla p^{(l)} \right] \\ + a_{ll} \Delta t \nabla \cdot \left[ \left( e(p^{(l)}, \rho^{(l)}) + \frac{p^{(l)}}{\rho^{(l)}} \right) \mathbf{m}^{(l)} \right] = \hat{e}^{(l)}, \quad (35) $$

where

$$ \mathbf{m}^{(l)} = \mathbf{q}^n - \frac{\Delta t}{M^2} \sum_{m=1}^{l-1} a_{lm} \nabla p^{(m)} - \Delta t \sum_{m=1}^{l-1} \tilde{a}_{lm} \nabla \cdot (\mathbf{q} \otimes \mathbf{u})^{(m)}, \quad (36a) $$

$$ \hat{e}^{(l)} = \mathcal{E}^{n} - \Delta t \sum_{m=1}^{l-1} a_{lm} \nabla \cdot (h\mathbf{q})^{(m)} - M^2 \Delta t \sum_{m=1}^{l-1} \tilde{a}_{lm} \nabla \cdot (k\mathbf{q})^{(m)}. \quad (36b) $$

Equation (35) is solved through a fixed point procedure [36, 71]. More specifically, setting $\xi^{(0)} = p^{(l-1)}, k^{(l,0)} = k^{(l-1)}$, one solves for $\tilde{k} = 0, \dots, L$ the equation

$$ \begin{aligned} \rho^{(l)} e(\xi^{(\tilde{k}+1)}, \rho^{(l)}) & - a_{ll}^2 \frac{\Delta t^2}{M^2} \nabla \cdot \left[ \left( e(\xi^{(\tilde{k})}, \rho^{(l)}) + \frac{\xi^{(\tilde{k})}}{\rho^{(l)}} \right) \nabla \xi^{(\tilde{k}+1)} \right] \\ & = \hat{e}^{(l)} - M^2 \rho^{(l)} k^{(l,\tilde{k})} - a_{ll} \Delta t \nabla \cdot \left[ \left( e(\xi^{(\tilde{k})}, \rho^{(l)}) + \frac{\xi^{(\tilde{k})}}{\rho^{(l)}} \right) \mathbf{m}^{(l)} \right] \end{aligned} $$

and then updates the velocity as

$$ \mathbf{u}^{(l, \tilde{k}+1)} + \frac{a_{ll} \Delta t}{\rho^{(l)} M^2} \nabla \xi^{(\tilde{k}+1)} = \frac{\mathbf{m}^{(l)}}{\rho^{(l)}}. $$

As already discussed in [36], solving directly (35) keeping a full implicit treatment of the enthalpy as in a classical Newton method yields a system strongly nonlinear and difficult to control. For this purpose, one adopts a Picard iteration technique in which the contribution of the enthalpy is computed at the previous fixed point iteration so as to reduce the nonlinearity of (35). Moreover, this choice is justified by the fact that two/three iterations are typically sufficient to obtain a satisfactory solution, as already observed in [26, 36] and as further confirmed by our numerical experiments (see in particular Section 4.2.1).
## 3.3 Semi-Implicit IMEX (SI-IMEX) time discretization for the Euler equations

The procedure outlined in the previous Section always requires the solution of a nonlinear system at each stage. In this Section, we outline the semi-implicit IMEX (SI-IMEX) time discretization for the Euler equations (3) [17], which provides similar results with less computational time. The governing partial differential equations (3) can be cast into a compact and general form as

$$ \frac{\partial \mathbf{U}}{\partial t} = \mathbf{H}(\mathbf{U}, \mathbf{U}). \qquad (37) $$

Here $\mathbf{U} = (\rho, \rho\mathbf{u}, \rho E)^{\top}$ is the vector of conserved variables, and

$$ \mathbf{H}(\mathbf{U}, \mathbf{U}) = -\nabla \cdot \begin{pmatrix} \rho \mathbf{u} \\ \rho \mathbf{u} \otimes \mathbf{u} \\ M^2 \rho k \mathbf{u} \end{pmatrix} - \nabla \cdot \begin{pmatrix} 0 \\ \frac{1}{M^2} p \\ h \rho \mathbf{u} \end{pmatrix}, \qquad (38) $$

where $\mathbf{H}: \mathbb{R}^{d+2} \times \mathbb{R}^{d+2} \rightarrow \mathbb{R}^{d+2}$ is a sufficiently regular mapping. Following [13], the governing partial differential equations (3) are written under the form of an autonomous system (37), for all $t > t_0$ with the initial condition $\mathbf{U}(t_0) = \mathbf{U}_0$.

We refer to semi-implicit schemes as numerical methods that address problems of the form (37), wherein the variable $\mathbf{U}$, appearing as first argument of $\mathbf{H}$, is treated explicitly and will be denoted by $\mathbf{U}_E$, while the variable $\mathbf{U}$ appearing as the second argument is treated implicitly and denoted by $\mathbf{U}_I$. Thus, we obtain a partitioned system of the form

$$ \begin{aligned} \frac{\partial \mathbf{U}_E}{\partial t} &= \mathbf{H}(\mathbf{U}_E, \mathbf{U}_I), \\ \frac{\partial \mathbf{U}_I}{\partial t} &= \mathbf{H}(\mathbf{U}_E, \mathbf{U}_I). \end{aligned} \qquad (39) $$

with

$$ \mathbf{H}(\mathbf{U}_E, \mathbf{U}_I) = -\nabla \cdot \mathbf{F}_E - \nabla \cdot \mathbf{F}_{SI}, \qquad (40) $$

and

$$ \mathbf{F}_E = \begin{pmatrix} (\rho \mathbf{u})_E \\ (\rho \mathbf{u} \otimes \mathbf{u})_E \\ M^2 (\rho k \mathbf{u})_E \end{pmatrix}, \quad \mathbf{F}_{SI} = \begin{pmatrix} 0 \\ \frac{1}{M^2} p_I \\ h_E (\rho \mathbf{u})_I \end{pmatrix}. \qquad (41) $$

Moreover, $\mathbf{U}_E = (\rho_E, (\rho \mathbf{u})_E, (\rho E)_E)^{\top}$ and $\mathbf{U}_I = (\rho_I, (\rho \mathbf{u})_I, (\rho E)_I)^{\top}$. Subscripts E and SI in (40) indicate the explicit and semi-implicit treatment of the first and the second term, respectively. Notice that the number of unknowns in (39) has been doubled. However, when specific time discretizations are chosen for autonomous systems, this doubling is only apparent [13]. Finally, the kinetic energy in the total energy definition splits into an explicit and an implicit contribution, namely:

$$ (\rho E)_I = (\rho e)_I + \frac{M^2}{2} \mathbf{u}_E \cdot (\rho \mathbf{u})_I. \qquad (42) $$

High-order time discretization is achieved making use of IMEX-RK schemes. More specifically, we adopt methods for which $\tilde{b}_l = b_l$, $l = 1, \dots, s$. We observe that, since $\tilde{b}_l = b_l$, the numerical solutions are the same, i.e., if $\mathbf{U}_E^0 = \mathbf{U}_I^0$, then $\mathbf{U}_E^n = \mathbf{U}_I^n$ for all $n > 0$. Hence, the duplication of the system is only apparent.

Under the assumption that system (39) is autonomous, a SI-IMEX-RK method is obtained as follows. First, set $\mathbf{U}_E^n = \mathbf{U}_I^n = \mathbf{U}^n$. Then, the internal stage values read

$$ \begin{aligned} \mathbf{U}_E^{(l)} &= \mathbf{U}_E^n + \Delta t \sum_{m=1}^{l-1} \tilde{a}_{lm} \mathbf{H} \left( \mathbf{U}_E^{(m)}, \mathbf{U}_I^{(m)} \right) \\ \mathbf{U}_I^{(l)} &= \mathbf{U}_I^n + \Delta t \sum_{m=1}^{l-1} a_{lm} \mathbf{H} \left( \mathbf{U}_E^{(m)}, \mathbf{U}_I^{(m)} \right) + \Delta t a_{ll} \mathbf{H} \left( \mathbf{U}_E^{(l)}, \mathbf{U}_I^{(l)} \right), \end{aligned} \qquad (43) $$
for $l = 1, \dots, s$. Finally, the numerical solution is updated with

$$ \mathbf{U}^{n+1} = \mathbf{U}^n + \Delta t \sum_{l=1}^{s} b_l \mathbf{H} \left( \mathbf{U}_E^{(l)}, \mathbf{U}_I^{(l)} \right). \quad (44) $$

For the sake of clarity in the notation, we denote

$$ \mathbf{q} = \rho \mathbf{u} \qquad \mathcal{E} = \rho E $$

as done in the previous Section. As an example, we present the first order semi-implicit scheme solving system (3) to compute the numerical solution $\mathbf{U}^{n+1} = (\rho^{n+1}, \mathbf{q}^{n+1}, \mathcal{E}^{n+1})^\top$. We focus on the time discretization, while keeping the space continuous. We consider the first order IMEX-RK scheme

$$ \begin{array}{c|c} 0 & 0 \\ \hline 1 & 1 \end{array} \quad \begin{array}{c|c} 1 & 1 \\ \hline 1 & 1 \end{array}. $$

Formally applying the above tableau to the partitioned system (39), it reads

$$ \begin{align} \mathbf{U}_E^{(1)} &= \mathbf{U}^n \nonumber \\ \mathbf{U}_I^{(1)} &= \mathbf{U}^n + \Delta t \mathbf{H} \left( \mathbf{U}^n, \mathbf{U}_I^{(1)} \right), \tag{45} \end{align} $$

$$ \mathbf{U}^{n+1} = \mathbf{U}_I^{(1)}, $$

and explicitly we get

$$ \mathbf{U}^{n+1} = \mathbf{U}^n - \Delta t \nabla \cdot \mathbf{F}_E(\mathbf{U}^n) - \Delta t \nabla \cdot \mathbf{F}_{SI}(\mathbf{U}^n, \mathbf{U}^{n+1}) \quad (46) $$

with

$$ \mathcal{E}^{n+1} = (\rho e)^{n+1} + \frac{M^2}{2} \mathbf{u}^n \cdot \mathbf{q}^{n+1}. $$

A generic stage $l$ for the Euler equations using a SI-IMEX-RK scheme reads therefore as follows.

**Explicit step.** Set

$$ \begin{align} \rho_E^{(l)} &= \rho^n - \Delta t \sum_{m=1}^{l-1} \tilde{a}_{lm} \nabla \cdot \mathbf{q}_E^{(m)} \nonumber \\ \mathbf{q}_E^{(l)} &= \mathbf{q}^n - \Delta t \sum_{m=1}^{l-1} \tilde{a}_{lm} \left[ \nabla \cdot (\mathbf{q} \otimes \mathbf{u})_E^{(m)} + \frac{1}{M^2} \nabla p_I^{(m)} \right] \tag{47} \end{align} $$

$$ \mathcal{E}_E^{(n,l)} = \mathcal{E}^n - \Delta t \sum_{m=1}^{l-1} \tilde{a}_{lm} \left[ M^2 \nabla \cdot (k\mathbf{q})_E^{(m)} + \nabla \cdot (h_E \mathbf{q}_I)^{(m)} \right]. $$

**Implicit step.** Solve $\mathbf{U}_I^{(l)}$:

$$ \begin{align} \rho_I^{(l)} &= \bar{\rho}_I^{(l)} - a_{ll} \Delta t \nabla \cdot \mathbf{q}_E^{(l)}, \nonumber \\ \mathbf{q}_I^{(l)} + \frac{1}{M^2} a_{ll} \Delta t \nabla p_I^{(l)} &= \mathbf{m}_I^{(l)} - a_{ll} \Delta t \nabla \cdot (\mathbf{q} \otimes \mathbf{u})_E^{(l)} \tag{48} \end{align} $$

$$ \mathcal{E}_I^{(l)} + a_{ll} \Delta t \nabla \cdot (h_E \mathbf{q}_I)^{(l)} = \hat{e}_I^{(l)} - a_{ll} \Delta t M^2 \nabla \cdot (k\mathbf{q})_E^{(l)}, $$
where

$$ \bar{\rho}_I^{(l)} = \rho^n - \Delta t \sum_{m=1}^{l-1} a_{lm} \nabla \cdot \mathbf{q}_E^{(m)} $$

$$ \mathbf{m}_I^{(l)} = \mathbf{q}^n - \Delta t \sum_{m=1}^{l-1} a_{lm} \left[ \nabla \cdot (\mathbf{q} \otimes \mathbf{u})_E^{(m)} + \frac{1}{M^2} \nabla p_I^{(m)} \right] \quad (49) $$

$$ \hat{e}_I^{(l)} = \mathcal{E}^n - \Delta t \sum_{m=1}^{l-1} a_{lm} \left[ M^2 \nabla \cdot (k\mathbf{q})_E^{(m)} - \nabla \cdot (h_E \mathbf{q}_I)^{(m)} \right]. $$

To solve system (48), we substitute $\mathbf{q}_I^{(l)}$ in the energy equation, so as to obtain an elliptic equation for $p_I^{(l)}$ [17], which reads as follows:

$$ \begin{aligned} \rho_I^{(l)} e_I(p_I^{(l)}, \rho_I^{(l)}) & - a_{ll}^2 \frac{\Delta t^2}{M^2} \nabla \cdot (h_E^{(l)} \nabla p_I^{(l)}) - \frac{1}{2} a_{ll} \Delta t \mathbf{u}_E^{(l)} \cdot \nabla p_I^{(l)} \\ & = \hat{e}_I^{(l)} - a_{ll} \Delta t M^2 \nabla \cdot (k \mathbf{q})_E^{(l)} - \frac{M^2}{2} \mathbf{u}_E^{(l)} \cdot \left[ \mathbf{m}_I^{(l)} - a_{ll} \Delta t \nabla \cdot (\mathbf{q} \otimes \mathbf{u})_E^{(l)} \right] \\ & - a_{ll} \Delta t \nabla \cdot \left[ h_E^{(l)} \left[ \mathbf{m}_I^{(l)} - a_{ll} \Delta t \nabla \cdot (\mathbf{q} \otimes \mathbf{u})_E^{(l)} \right] \right]. \end{aligned} \quad (50) $$

Next, one computes $\mathbf{q}_I^{(l)}$ and $\mathcal{E}_I^{(l)}$ from (48). Finally, one updates the numerical solution $\mathbf{U}^{n+1}$ from (44).

## 3.4 Impact of the EOS

Equations (37) and (50) need the relation between the internal energy and the pressure before being solved. In the case of the ideal gas law (4), since

$$ \rho e = \frac{1}{\gamma - 1} p, $$

equations (37) and (50) constitute a linear system for $\xi^{(\tilde{k}+1)}$ and $p_I^{(l)}$, respectively. Analogous considerations hold for the SG-EOS (5), since

$$ \rho e = \frac{p}{\gamma - 1} + \frac{\gamma \pi_{\infty}}{\gamma - 1} + \rho q_{\infty}. $$

Hence, only a supplementary term depending on the already updated density is present.

On the other hand, the use of a more general equation of state, such as the general cubic EOS (6), leads to a nonlinear relation between internal energy and pressure and therefore a nonlinear equation for the pressure should be solved [19]. We rewrite the term $\rho e$ as $\frac{\rho e}{p} p$, so that, following the discussion in [22], the nonlinear equation is solved by the following Picard iteration

$$ \begin{aligned} \frac{\rho^{(l)} e(\xi^{(\tilde{k})}, \rho^{(l)})}{\xi^{(\tilde{k})}} \xi^{(\tilde{k}+1)} & - a_{ll}^2 \frac{\Delta t^2}{M^2} \nabla \cdot \left[ \left( e(\xi^{(\tilde{k})}, \rho^{(l)}) + \frac{\xi^{(\tilde{k})}}{\rho^{(l)}} \right) \nabla \xi^{(\tilde{k}+1)} \right] \\ = \hat{e}^{(l)} & - M^2 \rho^{(l)} k^{(l, \tilde{k})} - a_{ll} \Delta t \nabla \cdot \left[ \left( e(\xi^{(\tilde{k})}, \rho^{(l)}) + \frac{\xi^{(\tilde{k})}}{\rho^{(l)}} \right) \mathbf{m}^{(l)} \right]. \end{aligned} \quad (51) $$

It is worth noting that in the case of the ideal gas law (4),

$$ \frac{\rho^{(l)} e(\xi^{(\tilde{k})}, \rho^{(l)})}{\xi^{(\tilde{k})}} \xi^{(\tilde{k}+1)} = \frac{1}{\gamma - 1} \xi^{(\tilde{k}+1)} = \rho^{(l)} e(\xi^{(\tilde{k}+1)}), \quad 13 $$
so that (51) reduces to (37). Notice also that (51) corresponds to a slightly different linearization with respect to the one proposed in [71], which was tailored for the general cubic EOS, while (51) is applicable to a general EOS.

Analogous considerations hold for (50): one can rewrite

$$ (\rho e)_I = \frac{(\rho e)_I}{p_I} p_I \qquad (52) $$

and solve the resulting mild nonlinear equation according to the Picard iteration described in [22]. Note that the fixed point procedure proposed in [22] corresponds to a Newton-type method which can also be applied to non-differentiable relations, like those that could be obtained for tabulated EOS, unlike the standard Newton method. We refer to [22] for a detailed description and analysis of the algorithm. Nevertheless, in the case of the semi-implicit time discretization, in order to avoid the solution of a nonlinear equation and the use of a fixed point loop, we approximate

$$ \frac{(\rho e)_I}{p_I} p_I \approx \frac{(\rho e)_E}{p_E} p_E, \qquad (53) $$

so that (50) modifies as

$$ \begin{aligned} \frac{(\rho e)_E^{(l)}}{p_E^{(l)}} p_I^{(l)} & - a_{ll}^2 \frac{\Delta t^2}{M^2} \nabla \cdot (h_E^{(l)} \nabla p_I^{(l)}) - \frac{1}{2} a_{ll} \Delta t \mathbf{u}_E^{(l)} \cdot \nabla p_I^{(l)} \\ & = \hat{e}_I^{(l)} - a_{ll} \Delta t M^2 \nabla \cdot (k \mathbf{q})_E^{(l)} - \frac{M^2}{2} \mathbf{u}_E^{(l)} \cdot [\mathbf{m}_I^{(l)} - a_{ll} \Delta t \nabla \cdot (\mathbf{q} \otimes \mathbf{u}_E^{(l)})] \\ & - a_{ll} \Delta t \nabla \cdot [h_E^{(l)} [\mathbf{m}_I^{(l)} - a_{ll} \Delta t \nabla \cdot (\mathbf{q} \otimes \mathbf{u}_E^{(l)})]]. \end{aligned} \qquad (54) $$

It is to be noted that, in the case of the ideal gas law (4),

$$ \frac{(\rho e)_E^{(l)}}{p_E^{(l)}} p_I^{(l)} = \frac{1}{\gamma - 1} p_I^{(l)} = (\rho e)_I^{(l)}, $$

so that (54) reduces to (50).

## 3.5 The spatial discretization strategy

In this Section, we briefly outline the spatial discretization adopted for (33) and (47)-(48), which is based on the Discontinuous Galerkin (DG) method [40] as implemented in the deal.II library [3, 9]. We consider a decomposition of the domain $\Omega$ into a family of quadrilaterals $T_h$ and denote each element by $K$. The skeleton $\mathcal{E} = \mathcal{E}^I \cup \mathcal{E}^B$ denotes the set of all the element faces, with $\mathcal{E}^I$ and $\mathcal{E}^B$ being the subset of interior and boundary faces, respectively. A face $\Gamma \in \mathcal{E}^I$ shares two elements, $K^+$ with outward unit normal $\mathbf{n}^+$, and $K^-$ with outward unit normal $\mathbf{n}^-$, while we simply denote by $\mathbf{n}$ the outward unit normal for a face $\Gamma \in \mathcal{E}^B$ (see Figure 1). For a scalar function $\varphi$, the jump is defined as

$$ \llbracket \varphi \rrbracket = \varphi^+ \mathbf{n}^+ + \varphi^- \mathbf{n}^- \quad \text{if } \Gamma \in \mathcal{E}^I \qquad \llbracket \varphi \rrbracket = \varphi \mathbf{n} \quad \text{if } \Gamma \in \mathcal{E}^B, \qquad (55) $$

while the average reads

$$ \{\{\varphi\}\} = \frac{1}{2}(\varphi^+ + \varphi^-) \quad \text{if } \Gamma \in \mathcal{E}^I \qquad \{\{\varphi\}\} = \varphi \quad \text{if } \Gamma \in \mathcal{E}^B. \qquad (56) $$

Analogous definitions apply for a vector function $\varphi$. More specifically, we define

$$ \llbracket \varphi \rrbracket = \varphi^+ \cdot \mathbf{n}^+ + \varphi^- \cdot \mathbf{n}^- \quad \text{if } \Gamma \in \mathcal{E}^I \qquad \llbracket \varphi \rrbracket = \varphi \cdot \mathbf{n} \quad \text{if } \Gamma \in \mathcal{E}^B \qquad (57a) $$

$$ \{\{\varphi\}\} = \frac{1}{2}(\varphi^+ + \varphi^-) \quad \text{if } \Gamma \in \mathcal{E}^I \qquad \{\{\varphi\}\} = \varphi \quad \text{if } \Gamma \in \mathcal{E}^B. \qquad (57b) $$
Finally, for vector functions, it is also useful to define a tensor jump as follows:

$$ \langle\langle\varphi\rangle\rangle = \varphi^+ \otimes \mathbf{n}^+ + \varphi^- \otimes \mathbf{n}^- \text{ if } \Gamma \in \mathcal{E}^I \qquad \langle\langle\varphi\rangle\rangle = \varphi \otimes \mathbf{n} \text{ if } \Gamma \in \mathcal{E}^B. \qquad (58) $$

We also introduce the following finite element spaces

$$ Q_r = \{v \in L^2(\Omega) : v|_{K \in \mathcal{Q}_r} \quad \forall K \in \mathcal{T}_H\} \qquad \mathbf{V}_r = [Q_r]^d, $$

where $\mathcal{Q}_r$ is the space of polynomials of degree $r$ in each coordinate direction. We then denote by $\varphi_i(\mathbf{x})$ the basis functions for the space $\mathbf{V}_r$ and by $\psi_i(\mathbf{x})$ the basis functions for the space $Q_r$, the finite element spaces chosen for the discretization of the velocity and of the pressure (as well as the density), respectively, so that

$$ \mathbf{u} \approx \sum_{j=1}^{|\mathcal{T}_h|(r+1)^d} u_j(t)\varphi_j(\mathbf{x}) \qquad p \approx \sum_{j=1}^{|\mathcal{T}_h|(r+1)^d} p_j(t)\psi_j(\mathbf{x}). $$

Here $|\mathcal{T}_h|$ denotes the number of elements of the computational mesh. Recall that $d$ denotes the space dimension. The number of degrees of freedom for scalar variable is indeed equal to $|\mathcal{T}_h|(r+1)^d$ [40]. The shape functions correspond to the products of Lagrange interpolation polynomials for the support points of $(r+1)$-order Gauss-Lobatto quadrature rule in each coordinate direction (Figure 1). In particular, we have grid points at the boundaries of the elements, where the solution can be discontinuous and this simplifies the evaluation of the integrals at the boundary itself [60]. Hence, for any given edge, the only shape functions with non-zero values are exactly those whose node points are located on that edge [60].

<figure><img src="image_2.png" /><figcaption>Figure 1: Example of two neighboring elements for a nodal DG formulation based on Lagrange polynomials. The nodes correspond to the support points of $(r+1)$-order Gauss-Lobatto quadrature rule (in the image $r=1$).</figcaption></figure>

Given these definitions, the weak formulation for the momentum equation at each stage (33) reads as follows [71, 74]:

$$ \mathbf{A}^{(l)}\mathbf{U}^{(l)} + \mathbf{B}^{(l)}\mathbf{P}^{(l)} = \mathbf{F}^{(n,l)}, \qquad (59) $$

with $\mathbf{U}^{(l)}$ denoting the vector of the degrees of freedom associated to the velocity field and $\mathbf{P}^{(l)}$
denoting the vector of the degrees of freedom associated to the pressure. Here we have set

$$A_{ij}^{(l)} = \sum_{K \in \mathcal{T}_h} \int_K \rho^{(l)} \varphi_j \cdot \varphi_i d\Omega \quad (60)$$

$$B_{ij}^{(l)} = \sum_{K \in \mathcal{T}_h} \int_K -a_{ll} \frac{\Delta t}{M^2} \nabla \cdot \varphi_i \Psi_j d\Omega + \sum_{\Gamma \in \mathcal{E}} \int_{\Gamma} a_{ll} \frac{\Delta t}{M^2} \{\{\Psi_j\}\} \llbracket \varphi_i \rrbracket d\Sigma \quad (61)$$

$$F_i^{(l)} = \sum_{K \in \mathcal{T}_h} \int_K \rho^n \mathbf{u}^n \cdot \varphi_i d\Omega \\ + \sum_{m=1}^{l-1} \sum_{K \in \mathcal{T}_h} \int_K \tilde{a}_{lm} \Delta t \left( \rho^{(m)} \mathbf{u}^{(m)} \otimes \mathbf{u}^{(m)} \right) : \nabla \varphi_i d\Omega + \sum_{m=1}^{l-1} \sum_{K \in \mathcal{T}_h} \int_K a_{lm} \frac{\Delta t}{M^2} p^{(m)} \nabla \cdot \varphi_i d\Omega \\ - \sum_{m=1}^{l-1} \sum_{\Gamma \in \mathcal{E}} \int_{\Gamma} \tilde{a}_{lm} \Delta t \left\{ \{ \rho^{(m)} \mathbf{u}^{(m)} \otimes \mathbf{u}^{(m)} \} \right\} : \langle\langle \varphi_i \rangle\rangle d\Sigma \\ - \sum_{m=1}^{l-1} \sum_{\Gamma \in \mathcal{E}} \int_{\Gamma} \tilde{a}_{lm} \Delta t \frac{\lambda^{(m)}}{2} \langle\langle \rho^{(m)} \mathbf{u}^{(m)} \rangle\rangle : \langle\langle \varphi_i \rangle\rangle d\Sigma \\ - \sum_{m=1}^{l-1} \sum_{\Gamma \in \mathcal{E}} \int_{\Gamma} a_{lm} \frac{\Delta t}{M^2} \left\{ \{ p^{(m)} \} \right\} \llbracket \varphi_i \rrbracket d\Sigma. \quad (62)$$

Following the discussion in [71, 74], one can notice that a centered flux is employed for the quantities defined implicitly, while an upwind-biased flux is adopted for the quantities computed explicitly. Moreover, following [1, 74], in order to obtain a numerical method effective for a wider range of Mach numbers, we take

$$\lambda^{(m)} = \max \left[ f \left( M_{\text{loc}}^{+, (m)} \right) \left( \left| \mathbf{u}^{+, (m)} \right| + \frac{1}{M} c^{+, (m)} \right) , f \left( M_{\text{loc}}^{-, (m)} \right) \left( \left| \mathbf{u}^{-, (m)} \right| + \frac{1}{M} c^{-, (m)} \right) \right], \quad (63)$$

with $M_{\text{loc}}^{\pm, (m)} = M_{c^{\pm, (m)}}^{\pm, (m)}$ and $f(M_{\text{loc}}) = \min(1, M_{\text{loc}})$. This choice corresponds to the convex combination between a centered flux and a Rusanov flux [80] as proposed in [1], so that, for $M_{\text{loc}} \ge 1$, we resort to a Rusanov flux, whereas for $M_{\text{loc}} \ll 1$, only the local fluid velocity is relevant for the numerical dissipation. Further considerations on the numerical flux will be discussed at the end of the Section. The numerical integration is based on the so-called over-integration or consistent integration, so as to guarantee exact integration. In particular, we employ $2r+1$ Gauss-Legendre quadrature points along each coordinate direction [73]. Analogously, the energy equation in (33) can be expressed as

$$\mathbf{C}^{(l)}\mathbf{U}^{(l)} + \mathbf{D}^{(l)}\mathbf{P}^{(l)} = \mathbf{G}^{(l)}. \quad (64)$$

For the sake of completeness, as well as to point out the contribution due to the novel strategy presented in Section 3.4 to handle a generic EOS, we report the expression of the components $\mathbf{C}^{(l)}$ and $\mathbf{D}^{(l)}$. The expression of $\mathbf{G}^{(l)}$ can be easily inferred from (33) and its definition entails that centered fluxes are employed for the quantities defined implicitly, while an upwind-biased flux is used for the quantities computed explicitly (see also [74]). Hence, we obtain

$$C_{ij}^{(l)} = \sum_{K \in \mathcal{T}_h} \int_K -a_{ll} \Delta t h^{(l)} \rho^{(l)} \varphi_j \cdot \nabla \Psi_i d\Omega + \sum_{\Gamma \in \mathcal{E}} \int_{\Gamma} a_{ll} \Delta t \left\{ \{h^{(l)} \rho^{(l)} \varphi_j\} \right\} \cdot \nabla \Psi_i d\Sigma \quad (65)$$

$$D_{ij}^{(l)} = \sum_{K \in \mathcal{T}_h} \int_K \frac{(\rho e)^{(l)}}{p^{(l)}} \Psi_j \Psi_i d\Omega \quad (66)$$

Formally, we can derive

$$\mathbf{U}^{(l)} = \left( \mathbf{A}^{(l)} \right)^{-1} \left( \mathbf{F}^{(l)} - \mathbf{B}^{(l)} \mathbf{P}^{(l)} \right), \quad (67)$$
so as to obtain

$$ \mathbf{D}^{(l)}\mathbf{P}^{(l)} + \mathbf{C}^{(l)}(\mathbf{A}^{(l)})^{-1}(\mathbf{F}^{(l)} - \mathbf{B}^{(l)}\mathbf{P}^{(l)}) = \mathbf{G}^{(l)}. \qquad (68) $$

The above system is then solved following the fixed point procedure described in [22, 36, 71]. More specifically, setting $\mathbf{U}^{(l,0)} = \mathbf{U}^{(l-1)}$, $\mathbf{P}^{(l,0)} = \mathbf{P}^{(l-1)}$, one solves for $\tilde{k} = 0, \dots, L$

$$ (\mathbf{D}^{(l,\tilde{k})} - \mathbf{C}^{(l,\tilde{k})}(\mathbf{A}^{(l)})^{-1}\mathbf{B}^{(l)})\mathbf{P}^{(l,\tilde{k}+1)} = \mathbf{G}^{(l,\tilde{k})} - \mathbf{C}^{(l,\tilde{k})}(\mathbf{A}^{(l)})^{-1}\mathbf{F}^{(l)} \quad (69) $$

and then updates the velocity solving

$$ \mathbf{A}^{(l)}\mathbf{U}^{(l,\tilde{k})} = \mathbf{F}^{(l)} - \mathbf{B}^{(l)}\mathbf{P}^{(l,\tilde{k}+1)}. \quad (70) $$

The algebraic formulation associated to (54) is obtained substituting the degrees of freedom of the velocity into the algebraic formulation of the energy equation. Relations (67) and (68) can be therefore employed to achieve this goal. For the sake of completeness, we report the new definitions of $\mathbf{C}^{(l)}$ and $\mathbf{D}^{(l)}$, while analogous modifications apply to the other variables. Hence, $\mathbf{C}^{(l)}$ and $\mathbf{D}^{(l)}$ now read as follows:

$$ \begin{aligned} C_{ij}^{(l)} = & \sum_{K \in \mathcal{T}_h} \int_K -a_{ll} \Delta t h_E^{(l)} \rho_I^{(l)} \varphi_j \cdot [[\Psi_i]] d\Omega \\ & + \sum_{\Gamma \in \mathcal{E}} \int_{\Gamma} a_{ll} \Delta t \{\{h_E^{(l)} \rho_I^{(l)} \varphi_j\}\} \cdot \nabla \Psi_i d\Sigma + \sum_{K \in \mathcal{T}_h} \int_K \frac{M^2}{2} \mathbf{u}_E^{(l)} \cdot \rho_I^{(l)} \varphi_j \Psi_i d\Omega \end{aligned} \quad (71) $$

$$ D_{ij}^{(l)} = \sum_{K \in \mathcal{T}_h} \int_K \frac{(\rho e)_E^{(l)}}{p_E^{(l)}} \Psi_j \Psi_i d\Omega. \quad (72) $$

A matrix-free approach is employed [3], meaning that no global sparse matrix is built and only the action of the linear operators on a vector is actually implemented. Matrices $\mathbf{A}^{(l)}$ and $\mathbf{D}^{(l)}$ are symmetric and positive definite, while the matrix $\mathbf{C}^{(l)}(\mathbf{A}^{(l)})^{-1}\mathbf{B}^{(l)}$ is not symmetric. We point out that if one directly discretizes (35) and (50), as done, e.g., in [19, 78], a symmetric positive definite linear system can be obtained. However, this approach implies the direct numerical solution of an elliptic equation and the discretization of a second order operator that, in the framework of a DG method, would require, e.g., the use of the Symmetric Interior Penalty method [4]. The use of a Schur complement type technique, as the one described in this work, allows one to employ only the standard numerical fluxes of hyperbolic problems (Rusanov and upwind-biased in this work), without defining and setting penalization constants typical of the aforementioned numerical strategy for elliptic equations. A comparison between the approach employed in this work and the direct solution of the Helmholtz-type equations (35) and (50) will be matter of future work. In view of these considerations, a preconditioned conjugate gradient method with a geometric multigrid preconditioner is employed to solve the symmetric positive definite linear systems. The GMRES solver with a Jacobi preconditioner is employed for the solution of the non-symmetric linear systems. In future developments, we aim to implement and employ multigrid preconditioners also for the non-symmetric linear systems in the matrix-free framework so as to further improve the performance of the solver.

The DG method naturally allows for high-order accuracy without the use of reconstructions which involve large stencils. However, as discussed in [53], its accuracy in the very low Mach regime depends on the numerical flux and on the shape of the elements. More specifically, a simplicial mesh is needed to establish low Mach accuracy. The lack of accuracy in the very low Mach limit can lead to a numerical scheme which is not convergent with a finite volume scheme, while an order reduction is observed in the case of the Discontinuous Galerkin method [53]. A low Mach fix for the Euler equations resolved with the finite volume method on Cartesian grids was proposed in [10]. Moreover, it is known that an upwind scheme fails to solve very subsonic
flows [46]. However, for moderate low Mach numbers, i.e. $M > 10^{-4}$, the convex combination (63) leads to a correct scaling of the pressure fluctuations [1, 74] and the low Mach number inaccuracy is counterbalanced by the high-order nature of the numerical scheme [53, 74].

A detailed analysis of the spatial discretization and the development of possible remedies for very subsonic flows goes beyond the scope of the present work and will subject of future developments (see also Sections 4.1 and 4.2). Since the main focus of this work is the comparison of the time discretization methods presented in Section 3.2 and 3.3, we believe that considering a minimum Mach number around $10^{-3} - 10^{-4}$ allows us to perform this analysis in the low Mach limit and compares well with the minimum Mach number chosen to analyze asymptotic-preserving schemes in the literature, see, e.g., [1, 20, 18, 32, 67, 85].

# 4 Numerical results

The numerical methods outlined in Section 3 are now validated in a number of relevant benchmarks. The implementation is carried out in the framework of the deal.II library [3, 9], that is a C++ open-source software supporting the creation of finite element codes. Several libraries based on deal.II have been developed in the last years [2, 45, 59]. All the simulations are performed in double precision. The employed time discretization schemes are reported in Appendix A. Discrete parameter choices are associated to two Courant numbers, one based on the speed of sound denoted by C, the so-called acoustic Courant number, and one based on the local velocity of the flow, the so-called advective Courant number, denoted by $C_u$:

$$C = \frac{1}{M} rc \frac{\Delta t}{\mathcal{H}} \sqrt{d} \quad C_u = ru \frac{\Delta t}{\mathcal{H}} \sqrt{d}. \qquad (73)$$

Here, $\mathcal{H} = \min\{\text{diam}(K) | K \in \mathcal{T}_h\}$, r is the polynomial degree employed for the spatial discretization, c is the speed of sound, and u is the magnitude of the flow velocity. Recall that d denotes the space dimension. For what concerns the tests with the ideal gas law, $\gamma = 1.4$ is employed in (4). Finally, following [71], the fixed point loop (37) is stopped at the iteration $\tilde{k}$ for which the maximum relative difference for the pressure is below a tolerance $\eta$, namely

$$\frac{\|\xi(\tilde{k}) - \xi(\tilde{k}-1)\|_{\infty}}{\|\xi(\tilde{k})\|_{\infty}} < \eta. \qquad (74)$$

## 4.1 Taylor-Green vortex

As a first benchmark to verify the scaling properties of the numerical methods with respect to the Mach number M, we consider the Taylor-Green vortex [29, 90], that represents an exact steady solution of the incompressible Euler equations. The initial condition in non-dimensional variables reads as follows:

$$\rho(\mathbf{x}, 0) = 1 \quad \mathbf{u}(\mathbf{x}, 0) = \begin{pmatrix} \sin(x) \cos(y) \\ -\cos(x) \sin(y) \end{pmatrix} \quad p = 1 + \frac{1}{4} M^2 (\cos(2x) + \cos(2y)). \qquad (75)$$

The computational domain is $\Omega = (0, 2\pi)^2$ endowed with periodic boundary conditions. The time step is such that the maximum advective Courant number is $C_u \approx 0.095$, while the maximum acoustic Courant number is $C \approx 0.11/M$. We employ $\eta = 10^{-10}$.

First, we consider the IMEX-RK(3,3,3) scheme of type II (Table 21) in combination with polynomial degree $r = 2$ and $N_{el} = 60$ elements along each coordinate direction. We employ the IMEX-DG method. One can easily notice that pressure fluctuations scale as $O(M^2)$ up to $M = 10^{-4}$, as expected, whereas the density fluctuations scale as $O(M^2)$ up to $M = 10^{-3}$ and a degradation is experienced at $M = 10^{-4}$. This degradation is likely related to the low Mach inaccuracy of the DG method on quadrilateral cells discussed in Section 3.5. For what
concerns the divergence of the velocity field, it scales as $O(H^2)$ (see the convergence analysis in Section 4.2), but it does not vanish as $M \to 0$ (Table 1). As discussed in [74], since the initial velocity field is solenoidal and the vortex is stationary, a quadratic convergence with respect to $M$ is expected for the divergence of the velocity field. This result is dependent on the spatial discretization. Indeed, since our method employs a standard nodal DG method, the divergence-free property is not imposed pointwise and the error associated to $\nabla \cdot \mathbf{u}$ is basically constant in time and it is therefore related to the interpolation of the initial datum into the employed finite element space (Figure 2). A quadratic convergence with respect to $M$ was recently obtained in [90] and preliminary results in our framework suggest that the use of Raviart-Thomas finite elements [5] for the velocity field improves the scaling properties as $M \to 0$. As already discussed at the end of Section 3.5, a more detailed analysis of the spatial discretization is currently under investigation, while the primary goal of the present work is to perform a quantitative comparison between two different IMEX-RK approaches for the Euler equations. However, further considerations about the spatial discretization will be added at the end of the upcoming Section 4.2.

Next, we consider the SI-IMEX-DG method. A stable solution is obtained up to $M = 10^{-3}$ and the density fluctuations scale as $O(M^2)$ only up to $M = 10^{-2}$ (Table 2). This degradation is again likely mainly related to an early manifestation of the low Mach inaccuracy of the DG method on quadrilateral cells. However, we point out that some issues in the low Mach regime employing schemes of type II for the semi-implicit time discretization were already experienced in [17]

<table><caption>Table 1: Mach number scaling of the density and pressure fluctuations and of the divergence of the velocity field for the Taylor-Green vortex test case. Here, and in the following Tables, $\delta\rho = \rho - 1$ and $\delta p = p - 1$. The results are obtained using the IMEX method with the IMEX-RK(3,3,3) scheme of type II in Table 21 together with polynomial degree $r = 2$ and $N_{el} = 60$.</caption><thead><tr><th>M</th><th>L<sup>2</sup> norm<br>∇·u</th><th>Rate ∇·u</th><th>L<sup>2</sup> norm δρ</th><th>Rate δρ</th><th>L<sup>2</sup> norm δp</th><th>Rate δp</th></tr></thead><tbody><tr><td>10<sup>-1</sup></td><td>5.76e-3</td><td></td><td>1.38e-3</td><td></td><td>1.57e-2</td><td></td></tr><tr><td>10<sup>-2</sup></td><td>1.82e-3</td><td>0.5</td><td>1.57e-5</td><td>1.9</td><td>1.57e-4</td><td>2.0</td></tr><tr><td>10<sup>-3</sup></td><td>1.82e-3</td><td>—</td><td>1.58e-7</td><td>2.0</td><td>1.57e-6</td><td>2.0</td></tr><tr><td>10<sup>-4</sup></td><td>1.82e-3</td><td>—</td><td>1.24e-8</td><td>1.1</td><td>1.62e-8</td><td>2.0</td></tr></tbody></table>

<table><caption>Table 2: Mach number scaling of the density and pressure fluctuations and of the divergence of the velocity field for the Taylor-Green vortex test case. The results are obtained using the SI-IMEX method with the IMEX-RK(3,3,3) scheme of type II in Table 21 together with polynomial degree $r = 2$ and $N_{el} = 60$.</caption><thead><tr><th>M</th><th>L<sup>2</sup> norm<br>∇·u</th><th>Rate ∇·u</th><th>L<sup>2</sup> norm δρ</th><th>Rate δρ</th><th>L<sup>2</sup> norm δp</th><th>Rate δp</th></tr></thead><tbody><tr><td>10<sup>-1</sup></td><td>5.76e-3</td><td></td><td>1.38e-3</td><td></td><td>1.57e-2</td><td></td></tr><tr><td>10<sup>-2</sup></td><td>1.82e-3</td><td>0.5</td><td>1.61e-5</td><td>1.9</td><td>1.57e-4</td><td>2.0</td></tr><tr><td>10<sup>-3</sup></td><td>1.82e-3</td><td>—</td><td>4.28e-6</td><td>0.6</td><td>1.57e-6</td><td>2.0</td></tr></tbody></table>

Next, we consider the IMEX-RK(4,3,3) scheme of type I (Table 20). For what concerns the IMEX method, analogous results with respect to the scheme of type II are obtained up to $M = 10^{-3}$ with an improvement of the scaling of the density fluctuations at $M = 10^{-4}$ (Table
<figure><img src="image_3.png" /><figcaption>Figure 2: Taylor-Green vortex test case, time evolution of the divergence of the velocity at $M = 10^{-3}$. The results are obtained using the IMEX method with the IMEX-RK(3,3,3) scheme of type II in Table 21 together with polynomial degree $r = 2$ and $N_{el} = 60$.</figcaption></figure>

3). For what concerns the SI-IMEX method, a stable solution is established at $M = 10^{-4}$ with a stagnation of the density fluctuations (Table 4). These results confirm the superior stability of schemes of type I with respect to schemes of type II for low Mach numbers flows already experienced in [17]. Moreover, we can infer that the low Mach accuracy is influenced also by the time discretization strategy and by the time discretization scheme.

<table><caption>Table 3: Mach number scaling of the density and pressure fluctuations and of the divergence of the velocity field for the Taylor-Green vortex test case. The results are obtained using the IMEX method with the IMEX-RK(4,3,3) scheme of type I in Table 20 together with polynomial degree $r = 2$ and $N_{el} = 60$.</caption><thead><tr><th>M</th><th>L<sup>2</sup> norm<br>∇·u</th><th>Rate ∇·u</th><th>L<sup>2</sup> norm δρ</th><th>Rate δρ</th><th>L<sup>2</sup> norm δp</th><th>Rate δp</th></tr></thead><tbody><tr><td>10<sup>-1</sup></td><td>5.76e-3</td><td></td><td>1.38e-3</td><td></td><td>1.57e-2</td><td></td></tr><tr><td>10<sup>-2</sup></td><td>1.82e-3</td><td>0.5</td><td>1.58e-5</td><td>1.9</td><td>1.57e-4</td><td>2.0</td></tr><tr><td>10<sup>-3</sup></td><td>1.82e-3</td><td>—</td><td>1.58e-7</td><td>2.0</td><td>1.57e-6</td><td>2.0</td></tr><tr><td>10<sup>-4</sup></td><td>1.82e-3</td><td>—</td><td>4.66e-9</td><td>1.5</td><td>1.63e-8</td><td>2.0</td></tr></tbody></table>

## 4.2 Traveling vortex at low Mach

Next, we consider for an ideal gas a two-dimensional traveling vortex inspired by the inviscid isentropic vortex studied, e.g., in [71, 84, 91]. For this test, a time-dependent analytic solution is available and the convergence properties of a numerical scheme can be therefore assessed.
<table><caption>Table 4: Mach number scaling of the density and pressure fluctuations and of the divergence of the velocity field for the Taylor-Green vortex test case. The results are obtained using the SI-IMEX method with the IMEX-RK(4,3,3) scheme of type I in Table 20 together with polynomial degree $r = 2$ and $N_{el} = 60$.</caption><thead><tr><th>M</th><th>L<sup>2</sup> norm ∇·u</th><th>Rate ∇·u</th><th>L<sup>2</sup> norm δρ</th><th>Rate δρ</th><th>L<sup>2</sup> norm δp</th><th>Rate δp</th></tr></thead><tbody><tr><td>10<sup>-1</sup></td><td>5.76e-3</td><td></td><td>1.38e-3</td><td></td><td>1.57e-2</td><td></td></tr><tr><td>10<sup>-2</sup></td><td>1.82e-3</td><td>0.5</td><td>1.58e-5</td><td>1.9</td><td>1.57e-4</td><td>2.0</td></tr><tr><td>10<sup>-3</sup></td><td>1.82e-3</td><td>—</td><td>1.69e-6</td><td>1.0</td><td>1.57e-6</td><td>2.0</td></tr><tr><td>10<sup>-4</sup></td><td>1.82e-3</td><td>—</td><td>1.68e-6</td><td>—</td><td>1.63e-8</td><td>2.0</td></tr></tbody></table>

The exact solution is indeed a propagation of the initial condition at the background velocity

$$ \rho(\mathbf{x}, t) = \rho(\mathbf{x} - \mathbf{u}_{\infty}t, 0) \quad \mathbf{u}(\mathbf{x}, t) = \mathbf{u}(\mathbf{x} - \mathbf{u}_{\infty}t, 0) \quad p(\mathbf{x}, t) = p(\mathbf{x} - \mathbf{u}_{\infty}t, 0). $$

Notice that a different version of this test case, which allows for a steady solution, was employed in [18, 19, 74]. Following [91], in order to emphasize the role of the Mach number $M$, we define the perturbation as

$$ \delta T = \frac{1 - \gamma}{8\gamma\pi^2} M^2 \beta^2 e^{1-\tilde{r}^2}, \qquad (76) $$

with

$$ \tilde{r}^2 = \frac{(x - x_0)^2 + (y - y_0)^2}{r_0^2} $$

denoting the scaled radial coordinate and $\beta$ being the vortex strength. We set $\beta = 10$, $r_0 = 1$, and

$$ \rho(\mathbf{x}, 0) = (1 + \delta T)^{\frac{1}{\gamma-1}} \quad p(\mathbf{x}, 0) = 1 + M^2 (1 + \delta T)^{\frac{\gamma}{\gamma-1}} = 1 + M^2 \rho^\gamma. \qquad (77) $$

For what concerns the velocity, we define its perturbation as

$$ \delta \mathbf{u} = \beta M \begin{pmatrix} -(y - y_0) \\ (x - x_0) \end{pmatrix} \frac{e^{\frac{1}{2}(1 - \tilde{r}^2)}}{2\pi} \qquad (78) $$

Finally, we set $\mathbf{u}_{\infty} = [1, 1]^\top$, whereas the final time is $T_f = 3$. To avoid problems related to the definition of boundary conditions, we choose a sufficiently large domain $\Omega = (-10, 10)^2$ and periodic boundary conditions. The behavior of the numerical methods is investigated at different Mach numbers. More specifically, we consider $M \in [10^{-4}, 10^{-1}]$. The tolerance in (74) is set to $\eta = 10^{-10}$.

First, we analyze the results obtained employing the IMEX-RK(2,2,2) scheme (Table 19) with polynomial degree $r = 1$ for the space discretization. The time step is chosen in such a way that the maximum advective Courant number is $C_u \approx 0.16$ and the maximum acoustic Courant number is $C \approx 0.12/M$. One can easily notice that IMEX and SI-IMEX time discretization methods yield the same level of accuracy (Tables 5-6). A degradation is experienced at $M = 10^{-2}$ and analogous results are obtained for lower values of $M$. Contour plots of the pressure perturbation show that the shape of the vortex is not preserved (Figure 3). This is likely related to well known issues of collocated finite element type discretization on quadrilateral meshes in the low Mach regime [53, 79] (see also the discussion at the end of Section 3.5). A simple workaround consists of increasing the polynomial degree of the finite element space employed for the discretization of the velocity field. We also refer to [52], where a space enrichment for the velocity has been used in the framework of a finite volume scheme. One can easily notice that, if we take $r = 2$ for the finite element space of the velocity, the shape of the vortex is preserved (Figure 4) and the expected second order convergence is established (Table 7).
<table><caption>Table 5: Convergence analysis for the traveling vortex test case using the IMEX method with the IMEX-RK(2,2,2) scheme (Table 19) and polynomial degree $r = 1$. Here, and in the following Tables, $N_{el}$ denotes the number of elements along each direction.</caption><thead><tr><th rowspan="2">N<sub>el</sub></th><th colspan="2">M = 10<sup>-1</sup></th><th colspan="2">M = 10<sup>-2</sup></th></tr><tr><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate δu</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate δu</th></tr></thead><tbody><tr><td>40</td><td>1.28e-1</td><td></td><td>1.33e-1</td><td></td></tr><tr><td>80</td><td>3.99e-2</td><td>1.7</td><td>5.18e-2</td><td>1.3</td></tr><tr><td>160</td><td>1.02e-2</td><td>2.1</td><td>2.27e-2</td><td><b>1.2</b></td></tr><tr><td>320</td><td>2.03e-3</td><td><b>2.1</b></td><td></td><td></td></tr></tbody></table>

<table><caption>Table 6: Convergence analysis for the traveling vortex test case using the SI-IMEX method with the IMEX-RK(2,2,2) scheme (Table 19) and polynomial degree $r = 1$.</caption><thead><tr><th rowspan="2">N<sub>el</sub></th><th colspan="2">M = 10<sup>-1</sup></th><th colspan="2">M = 10<sup>-2</sup></th></tr><tr><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate δu</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate δu</th></tr></thead><tbody><tr><td>40</td><td>1.28e-1</td><td></td><td>1.32e-1</td><td></td></tr><tr><td>80</td><td>3.99e-2</td><td>1.7</td><td>5.18e-2</td><td>1.3</td></tr><tr><td>160</td><td>1.02e-2</td><td>2.0</td><td>2.27e-2</td><td><b>1.2</b></td></tr><tr><td>320</td><td>2.03e-3</td><td><b>2.3</b></td><td></td><td></td></tr></tbody></table>

<table><caption>Table 7: Convergence analysis for the traveling vortex test case using the IMEX method with the IMEX-RK(2,2,2) scheme (Table 19), polynomial degree $r = 2$ for the velocity and $r = 1$ for the remaining variables.</caption><thead><tr><th rowspan="2">N<sub>el</sub></th><th colspan="2">M = 10<sup>-2</sup></th></tr><tr><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate δu</th></tr></thead><tbody><tr><td>20</td><td>2.48e-2</td><td></td></tr><tr><td>40</td><td>3.08e-3</td><td>3.0</td></tr><tr><td>80</td><td>4.86e-4</td><td>2.7</td></tr><tr><td>160</td><td>9.30e-5</td><td><b>2.4</b></td></tr></tbody></table>

Next, we focus on third order time discretization schemes (Tables 20-21) in combination with polynomial degree $r = 2$. The time step is such that the maximum advective Courant number is $C_u \approx 0.094$, while the maximum acoustic Courant number is $C \approx 0.07/M$. All the schemes provide a similar level of accuracy and the expected third order convergence rate is established (Tables 8-11) up to $M = 10^{-3}$, except for the SI-IMEX method in combination with the scheme of type II at $M = 10^{-3}$. As already remarked in Section 4.1, some issues in the low Mach regime employing schemes of type II for the semi-implicit time discretization were already experienced in [17]. Moreover, for $M = 10^{-4}$, a stable solution for the SI-IMEX method is obtained only employing the scheme of type I. This further confirms the superior stability of schemes of type I with respect to schemes of type II for low Mach numbers when using the SI-IMEX method.

A saturation of the error is experienced at $M = 10^{-4}$. Analogous results with a saturation of the error at $M = 10^{-4}$ were obtained in [91]. Since in most low Mach number applications the
<figure><img src="image_4.png" /><figcaption>Figure 3: Traveling vortex test case with the IMEX-RK(2,2,2) scheme (Table 19) and polynomial degree $r = 1$, contour plots of the pressure perturbation $p - (1 - M^2)$ at $M = 10^{-2}$ with $N_{el} = 160$. Left: initial field. Right: field at final time $t = T_f = 3$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Traveling vortex test case with the IMEX-RK(2,2,2) scheme (Table 19), $r = 2$ for the velocity field and $r = 1$ for the remaining variables, contour plots of the pressure perturbation $p - (1 + M^2)$ at $M = 10^{-2}$ with $N_{el} = 160$. Left: initial field. Right: field at final time $t = T_f = 3$.</figcaption></figure>

main interest lies on the velocity field rather than on the acoustics, we notice, similarly to [91], that no visible spurious effect arises in the velocity field in spite of the order reduction (Figure 5). As already discussed in Section 3.5, the DG method on quadrilateral cells becomes low Mach inaccurate for $M < 10^{-4}$ and therefore this order reduction is likely related to a manifestation of this phenomenon. Moreover, round-off errors are dominant in this configuration [91] and the use of quadruple precision is crucial to reach really very low Mach numbers [21]. However, it is worth to notice that in the case of the steady vortex, a correct scaling was established at $M = 10^{-4}$ using the third order method of type II in [74]. Hence, the saturation of the error could be also related to the well-known order reduction phenomenon experienced for stiff ODE's problems [88]. In future work, as already mentioned, we aim to analyze the behavior employing a spatial discretization based on compatible finite elements, as recently done in [90], so as to guarantee that the initial condition of the velocity field is set up as the discrete derivative of a vector potential, and to use simplicial meshes [53, 90] or Voronoi meshes [18], which have been
shown to be low Mach accurate for a steady vortex.

<table><caption>Table 8: Convergence analysis for the traveling vortex test case using the IMEX method with the IMEX-RK(3,3,3) scheme of type II (Table 21) and polynomial degree $r = 2$.</caption><thead><tr><th rowspan="2">N<sub>el</sub></th><th colspan="2">M = 10<sup>-1</sup></th><th colspan="2">M = 10<sup>-2</sup></th><th colspan="2">M = 10<sup>-3</sup></th><th colspan="2">M = 10<sup>-4</sup></th></tr><tr><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th></tr></thead><tbody><tr><td>20</td><td>3.12e-2</td><td></td><td>3.55e-2</td><td></td><td>3.55e-2</td><td></td><td>3.59e-2</td><td></td></tr><tr><td>40</td><td>2.40e-3</td><td>3.7</td><td>2.32e-3</td><td>3.9</td><td>2.32e-3</td><td>3.9</td><td>5.51e-3</td><td>2.7</td></tr><tr><td>80</td><td>3.17e-4</td><td>2.9</td><td>2.92e-4</td><td>3.0</td><td>2.90e-4</td><td>3.0</td><td>5.71e-3</td><td>-</td></tr><tr><td>160</td><td>4.13e-5</td><td><b>2.9</b></td><td>3.73e-5</td><td><b>3.0</b></td><td>4.61e-5</td><td><b>2.7</b></td><td></td><td></td></tr></tbody></table>

<table><caption>Table 9: Convergence analysis for the traveling vortex test case using the SI-IMEX method with the IMEX-RK(3,3,3) scheme of type II (Table 21) and polynomial degree $r = 2$.</caption><thead><tr><th rowspan="2">N<sub>el</sub></th><th colspan="2">M = 10<sup>-1</sup></th><th colspan="2">M = 10<sup>-2</sup></th><th colspan="2">M = 10<sup>-3</sup></th></tr><tr><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th></tr></thead><tbody><tr><td>20</td><td>3.12e-2</td><td></td><td>3.56e-2</td><td></td><td>3.57e-2</td><td></td></tr><tr><td>40</td><td>2.40e-3</td><td>3.7</td><td>2.33e-3</td><td>3.9</td><td>2.33e-3</td><td>3.9</td></tr><tr><td>80</td><td>3.17e-4</td><td>2.9</td><td>2.92e-4</td><td>3.0</td><td>3.85e-4</td><td>2.6</td></tr><tr><td>160</td><td>4.13e-5</td><td><b>2.9</b></td><td>3.73e-5</td><td><b>3.0</b></td><td>2.42e-4</td><td>0.7</td></tr></tbody></table>

<table><caption>Table 10: Convergence analysis for the traveling vortex test case using the IMEX method with the IMEX-RK(4,3,3) scheme of type I (Table 20) and polynomial degree $r = 2$.</caption><thead><tr><th rowspan="2">N<sub>el</sub></th><th colspan="2">M = 10<sup>-1</sup></th><th colspan="2">M = 10<sup>-2</sup></th><th colspan="2">M = 10<sup>-3</sup></th><th colspan="2">M = 10<sup>-4</sup></th></tr><tr><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error u</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th></tr></thead><tbody><tr><td>20</td><td>3.12e-2</td><td></td><td>3.55e-2</td><td></td><td>3.55e-2</td><td></td><td>3.63e-2</td><td></td></tr><tr><td>40</td><td>2.40e-3</td><td>3.7</td><td>2.32e-3</td><td>3.9</td><td>2.31e-3</td><td>3.9</td><td>3.57e-3</td><td>3.3</td></tr><tr><td>80</td><td>3.17e-4</td><td>2.9</td><td>2.92e-4</td><td>3.0</td><td>2.90e-4</td><td>3.0</td><td>4.50e-3</td><td>-</td></tr><tr><td>160</td><td>4.13e-5</td><td><b>2.9</b></td><td>3.73e-5</td><td><b>3.0</b></td><td>4.61e-5</td><td><b>2.7</b></td><td></td><td></td></tr></tbody></table>

The IMEX method described in Section 3.2 allows the use of a larger time step ensuring the same level of accuracy, because of the superior stability provided by the fixed point loop [36]. In particular, the time step can be doubled, yielding a maximum advective Courant number $C_u \approx 0.19$. In spite of the use of a smaller time step, the SI-IMEX method and, more specifically, the SI-IMEX time discretization using the scheme of type II, provides in general better computational performance (Figure 6). As already discussed, the scheme of type I is
<table><caption>Table 11: Convergence analysis for the traveling vortex test case using the SI-IMEX method with IMEX-RK(4,3,3) scheme of type I (Table 20) and polynomial degree $r = 2$</caption><thead><tr><th rowspan="2">N<sub>el</sub></th><th colspan="2">M = 10<sup>-1</sup></th><th colspan="2">M = 10<sup>-2</sup></th><th colspan="2">M = 10<sup>-3</sup></th><th colspan="2">M = 10<sup>-4</sup></th></tr><tr><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th></tr></thead><tbody><tr><td>20</td><td>3.12e-2</td><td></td><td>3.55e-2</td><td></td><td>3.56e-2</td><td></td><td>4.70e-2</td><td></td></tr><tr><td>40</td><td>2.40e-3</td><td>3.7</td><td>2.32e-3</td><td>3.9</td><td>2.32e-3</td><td>3.9</td><td>2.12e-2</td><td>1.1</td></tr><tr><td>80</td><td>3.17e-4</td><td>2.9</td><td>2.92e-4</td><td>3.0</td><td>2.90e-4</td><td>3.0</td><td>1.01e-2</td><td>1.1</td></tr><tr><td>160</td><td>4.13e-5</td><td><b>2.9</b></td><td>3.73e-5</td><td><b>2.9</b></td><td>4.69e-5</td><td><b>2.7</b></td><td></td><td></td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 5: Traveling vortex test case, contour plot of the velocity perturbation at $M = 10^{-4}$. Left: initial field. Right: field at final time $t = T_f = 3$. The results are obtained using the IMEX method with the IMEX-RK(3,3,3) scheme (Table 21) and polynomial degree $r = 2$ with $N_{el} = 40$.</figcaption></figure>

more robust for $M = 10^{-3}$ and it becomes also more efficient as the spatial resolution increases. In the low Mach regime, differences in terms of computational cost between the IMEX method and the SI-IMEX method are significantly smaller: the efficiency gain reduces to around 5% for $M = 10^{-3}$, meaning that the two methods are essentially equivalent in this respect. We will further discuss this point in Section 4.2.1.

Next, we employ the fourth order time discretization schemes of type ARS (Table 22) and of type II (Table 23) in combination with polynomial degree $r = 3$. The time step is such that the maximum advective Courant number is $C_u \approx 0.05$ and the maximum acoustic Courant number is $C \approx 0.035/M$. The SI-IMEX method requires a smaller time step for $M \le 10^{-2}$. More specifically, an advective Courant number $C_u \approx 0.025$ is required to achieve a stable solution. One can easily notice that the expected convergence rates are established up to $M = 10^{-3}$ (Tables 12-13), except for the IMEX method at $M = 10^{-2}$, for which an order reduction is experienced. Analogous considerations to those reported for the third order discretization schemes for $M \le 10^{-4}$ are valid. In particular, the correct scaling at $M = 10^{-4}$ for the steady vortex using the IMEX method with the scheme of type ARS was obtained in [74], so that the experienced order reduction is likely dependent on both the low Mach inaccuracy of the DG method on quadrilateral cells and the order reduction phenomenon typical of stiff problems [88].

For what concerns the time discretization method of type II (Table 23), the expected conver-
<figure><img src="image_7.png" /><figcaption>Figure 6: Traveling vortex test case, comparison of wall-clock times employing third order time discretization schemes and polynomial degree $r = 2$. Top-left: $M = 10^{-1}$. Top-right: $M = 10^{-2}$. Bottom: $M = 10^{-3}$. An advective Courant number $C_u \approx 0.19$ is employed for the IMEX method. The solid black lines denote the results obtained with the IMEX method with the scheme of type II (Table 21), the solid blue lines show the results with the SI-IMEX method with the scheme of type II (Table 21), the dashed red lines report the results achieved with the IMEX method with the scheme of type I (Table 20), whereas the dashed green lines represent the results established with the SI-IMEX method with the scheme of type I (Table 20).</figcaption></figure>

gence rates are established for the IMEX method and the accuracy is preserved up to $M = 10^{-3}$ (Table 14). On the contrary, the SI-IMEX method shows the expected convergence rates for $M \ge 10^{-2}$, whereas severe issues start from $M = 10^{-3}$ (Table 15), unlike the scheme of type ARS. As already mentioned, schemes of type II show some issues in the low Mach regime for the semi-implicit time discretization [17]. Moreover, following also the results in Section 4.1, we can further infer that the activation of spurious modes (hence the low Mach accuracy) is influenced by the time discretization strategy and by the time discretization scheme.

### 4.2.1 Investigation of fixed point iterations

In this Section, we analyze the impact of the value of the tolerance $\eta$ in the fixed point loop on the overall performance of the IMEX method depicted in Section 3.2. For the sake of brevity, we focus on $M = 10^{-1}$ and on $M = 10^{-3}$. We consider the third order scheme of type II (Table 21), but analogous considerations hold for the other schemes. For moderate values of the Mach number, a sufficiently low value of the tolerance $\eta$ has to be chosen in order to achieve full convergence and the number of fixed point iterations depends on the value of $\eta$ (Table 16). For what concerns low values of the Mach number, one fixed point iteration is sufficient to achieve full convergence, independently of $\eta$ (Table 16). This explains why the IMEX method and
<table><caption>Table 12: Convergence analysis for the traveling vortex test case using the IMEX method with the IMEX-RK(5,5,4) scheme of type ARS (Table 22) and polynomial degree $r = 3$.</caption><thead><tr><th rowspan="2">N<sub>el</sub></th><th colspan="2">M = 10<sup>-1</sup></th><th colspan="2">M = 10<sup>-2</sup></th><th colspan="2">M = 10<sup>-3</sup></th><th colspan="2">M = 10<sup>-4</sup></th></tr><tr><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th></tr></thead><tbody><tr><td>10</td><td>4.85e-2</td><td></td><td>4.88e-2</td><td></td><td>4.88e-2</td><td></td><td>9.15e-2</td><td></td></tr><tr><td>20</td><td>2.88e-3</td><td>4.1</td><td>2.73e-3</td><td>4.2</td><td>2.71e-3</td><td>4.2</td><td>6.37e-2</td><td>0.5</td></tr><tr><td>40</td><td>1.23e-4</td><td>4.5</td><td>1.30e-4</td><td>4.4</td><td>1.31e-4</td><td>4.4</td><td>2.23e-2</td><td>1.5</td></tr><tr><td>80</td><td>8.47e-6</td><td><b>3.9</b></td><td>6.41e-5</td><td>1.0</td><td>1.22e-5</td><td><b>3.4</b></td><td></td><td></td></tr></tbody></table>

<table><caption>Table 13: Convergence analysis for the traveling vortex test case using the SI-IMEX method with the IMEX-RK(5,5,4) scheme of type ARS (Table 22) and polynomial degree $r = 3$. Note that for $M = 10^{-1}$ the advective Courant number is $C_u \approx 0.05$, while for $M \le 10^{-2}$ the advective Courant number is $C_u \approx 0.025$ to achieve a stable solution.</caption><thead><tr><th rowspan="2">N<sub>el</sub></th><th colspan="2">M = 10<sup>-1</sup></th><th colspan="2">M = 10<sup>-2</sup></th><th colspan="2">M = 10<sup>-3</sup></th><th colspan="2">M = 10<sup>-4</sup></th></tr><tr><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate u</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th></tr></thead><tbody><tr><td>10</td><td>4.85e-2</td><td></td><td>4.89e-2</td><td></td><td>4.88e-2</td><td></td><td>7.38e-2</td><td></td></tr><tr><td>20</td><td>2.88e-3</td><td>4.1</td><td>2.74e-3</td><td>4.2</td><td>2.71e-3</td><td>4.2</td><td>6.52e-2</td><td>0.2</td></tr><tr><td>40</td><td>1.23e-4</td><td>4.5</td><td>1.30e-4</td><td>4.4</td><td>1.31e-4</td><td>4.4</td><td>6.27e-2</td><td>-</td></tr><tr><td>80</td><td>8.47e-6</td><td><b>3.9</b></td><td>1.10e-5</td><td><b>3.6</b></td><td>1.24e-5</td><td><b>3.4</b></td><td></td><td></td></tr></tbody></table>

<table><caption>Table 14: Convergence analysis for the traveling vortex test case using the IMEX method with the IMEX-RK(6,6,4) scheme of type II (Table 23) and polynomial degree $r = 3$.</caption><thead><tr><th rowspan="2">N<sub>el</sub></th><th colspan="2">M = 10<sup>-1</sup></th><th colspan="2">M = 10<sup>-2</sup></th><th colspan="2">M = 10<sup>-3</sup></th><th colspan="2">M = 10<sup>-4</sup></th></tr><tr><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate</th></tr></thead><tbody><tr><td>10</td><td>4.85e-2</td><td></td><td>4.88e-2</td><td></td><td>4.88e-2</td><td></td><td>5.90e-2</td><td></td></tr><tr><td>20</td><td>2.88e-3</td><td>4.1</td><td>2.73e-3</td><td>4.2</td><td>2.71e-3</td><td>4.2</td><td>1.45e-2</td><td>2.0</td></tr><tr><td>40</td><td>1.23e-4</td><td>4.5</td><td>1.30e-4</td><td>4.4</td><td>1.30e-4</td><td>4.4</td><td>5.56e-3</td><td>1.4</td></tr><tr><td>80</td><td>8.47e-6</td><td><b>3.9</b></td><td>1.10e-5</td><td><b>3.6</b></td><td>1.26e-5</td><td><b>3.4</b></td><td></td><td></td></tr></tbody></table>

the SI-IMEX method behave similarly in this regime in terms of computational cost. These considerations are further confirmed by the experimental contraction rate (ECR), as defined in [67], i.e.

$$ECR(\tilde{k}) = \frac{\|p_{\tilde{k}+1} - p_{\tilde{k}}\|_{\infty}}{\|p_{\tilde{k}} - p_{\tilde{k}-1}\|_{\infty}}, \quad \tilde{k} > 0, \qquad (79)$$
<table><caption>Table 15: Convergence analysis for the traveling vortex test case using the SI-IMEX method with the IMEX-RK(6,6,4) scheme of type II (Table 23) and polynomial degree $r = 3$. Note that for $M = 10^{-1}$ the advective Courant number is $C_u \approx 0.05$, while for $M = 10^{-2}$ the advective Courant number is $C_u \approx 0.025$ to achieve a stable solution.</caption><thead><tr><th rowspan="2">N<sub>el</sub></th><th colspan="2">M = 10<sup>-1</sup></th><th colspan="2">M = 10<sup>-2</sup></th></tr><tr><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate δu</th><th>L<sup>2</sup> relative error δu</th><th>L<sup>2</sup> rate δu</th></tr></thead><tbody><tr><td>10</td><td>4.85e-2</td><td></td><td>4.88e-2</td><td></td></tr><tr><td>20</td><td>2.88e-3</td><td>4.1</td><td>2.74e-3</td><td>4.2</td></tr><tr><td>40</td><td>1.23e-4</td><td>4.5</td><td>1.30e-4</td><td>4.4</td></tr><tr><td>80</td><td>8.47e-6</td><td><b>3.9</b></td><td>1.10e-5</td><td><b>3.6</b></td></tr></tbody></table>

which is reported in Table 17. One can easily notice that the error is significantly reduced after three fixed point iterations as evident by the small value of the contraction rate. This results confirms that very few iterations of the fixed point loop are sufficient to obtain a satisfactory solution. We also monitor the experimental order of convergence of the fixed point method (EOC), defined as

$$EOC(\tilde{k}) = \frac{\log\left(\frac{\|p_{\tilde{k}+1} - p_{\tilde{k}}\|_{\infty}}{\|p_{\tilde{k}} - p_{\tilde{k}-1}\|_{\infty}}\right)}{\log\left(\frac{\|p_{\tilde{k}} - p_{\tilde{k}-1}\|_{\infty}}{\|p_{\tilde{k}-1} - p_{\tilde{k}-2}\|_{\infty}}\right)}, \quad \tilde{k} > 1 \qquad (80)$$

and we notice that, as expected, the fixed point method is of order 1 (Table 17). From now on, we set $\eta = 10^{-6}$ in the following test cases.

<table><caption>Table 16: Traveling vortex test case, investigation of the fixed point loop using the third order scheme of type II (Table 21) and polynomial degree $r = 2$. Here, $\eta$ denotes the tolerance for the stopping criterion of the fixed point loop (74), while $N$. iters represents the average number of iterations in the fixed point loop. The computational mesh of the reported results is composed by 80 × 80 = 6400 elements.</caption><thead><tr><th rowspan="2">η</th><th colspan="2">M = 10<sup>-1</sup></th><th colspan="2">M = 10<sup>-4</sup></th></tr><tr><th>L<sup>2</sup> relative error δu</th><th>N. iters</th><th>L<sup>2</sup> relative error δu</th><th>N. iters</th></tr></thead><tbody><tr><td>10<sup>-4</sup></td><td>3.35e-4</td><td>1</td><td>2.90e-4</td><td>1</td></tr><tr><td>10<sup>-6</sup></td><td>3.17e-4</td><td>3</td><td>2.90e-4</td><td>1</td></tr><tr><td>10<sup>-10</sup></td><td>3.17e-4</td><td>5</td><td>2.90e-4</td><td>1</td></tr></tbody></table>

## 4.3 Sod shock tube for the Peng-Robinson EOS

The main aim of the proposed time discretization methods is to use them for low Mach number flows so as to avoid the acoustic CFL restriction. In the case of high Mach number flows, acoustic waves are not negligible. Hence, one is interested in resolving both acoustic and material waves and explicit time discretization schemes are well suited to achieve this goal. However, we show that, when coupled with a monotonicity preserving spatial discretization, both IMEX and SI-IMEX methods can be effective also for high Mach number flows. We consider the Sod shock tube problem [83] for the Peng-Robinson EOS (9). The computational domain is
<table><caption>Table 17: Traveling vortex test case at $M = 10^{-1}$, investigation of the fixed point loop using the third order scheme of type II (Table 21) and polynomial degree $r = 2$. Here $k$ is the index of the fixed point iteration, ECR represents the experimental convergence rate (79), whereas EOC denotes the experimental order of convergence (80). The computational mesh of the reported results is composed by $80 \times 80 = 6400$ elements.</caption><thead><tr><th rowspan="2">k&#773;</th><th colspan="2">1st time step</th><th colspan="2">200th time step</th></tr><tr><th>||p&#773;<sub>k+1</sub> - p&#773;<sub>k</sub>||<sub>&infin;</sub> ECR</th><th>EOC</th><th>||p&#773;<sub>k+1</sub> - p&#773;<sub>k</sub>||<sub>&infin;</sub> ECR</th><th>EOC</th></tr></thead><tbody><tr><td>0</td><td>1.02e-5</td><td>-</td><td>5.84e-6</td><td>-</td></tr><tr><td>1</td><td>5.42e-6</td><td>0.53</td><td>5.34e-7</td><td>0.91</td></tr><tr><td>2</td><td>6.31e-8</td><td>0.012</td><td>4.25e-8</td><td>0.0080</td></tr><tr><td>3</td><td>1.10e-9</td><td>0.017</td><td>3.03e-10</td><td>0.0071</td></tr><tr><td>4</td><td>3.42e-11</td><td>0.031</td><td>4.68e-12</td><td>0.015</td></tr><tr><td>5</td><td>2.75e-12</td><td>0.08</td><td>1.42e-13</td><td>0.03</td></tr><tr><td>6</td><td>1.86e-13</td><td>0.07</td><td>7.55e-15</td><td>0.05</td></tr></tbody></table>

$\Omega = (-0.5, 0.5)$, whereas the final time is $T_f = 0.2$. The initial condition reads as follows:

$$ (\rho, u, p)(x, 0) = \begin{cases} (1, 0, 1) & \text{if } x < 0 \\ (0.125, 0, 0.1) & \text{if } x > 0. \end{cases} \qquad (81) $$

Dirichlet boundary conditions are imposed. Following [36], we take $c_v = 1$ and $R_g = 0.4$. Moreover, we take $b = 0.5$ and $a(T) = 0.5/\sqrt{T}$. We employ the IMEX-RK(2,2,2) scheme (Table 19) with polynomial degree $r = 0$ so as to avoid the oscillations that arise in the case of discontinuous solution when using high-order space discretization methods. The computational mesh is composed by $N_{el} = 500$ elements and the time step is $\Delta t \approx 1.33 \times 10^{-3}$, yielding a maximum acoustic Courant number $C \approx 1.23$ and $C_u \approx 0.47$. In particular, we point out that the acoustic Courant number is greater than one. A reference solution is computed using the optimal third order explicit strong stability-preserving scheme presented in [43] with $N_{el} = 8000$. An excellent agreement is established between the IMEX and the SI-IMEX method and a good agreement is established with the reference solution (Figure 7). Moreover, one can start appreciating the effectiveness of the linearization presented in Section 3.4. We will further discuss this point in Section 4.5.

Finally, we point out that the employed spatial discretization is not TVD for $r > 0$. Hence, spurious oscillations arise in the case of discontinuous solutions. A detailed discussion of possible approaches to overcome this issue is not in the scope of the present work. However, a number of approaches have been proposed in the literature to obtain essentially monotone schemes using high order DG methods, see, e.g., [37, 69].

## 4.4 Flow in an open tube

Next, we consider the test case III originally proposed in [58] for an ideal gas and also employed in [74]. The domain is $\Omega = (0, 10)$ with a time-dependent density and velocity prescribed at the left-end, while a time-dependent outflow pressure with a large amplitude variation is imposed at the right-end. More specifically, the initial conditions read as follows:

$$ (\rho, u, p)(x, 0) = (1, 1, 1), \qquad (82) $$
<figure><img src="image_8.png" /><figcaption>Figure 7: Sod shock tube for the Peng-Robinson EOS (9) with polynomial degree $r = 0$, a) density, b) velocity, c) pressure. The continuous black lines represent the reference solution computed using the third order optimal explicit strong stability-preserving scheme in [43], the continuous blue show the results obtained with the IMEX method using the IMEX-RK(2,2,2) scheme in Table 19, whereas the red dots report the numerical results obtained using the SI-IMEX method.</figcaption></figure>

while the boundary conditions are

$$ \rho(0,t) = 1 + \frac{3}{10} \sin(4t) \quad (83a) $$

$$ u(0,t) = 1 + \frac{1}{2} \sin(2t) \quad (83b) $$

$$ p(L,t) = 1 + \frac{1}{4} \sin(3t), \quad (83c) $$

with $L = 10$. The final time is $T_f = 7.47$, whereas the Mach number is set to $M = 10^{-4}$. Since (24) reduces to

$$ \nabla \cdot \bar{\mathbf{u}} = -\frac{1}{\gamma \bar{p}} \frac{d\bar{p}}{dt}, \quad (84) $$

the velocity field is not solenoidal as $M \to 0$, and, in one space dimension, it is a linear function of the space with a given time-dependent slope and boundary value at $x = 0$ [58, 74]. In particular, the leading order solution for the velocity reads as follows:

$$ \bar{u}(x,t) = u(0,t) - \frac{1}{\gamma p(L,t)} \frac{dp(L,t)}{dt} x = 1 + \frac{1}{2} \sin(2t) - \frac{3 \cos(3t)}{4\gamma \left(1 + \frac{1}{4} \sin(3t)\right)} x. \quad (85) $$
Moreover, the leading order relation (13) reduces to

$$ \frac{D \log \bar{\rho}}{Dt} = -\frac{\partial \bar{u}}{\partial x} = \frac{1}{\gamma p(L,t)} \frac{dp(L,t)}{dt} = \frac{3 \cos(3t)}{4\gamma \left(1 + \frac{1}{4} \sin(3t)\right)}, \quad (86) $$

where $\frac{D \log \bar{\rho}}{Dt} = \frac{\partial \log \bar{\rho}}{\partial t} + \bar{u} \frac{\partial \log \bar{\rho}}{\partial x}$.

First, we employ the second order IMER-RK(2,2,2) (Table 19) with polynomial degree $r = 1$. We consider a number of elements $N_{el} = 50$, whereas the time step is $\Delta t = 3.735 \times 10^{-3}$, leading to a maximum advective Courant number $C_u \approx 0.13$ and a maximum acoustic Courant number $C \approx 310$. The results at $t = \frac{T_f}{2}$ and at $t = T_f$ obtained with the IMEX method are those expected by the asymptotic analysis [58, 74] for both the density and velocity profiles (Figure 8), while the SI-IMEX method does not achieve a stable solution employing this time step. We also notice that no significant issue seems to arise in the low Mach regime.

The SI-IMEX method requires a smaller time step to achieve a stable solution. More specifically, we take $\Delta t = 1.8675 \times 10^{-3}$, which leads to a maximum acoustic Courant number $C \approx 155$ and a maximum advective Courant number $C_u \approx 0.065$. The use of the SI-IMEX method yields a computational time saving of around 37% at fixed time step. One can easily notice that a good agreement with the leading order solution is established for both the density and the velocity (Figure 8).

<figure><img src="image_9.png" /><figcaption>Figure 8: Open tube test case with the ideal gas law (4), comparison between the IMEX method and the SI-IMEX method using the second order scheme (Table 19) and polynomial degree $r = 1$ with $N_{el} = 50$. Top: results at $t = \frac{T_f}{2} = 3.735$. Bottom: results at $t = T_f = 7.47$. Left: density. Right: velocity. The continuous black lines show the leading order solution as $M \to 0$, the blue dots report the numerical results with the IMEX method and $\Delta t = 3.735 \times 10^{-3}$, whereas the red crosses represent the results with the SI-IMEX method and $\Delta t = 1.8675 \times 10^{-3}$.</figcaption></figure>

Next, we consider the third order scheme IMEX-RK(4,3,3) of type I (Table 20) and the third order scheme IMEX-RK(3,3,3) of type II (Table 21) using polynomial degree $r = 2$ with
$N_{el} = 50$. We take $\Delta t = 1.8675 \times 10^{-3}$, yielding a maximum advective Courant number $C_u \approx 0.14$ and a maximum acoustic Courant number $C \approx 292$. Similarly to the second order scheme, a good agreement with the leading order solution is established for both the time discretization schemes, in particular for the scheme of type I, while the SI-IMEX method does not achieve a stable solution. We take $\Delta t = 4.66875 \times 10^{-4}$, for which a stable solution for the scheme of type II is obtained. On the contrary, a stable solution is not achieved for the scheme of type I. This is likely related to the employed boundary conditions. The AP property and the AA property of the SI-IMEX method indeed has been proven considering periodic or no-flux boundary conditions [6, 17, 50], while we are considering time-dependent Dirichlet boundary conditions. Hence, the system is not autonomous and therefore it is outside of the theoretical framework of the SI-IMEX method depicted in Section 3.3. Since $c \neq \tilde{c}$ for schemes of type I, $U_E$ and $U_I$ are computed at different time instants, which seems to cause issues for time-dependent Dirichlet boundary conditions. An excellent agreement with the leading order solution is established for both the density and velocity field for $\Delta t = 4.66875 \times 10^{-4}$. Hence, for configurations which involve large variations of density and pressure, and time-dependent boundary conditions, the IMEX method is globally more robust and allows for sizeable larger time steps.

<figure><img src="image_10.png" /><figcaption>Figure 9: Open tube test case with the ideal gas law (4), comparison between the IMEX method and the SI-IMEX method using the third order schemes (Tables 20-21) in combination with polynomial degree $r = 2$ with $N_{el} = 50$. Top: results at $t = \frac{T_f}{2} = 3.735$. Bottom: results at $t = T_f = 7.47$. Left: density. Right: velocity. The continuous black lines show the leading order solution as $M \to 0$, the green squares report the numerical results with the IMEX method of type II and $\Delta t = 2.334375 \times 10^{-4}$, the orange circles show the numerical results with the IMEX method of type I and $\Delta t = 2.334375 \times 10^{-4}$, the red crosses represent the results with the SI-IMEX method of type II and $\Delta t = 2.334375 \times 10^{-4}$, the blue triangles show the results with the IMEX method of type II and $\Delta t = 1.245 \times 10^{-3}$, while the magenta diamonds report the results with the IMEX method of type I and $\Delta t = 1.245 \times 10^{-3}$.</figcaption></figure>
Next, we consider an extension of this test case using the stiffened gas equation of state (SG-EOS) (5). We take $\gamma = 4.4$, $\pi_{\infty} = 6.8 \times 10^2$, and $q_{\infty} = 0$ in (5). We employ the third-order scheme of type II (Table 21) in combination with $r = 2$. We take $\Delta t = 2.334375 \times 10^{-4}$, yielding a maximum acoustic Courant number $C \approx 1530$. An excellent agreement is established between the IMEX method and the SI-IMEX method (Figure 10). Moreover, one can easily notice that the leading order solution changes with the equation of state and modifying its parameters [74] (Figure 10). In particular, relation (24) for the SG-EOS reduces to

$$ \nabla \cdot \bar{\mathbf{u}} = - \frac{1}{\gamma (\bar{p} + \pi_{\infty})} \frac{d\bar{p}}{dt}. \qquad (87) $$

and therefore the leading order solution for the velocity reads as follows:

$$ \bar{u}(x, t) = 1 + \frac{1}{2} \sin(2t) - \frac{3 \cos(3t)}{4(\gamma + \pi_{\infty})(1 + \frac{1}{4} \sin(3t))} x. \qquad (88) $$

Hence, since from (83c) $\bar{p} \ge \frac{3}{4}$ and $|\frac{d\bar{p}}{dt}| \le \frac{3}{4}$, we obtain

$$ |\nabla \cdot \bar{\mathbf{u}}| \le \frac{1}{\gamma \left(\frac{3}{4} + \pi_{\infty}\right)} \frac{3}{4} \approx 2.5 \times 10^{-4}, \qquad (89) $$

meaning that the velocity field is almost constant.

<figure><img src="image_11.png" /><figcaption>Figure 10: Open tube test case with the SG-EOS (5), comparison between the IMEX method and the SI-IMEX method using the third order schemes (Table 21) and polynomial degree $r = 2$ with $N_{el} = 50$. Results at $t = T_f = 7.47$. Left: density. Right: velocity. The continuous black lines show the leading order solution as $M \to 0$, the continuous blue lines report the leading order solution as $M \to 0$ for the ideal gas with $\gamma = 1.4$, the green squares report the numerical results with the IMEX method, while the red crosses represent the results with the SI-IMEX method.</figcaption></figure>

## 4.5 Kelvin-Helmholtz instability at low Mach

In a final test, we consider the Kelvin-Helmholtz instability at low Mach studied, e.g., in [90], which we briefly recall here for the convenience of the reader. The computational domain is the square $\Omega = (-1, 1)^2$ endowed with periodic boundary conditions, while the final time is $T_f = 5$.
The initial conditions are

$$ \rho(\mathbf{x}, 0) = 1 - \frac{1}{4} \tanh \left[ 25 \left( |y| - \frac{1}{2} \right) \right] \qquad (90a) $$

$$ p(\mathbf{x}, 0) = \frac{10^4}{\gamma} \qquad (90b) $$

$$ u(\mathbf{x}, 0) = -\frac{1}{2} \tanh \left[ 25 \left( |y| - \frac{1}{2} \right) \right] \qquad (90c) $$

$$ v(\mathbf{x}, 0) = \frac{1}{100} \sin(2\pi x) \cos(2\pi y), \qquad (90d) $$

with $\gamma = 1.4$. Notice that, as discussed in [90], this configuration is in the incompressible regime, but the density is not constant. This shows that assuming a constant background density in (10), as done in some contributions, is not always a valid assumption, even in the incompressible regime. First, we employ the ideal gas law (4). We consider the IMEX-RK(2,2,2) scheme (Table 19) with polynomial degree $r = 1$ for the space discretization. The computational mesh is composed by $N_{el} = 120 \times 120 = 14400$ elements along each direction, whereas the time step is $\Delta t = 2 \times 10^{-3}$, yielding a maximum acoustic Courant number $C \approx 15$ and a maximum advective Courant number $C_u \approx 0.10$. The contours of density at $t = 2$ and at $t = 5$ are in good agreement with the reference results reported in [90] (Figure 11). Moreover, since we are analyzing a fluid mechanic instability, every small variation in the flow can lead to large variations [72, 70], and the excellent agreement obtained between the IMEX method and the SI-IMEX further confirms the correctness of our implementation and the properties of both methods.

Next, we employ the third order scheme of type II reported (Table 21) using polynomial degree $r = 2$ and $N_{el} = 80 \times 80 = 6400$ elements, so that the number of degrees of freedom does not change. Recall indeed that the total number of degrees of freedom per scalar variable is $N_{el} (r + 1)^2$. One can easily notice that the use of higher order methods is beneficial to resolve the fine details of the solution (Figure 12). An excellent agreement is once more established between the IMEX method and the SI-IMEX method (Figure 12).

Next, we consider an extension of this test case employing the general cubic EOS (8). Following [19], we set

$$ a(T) = \frac{1}{2\sqrt{T}} \qquad (91) $$

for the attraction term in (8). Moreover, we set $b = 5 \times 10^{-3}$. Finally, we assume that (9) is valid, with $c_v \approx 742.0$, and $R_g \approx 296.8$ in (6), which correspond to the thermodynamic properties of the nitrogen. An excellent agreement is obtained between the IMEX method and the SI-IMEX method, in spite of the linearization described in (54) (Figure 13). For the sake of completeness, we have also computed the solution obtained with the SI-IMEX method without linearizing the relation between internal energy and pressure, solving therefore a nonlinear equation for the pressure. No difference arises in the development of the Kelvin-Helmholtz instability (Figure 13) and a computational time saving of around 70% is obtained thanks to the linearization proposed in (54).

Finally, we consider two different configurations. First, we modify the initial condition, taking

$$ p(\mathbf{x}, 0) = \frac{10^5}{\gamma}, \qquad (92) $$

so that a more realistic maximum temperature $T \approx 320$ is obtained. In this case, we need to decrease the time step of the SI-IMEX method to $\Delta t = 5 \times 10^{-4}$ in order to achieve a stable solution. Next, we consider the proper expression of the attraction term $a(T)$ and of the
<figure><img src="image_12.png" /><figcaption>Figure 11: Kelvin-Helmholtz instability with the ideal gas law (4), results using the second order scheme (Table 19) and polynomial degree $r = 1$ with $N_{el} = 120$. Top: results at $t = 2$. Bottom: results at $t = T_f = 5$. Left: contour plots of the density field obtained with the IMEX method. Right: comparison between the IMEX method and the SI-IMEX method for the isoline equal to 1. The continuous black lines show the results with IMEX method, while the dashed red lines represent the results with the SI-IMEX.</figcaption></figure>

co-volume b for the Peng-Robinson EOS [39], [81, p. 263]:

$$ 
\left\{ 
\begin{aligned}
a(T) &= 0.45724 \frac{R_g^2 T_c^2}{p_c} \alpha(T)^2 \\
\alpha(T) &= 1 + \Gamma \left( 1 - \sqrt{\frac{T}{T_c}} \right) \\
\Gamma &= 0.37464 + 1.54226\omega - 0.26992\omega^2 \\
b &= 0.07780 \frac{R_g T_c}{p_c},
\end{aligned}
\right. 
\qquad (93) 
$$

where $T_c$ denotes the critical temperature, $p_c$ the critical pressure, and $\omega$ the acentric factor. For what concerns the nitrogen, we find $T_c = 126.19$ K, $p_c = 3.3978 \cdot 10^6$ Pa, and $\omega = 0.0372$ [64, 51]. Moreover, we consider the following relation for $e^\#(T)$ in (8) [28, 64]:

$$ e^\#(T) = \left[ A \frac{T}{1000} + \frac{1}{2} B \left( \frac{T}{1000} \right)^2 + \frac{1}{3} C \left( \frac{T}{1000} \right)^3 + \frac{1}{4} D \left( \frac{T}{1000} \right)^4 - E \frac{1000}{T} \right] \frac{10^6}{M_w} - R_g \omega $$
<figure><img src="image_13.png" /><figcaption>Figure 12: Kelvin-Helmholtz instability with the ideal gas law (4), results using the third order scheme of type II (Table 21) and polynomial degree $r = 2$ with $N_{el} = 80$ at $t = T_f = 5$. Left: contour plots of the density field obtained with the IMEX method. Right: comparison between the IMEX method, the SI-IMEX method, and the second order IMEX method for the isoline equal to 1. The continuous black line shows the results with third order IMEX method, the dashed red line represents the results with the third order SI-IMEX method, whereas the dashed-dotted blue line reports the results with the second order IMEX method.</figcaption></figure>

with $M_w = 28.0134 \text{ g mol}^{-1}$ and $A, B, C, D, E$ denoting suitable coefficients whose values are reported in Table 18. Notice that the polynomial expansion employed in [64] provides results expressed in kJ mol$^{-1}$. Hence, a proper conversion to obtain results in kJ kg$^{-1}$ has to be applied. The same consideration holds for the factor 1000, since the argument of the polynomial is expressed in thousandths of Kelvin. The IMEX method can achieve a stable solution for this configuration, whereas severe issues arise for the SI-IMEX method in the development of the instability.

<table><caption>Table 18: Values of the coefficients in (94)</caption><tr><td>A</td><td>28.98641</td></tr><tr><td>B</td><td>1.853978</td></tr><tr><td>C</td><td>-9.647459</td></tr><tr><td>D</td><td>16.63537</td></tr><tr><td>E</td><td>0.000117</td></tr></table>

# 5 Conclusions

Based on the experience of [17] and [74], we have performed a quantitative comparison between two different Implicit-Explicit Runge-Kutta (IMEX-RK) approaches for the Euler equations of gas dynamics. The two methods are particularly well suited for low Mach number flows, but keep their full accuracy for moderate values of the Mach number. The spatial discretization is based on the Discontinuous Galerkin (DG) method, which naturally allows for high-order accuracy, even though it is characterized by some limitations in the very low Mach limit on quadrilateral cells. The two schemes, namely the IMEX-DG method and the SI-IMEX-DG method, have been compared in a number of relevant benchmarks for ideal gases and on their
<figure><img src="image_14.png" /><figcaption>Figure 13: Kelvin-Helmholtz instability with the Peng-Robinson EOS (4), results using the third order scheme of type II (Table 21) and polynomial degree $r = 2$ with $N_{el} = 80$. Top-left: contour plots of the density field obtained with the IMEX method at $t = T_f = 5$. Top-right: comparison between the IMEX method and the SI-IMEX method for the isoline equal to 1. The continuous black lines show the results with IMEX method, while the dashed red lines represent the results with the SI-IMEX method. Bottom: comparison between the SI-IMEX method solving the nonlinear pressure equation and the SI-IMEX method using the linearization proposed in Section 3.4 for the isoline equal to 1. The continuous black lines show the results with SI-IMEX method, while the dashed red lines represent the results with the linearized SI-IMEX.</figcaption></figure>

non-trivial extension for non-ideal gases. The stiff dependence has been carefully analyzed in order to avoid the solution of a nonlinear pressure equation for a general class of equations of state (EOS).

First, we have assessed the convergence properties of the two methods. We have shown that they are asymptotic-preserving (AP) and asymptotically-accurate (AA) in the range of low Mach accuracy established by the spatial discretization. The SI-IMEX method provides a sizeable computational time saving ensuring the same level of accuracy, in particular for moderate values of the Mach number. Moreover, we have noticed an impact of the time discretization strategy and also of the specific time discretization scheme in the activation of spurious modes for low Mach numbers. More specifically, schemes of type I provide a superior stability for low values of the Mach number. A detailed analysis of the spatial discretization will be matter of future
work.

Next, we have considered non-trivial and less standard configurations. First, we have analyzed the case in which a time-dependent pressure is imposed at the boundary, for which the asymptotic limit does not coincide with the incompressible Euler equations. Notice that this configuration is outside of the theoretical framework of the SI-IMEX method and, indeed, only some IMEX-RK schemes of type II are suitable. The SI-IMEX method requires a significant smaller time step with respect to that needed by the IMEX method to achieve a stable solution.

Finally, we have focused on a Kelvin-Helmholtz instability at low Mach, that is in the incompressible regime, but for which the density is not constant. Here we have tested for non-ideal gases a novel linearization in the relation between internal energy and pressure, so as to avoid the solution of a nonlinear equation for the pressure. The proposed linearization is applicable to a general EOS and automatically recovers the linear systems obtained using the ideal gas law. No evident loss of accuracy occurs and a significant computational time saving is established.

In future work, as already mentioned, we aim to employ a spatial discretization based on simplices or Voronoi meshes that has been shown to be low Mach accurate for steady flows. Moreover, we aim to employ a spatial discretization based on compatible finite elements so as to improve the scaling properties with respect to the Mach number $M$. Finally, we aim to further analyze the stability properties of the two methods and to consider an extension of these approaches for the compressible Navier-Stokes equations and for two-phase flows.

## Acknowledgements

G. Orlando would like to acknowledge Vincent Perrier for useful discussions on related topics. G. Orlando, S. Boscarino, and G. Russo are part of the INdAM-GNCS National Research Group.

The simulations have been partly run at CINECA thanks to the computational resources made available through the ISCRA-C projects FEM-GPU - HP10CQYKJ1 and DGNWP - HP10C121HQ. We acknowledge the CINECA award, for the availability of high-performance computing resources and support.

G. Russo and S. Boscarino thank the Italian Ministry of Instruction, University and Research (MIUR) to support this research with funds coming from PRIN Project 2022, 2022KA3JBA, entitled "Advanced numerical methods for time dependent parametric partial differential equations and applications" and from PRIN 2022 PNRR "FIN4GEO: Forward and Inverse Numerical Modeling of hydrothermal systems in volcanic regions with application to geothermal energy exploitation.", (No. P2022BNB97). Both authors also have been supported for this work by the Spoke 1"FutureHPC&BigData" of the Italian Research Center on High-Performance Computing, Big Data and Quantum Computing (ICSC) funded by MUR Missione 4 Componente 2 Investimento 1.4: Potenziamento strutture di ricerca e creazione di "campioni nazionali di R&S (M4C2-19)" - Next Generation EU (NGEU).

## A Coefficients of employed IMEX-RK schemes

We report here for the convenience of the reader some information concerning the IMEX-RK schemes employed in the numerical simulations. We consider the second order IMEX-RK scheme proposed in [41] and also employed in [71, 74], whose coefficients are reported in the following Butcher tableaux
<table><caption>Table 19: *Butcher tableaux of the IMEX-RK(2,2,2) scheme in [41].* Left: explicit method. Right: implicit method.</caption><thead><tr><th>0</th><th>0</th><th>0</th><th>0</th><th>0</th><th>0</th><th>0</th><th>0</th></tr><tr><th>2 - &radic;2</th><th>2 - &radic;2</th><th>0</th><th>0</th><th>2 - &radic;2</th><th>1 - &radic;2/2</th><th>1 - &radic;2/2</th><th>0</th></tr><tr><th>1</th><th>1/2</th><th>1/2</th><th>0</th><th>1</th><th>&radic;2/4</th><th>&radic;2/4</th><th>1 - &radic;2/2</th></tr></thead><tbody><tr><td></td><td>&radic;2/4</td><td>&radic;2/4</td><td>1 - &radic;2/2</td><td></td><td>&radic;2/4</td><td>&radic;2/4</td><td>1 - &radic;2/2</td></tr></tbody></table>

For what concerns third order time discretization methods, we consider the following method of type I [13]:

<table><caption>Table 20: *Butcher tableaux of the IMEX-RK(4,3,3) scheme in [17].* Top: explicit method. Bottom: implicit method.</caption><tbody><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>0.435866521508</td><td>0.435866521508</td><td>0</td><td>0</td><td>0</td></tr><tr><td>0.717933260754</td><td>0.435866521508</td><td>0.282066739245</td><td>0</td><td>0</td></tr><tr><td>1</td><td>-0.733534082748750</td><td>2.150527381100</td><td>-0.416993298352</td><td>0</td></tr><tr><td></td><td>0</td><td>1.208496649176</td><td>-0.644363170684</td><td>0.435866521508</td></tr><tr><td>0.435866521508</td><td>0.435866521508</td><td>0</td><td>0</td><td>0</td></tr><tr><td>0.435866521508</td><td>0</td><td>0.435866521508</td><td>0</td><td>0</td></tr><tr><td>0.717933260754</td><td>0</td><td>0.282066739245</td><td>0.435866521508</td><td>0</td></tr><tr><td>1</td><td>0</td><td>1.208496649176</td><td>-0.644363170684</td><td>0.435866521508</td></tr><tr><td></td><td>0</td><td>1.208496649176</td><td>-0.644363170684</td><td>0.435866521508</td></tr></tbody></table>

and the following scheme of type II [55]:

<table><caption>Table 21: *Butcher tableaux of the IMEX-RK(3,3,3) scheme.* Top: explicit method. Bottom: implicit method.</caption><tbody><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>1767732205903</td><td>1767732205903</td><td>0</td><td>0</td><td>0</td></tr><tr><td>2027836641118</td><td>2027836641118</td><td>0</td><td>0</td><td>0</td></tr><tr><td>3/5</td><td>5535828885825</td><td>788022342437</td><td>0</td><td>0</td></tr><tr><td>1</td><td>10492691773637</td><td>10882634858940</td><td>0</td><td>0</td></tr><tr><td></td><td>6485989280629</td><td>-4246266847089</td><td>10755448449292</td><td>0</td></tr><tr><td></td><td>16251701735622</td><td>9704473918619</td><td>10357097424841</td><td>0</td></tr><tr><td></td><td>1471266399579</td><td>-4482444167858</td><td>11266239266428</td><td>1767732205903</td></tr><tr><td></td><td>7840856788654</td><td>7529755066697</td><td>11593286722821</td><td>4055673282236</td></tr><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>1767732205903</td><td>1767732205903</td><td>1767732205903</td><td>0</td><td>0</td></tr><tr><td>2027836641118</td><td>2027836641118</td><td>2027836641118</td><td>0</td><td>0</td></tr><tr><td>3/5</td><td>2746238789719</td><td>640167445237</td><td>1767732205903</td><td>0</td></tr><tr><td>1</td><td>10658868560708</td><td>6845629431997</td><td>4055673282236</td><td>0</td></tr><tr><td></td><td>1471266399579</td><td>-4482444167858</td><td>11266239266428</td><td>1767732205903</td></tr><tr><td></td><td>7840856788654</td><td>7529755066697</td><td>11593286722821</td><td>4055673282236</td></tr><tr><td></td><td>1471266399579</td><td>-4482444167858</td><td>11266239266428</td><td>1767732205903</td></tr><tr><td></td><td>7840856788654</td><td>7529755066697</td><td>11593286722821</td><td>4055673282236</td></tr></tbody></table>

Finally, we employ the fourth order time discretization method of type ARS proposed in [24]
<table><caption>Table 22: Butcher tableaux of the IMEX-RK(5,5,4) scheme. Top: explicit method. Bottom: implicit method.</caption><tbody><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>1/4</td><td>1/4</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>3/4</td><td>-1/4</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>11/20</td><td>-13/100</td><td>43/75</td><td>8/75</td><td>0</td><td>0</td><td>0</td></tr><tr><td>1/2</td><td>-6/85</td><td>42/85</td><td>179/1360</td><td>-15/272</td><td>0</td><td>0</td></tr><tr><td>1</td><td>0</td><td>79/24</td><td>-5/8</td><td>25/2</td><td>-85/6</td><td>0</td></tr><tr><td></td><td>0</td><td>25/24</td><td>-49/48</td><td>125/16</td><td>-85/12</td><td>1/4</td></tr><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>1/4</td><td>0</td><td>1/4</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>3/4</td><td>0</td><td>1/2</td><td>1/4</td><td>0</td><td>0</td><td>0</td></tr><tr><td>11/20</td><td>0</td><td>17/50</td><td>-1/25</td><td>1/4</td><td>0</td><td>0</td></tr><tr><td>1/2</td><td>0</td><td>371/1360</td><td>-137/2720</td><td>15/544</td><td>1/4</td><td>0</td></tr><tr><td>1</td><td>0</td><td>25/24</td><td>-49/48</td><td>125/16</td><td>-85/12</td><td>1/4</td></tr><tr><td></td><td>0</td><td>25/24</td><td>-49/48</td><td>125/16</td><td>-85/12</td><td>1/4</td></tr></tbody></table>

and the fourth order time discretization method of type II presented in [56]

<table><caption>Table 23: Butcher tableaux of the IMEX-RK(6, 6, 4) scheme. Top: explicit method. Bottom: implicit method.</caption><tbody><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>2γ</td><td>2γ</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>c3</td><td>ã31</td><td>ã32</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>c4</td><td>ã41</td><td>ã42</td><td>ã43</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>c5</td><td>ã51</td><td>ã52</td><td>ã53</td><td>ã54</td><td>0</td><td>0</td><td>0</td></tr><tr><td>c6</td><td>ã61</td><td>ã62</td><td>ã63</td><td>ã64</td><td>ã65</td><td>0</td><td>0</td></tr><tr><td>1</td><td>ã71</td><td>ã72</td><td>ã73</td><td>ã74</td><td>ã75</td><td>ã76</td><td>0</td></tr><tr><td></td><td>0</td><td>0</td><td>b3</td><td>b4</td><td>b5</td><td>b6</td><td>γ</td></tr><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>2γ</td><td>γ</td><td>γ</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>c3</td><td>a32</td><td>a32</td><td>γ</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>c4</td><td>a42</td><td>a42</td><td>a43</td><td>γ</td><td>0</td><td>0</td><td>0</td></tr><tr><td>c5</td><td>a52</td><td>a52</td><td>a53</td><td>a54</td><td>γ</td><td>0</td><td>0</td></tr><tr><td>c6</td><td>a62</td><td>a62</td><td>a63</td><td>a64</td><td>a65</td><td>γ</td><td>0</td></tr><tr><td>1</td><td>0</td><td>0</td><td>b3</td><td>b4</td><td>b5</td><td>b6</td><td>γ</td></tr><tr><td></td><td>0</td><td>0</td><td>b3</td><td>b4</td><td>b5</td><td>b6</td><td>γ</td></tr></tbody></table>
with

$$
\begin{align*}
\gamma &= 247/2000, & c_2 &= 2\gamma, \\
c_3 &= (2 + \sqrt{2})\gamma, & c_4 &= 67/200, \\
c_5 &= 3/40, & c_6 &= 7/10,
\end{align*}
$$

$$
\begin{align*}
b_3 &= 9164257142617/17756377923965, & b_4 &= -10812980402763/74029279521829, \\
b_5 &= 1335994250573/5691609445217, & b_6 &= 2273837961795/8368240463276, \\
b_7 &= 247/2000, \\
a_{32} &= 624185399699/4186980696204, & & \\
a_{42} &= 1258591069120/10082082980243, & a_{43} &= -322722984531/8455138723562, \\
a_{52} &= -436103496990/5971407786587, & a_{53} &= -2689175662187/11046760208243, \\
a_{54} &= 4431412449334/12995360898505, & & \\
a_{62} &= -2207373168298/14430576638973, & a_{63} &= 242511121179/3358618340039, \\
a_{64} &= 3145666661981/7780404714551, & a_{65} &= 5882073923981/14490790706663.
\end{align*}
$$

and

$$
\begin{align*}
\tilde{a}_{31} &= 247/4000, & \tilde{a}_{32} &= 2694949928731/7487940209513 \\
\tilde{a}_{41} &= 464650059369/8764239774964, & \tilde{a}_{42} &= 878889893998/2444806327765, \\
\tilde{a}_{43} &= -952945855348/12294611323341, & & \\
\tilde{a}_{51} &= 476636172619/8159180917465, & \tilde{a}_{52} &= -1271469283451/7793814740893, \\
\tilde{a}_{53} &= -859560642026/4356155882851, & \tilde{a}_{54} &= 1723805262919/4571918432560, \\
\tilde{a}_{61} &= 6338158500785/11769362343261, & \tilde{a}_{62} &= -4970555480458/10924838743837, \\
\tilde{a}_{63} &= 3326578051521/2647936831840, & \tilde{a}_{64} &= -880713585975/1841400956686, \\
\tilde{a}_{65} &= -1428733748635/8843423958496, & & \\
\tilde{a}_{71} &= 760814592956/3276306540349, & \tilde{a}_{72} &= 760814592956/3276306540349, \\
\tilde{a}_{73} &= -47223648122716/6934462133451, & \tilde{a}_{74} &= 71187472546993/9669769126921, \\
\tilde{a}_{75} &= -13330509492149/9695768672337, & \tilde{a}_{76} &= 11565764226357/8513123442827.
\end{align*}
$$

References

[1] E. Abbate, A. Iollo, and G. Puppo. “An asymptotic-preserving all-speed scheme for fluid dynamics and nonlinear elasticity”. *SIAM Journal on Scientific Computing* 41.5 (2019), A2850–A2879.

[2] P.C. Africa. “lifex: A flexible, high performance library for the numerical solution of complex finite element problems”. *SoftwareX* 20 (2022), p. 101252.

[3] D. Arndt et al. “The deal.II library, version 9.5”. *Journal of Numerical Mathematics* 31.3 (2023), pp. 231–246.

[4] D.N. Arnold. “An Interior Penalty Finite Element Method with Discontinuous Elements”. *SIAM Journal of Numerical Analysis* 19 (1982), pp. 742–760.

[5] D.N. Arnold, D. Boffi, and R.S. Falk. “Quadrilateral H(div) finite elements”. *SIAM Journal on Numerical Analysis* 42.6 (2005), pp. 2429–2451.

[6] K.R. Arun, A.J.D Gupta, and S. Samantaray. “Analysis of an asymptotic preserving low Mach number accurate IMEX-RK scheme for the wave equation system”. *Applied Mathematics and Computation* 411 (2021), p. 126469.

[7] K.R. Arun and S. Samantaray. “Asymptotic preserving low Mach number accurate IMEX finite volume schemes for the isentropic Euler equations”. *Journal of Scientific Computing* 82 (2020), pp. 1–32.
[8] U.M. Ascher, S.J. Ruuth, and R.J. Spiteri. "Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations". *Applied Numerical Mathematics* 25.2-3 (1997), pp. 151-167.

[9] W. Bangerth, R. Hartmann, and G. Kanschat. "deal.II: a general-purpose object-oriented finite element library". *ACM Transactions on Mathematical Software (TOMS)* 33 (2007), pp. 24-51.

[10] W. Barsukow. "Truly multi-dimensional all-speed schemes for the Euler equations on Cartesian grids". *Journal of Computational Physics* 435 (2021), p. 110216.

[11] F. Bassi and S. Rebay. "High-order accurate discontinuous finite element solution of the 2d Euler equations". *Journal of Computational Physics* 138 (1997), pp. 251-285.

[12] S. Boscarino. "Error analysis of IMEX Runge-Kutta methods derived from differential-algebraic systems". *SIAM Journal on Numerical Analysis* 45.4 (2007), pp. 1600-1621.

[13] S. Boscarino, F. Filbet, and G. Russo. "High order semi-implicit schemes for time dependent partial differential equations". *Journal of Scientific Computing* 68 (2016), pp. 975-1001.

[14] S. Boscarino, L. Pareschi, and G. Russo. *Implicit-explicit methods for evolutionary partial differential equations*. 1st ed. Society for Industrial and Applied Mathematics, 2024. DOI: 10.1137/1.9781611978209.

[15] S. Boscarino and G. Russo. "On a class of uniformly accurate IMEX Runge-Kutta schemes and applications to hyperbolic systems with relaxation". *SIAM Journal on Scientific Computing* 31.3 (2009), pp. 1926-1945.

[16] S. Boscarino and G. Russo. "Asymptotic preserving methods for quasilinear hyperbolic systems with stiff relaxation: a review". *SeMA Journal* 81.1 (2024), pp. 3-49.

[17] S. Boscarino et al. "High order semi-implicit WENO schemes for all-Mach full Euler system of gas dynamics". *SIAM Journal on Scientific Computing* 44 (2022), B368-B394.

[18] W. Boscheri, S. Busto, and M. Dumbser. "An all Mach number semi-implicit hybrid Finite Volume/Virtual Element method for compressible viscous flows on Voronoi meshes". *Computer Methods in Applied Mechanics and Engineering* 433 (2025), p. 117502.

[19] W. Boscheri and L. Pareschi. "High order pressure-based semi-implicit IMEX schemes for the 3D Navier-Stokes equations at all Mach numbers". *Journal of Computational Physics* 434 (2021), p. 110206.

[20] W. Boscheri et al. "A second order all Mach number IMEX finite volume solver for the three dimensional Euler equations". *Journal of Computational Physics* 415 (2020), p. 109486.

[21] W. Boscheri et al. "A structure-preserving staggered semi-implicit finite volume scheme for continuum mechanics". *Journal of Computational Physics* 424 (2021), p. 109866.
[22] L. Brugnano and V. Casulli. "Iterative solution of piecewise linear systems". *SIAM Journal on Scientific Computing* 30.1 (2008), pp. 463-472.

[23] J. Butcher. *Numerical Methods for Ordinary Differential Equations*. 2nd ed. Wiley, 2008.

[24] M.P. Calvo, J. De Frutos, and J. Novo. "Linearly implicit Runge-Kutta methods for advection-reaction-diffusion equations". *Applied Numerical Mathematics* 37.4 (2001), pp. 535-549.

[25] V. Casulli and D. Greenspan. "Pressure method for the numerical solution of transient, compressible fluid flows". *International Journal for Numerical Methods in Fluids* 4 (1984), pp. 1001-1012.

[26] V. Casulli and P. Zanolli. "A nested Newton-type algorithm for finite volume methods solving Richards' equation in mixed form". *SIAM Journal on Scientific Computing* 32.4 (2010), pp. 2255-2273.

[27] C. Chalons, M. Girardin, and S. Kokh. "Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms". *SIAM Journal on Scientific Computing* 35.6 (2013), A2874-A2902.

[28] M.W. Chase. "NIST-JANAF thermochemical tables for oxygen fluorides". *Journal of Physical and Chemical Reference Data* 25.2 (1996), pp. 551-603.

[29] A.J. Chorin. "Numerical solution of the Navier-Stokes equations". *Mathematics of Computation* 22 (1968), pp. 745-762.

[30] B. Cockburn, S.Y. Lin, and C.W. Shu. "TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. III. One-dimensional systems". *Journal of Computational Physics* 84 (1989), pp. 90-113.

[31] B. Cockburn and C.W. Shu. "TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. II. General framework". *Mathematics of Computation* 52 (1989), pp. 411-435.

[32] F. Cordier, P. Degond, and A. Kumbaro. "An asymptotic-preserving all-speed scheme for the Euler and Navier-Stokes equations". *Journal of Computational Physics* 231.17 (2012), pp. 5685-5704.

[33] R. Courant, K. Friedrichs, and H. Lewy. "Über die partiellen Differenzengleichungen der mathematischen Physik". *Mathematische annalen* 100.1 (1928), pp. 32-74.

[34] S. Dellacherie. "Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number". *Journal of Computational Physics* 229.4 (2010), pp. 978-1016.

[35] G. Dimarco, R. Loubère, and M.-H. Vignal. "Study of a new asymptotic preserving scheme for the Euler system in the low Mach number limit". *SIAM journal on Scientific Computing* 39.5 (2017), A2099-A2128.

[36] M. Dumbser and V. Casulli. "A conservative, weakly nonlinear semi-implicit finite volume scheme for the compressible Navier-Stokes equations with general equation of state". *Applied Mathematics and Computation* 272 (2016), pp. 479-497.
[37] M. Dumbser and R. Loubère. "A simple robust and accurate a posteriori subcell finite volume limiter for the discontinuous Galerkin method on unstructured meshes". *Journal of Computational Physics* 319 (2016), pp. 163-199.

[38] E. Feireisl et al. "On singular limits arising in the scale analysis of stratified fluid flows". *Mathematical Models and Methods in Applied Sciences* 26.03 (2016), pp. 419-443.

[39] M. Fernandez. "Propellant tank pressurization modeling for a hybrid rocket". MA thesis. Rochester Institute of Technology, 2009.

[40] F.X. Giraldo. *An Introduction to Element-Based Galerkin Methods on Tensor-Product Bases*. Springer Nature, 2020.

[41] F.X. Giraldo, J.F. Kelly, and E.M. Constantinescu. "Implicit-Explicit Formula-tions Of A Three-Dimensional Nonhydrostatic Unified Model Of The Atmosphere (NUMA)". *SIAM Journal of Scientific Computing* 35 (2013), pp. 1162-1194.

[42] E. Godlewski and P.-A. Raviart. *Numerical approximation of hyperbolic systems of conservation laws*. Vol. 118. Springer Science & Business Media, 2013.

[43] S. Gottlieb, C.-W. Shu, and E. Tadmor. "Strong stability-preserving high-order time discretization methods". *SIAM review* 43.1 (2001), pp. 89-112.

[44] W. Graebel. *Engineering fluid mechanics*. CRC Press, 2018.

[45] J.-L. Guermond et al. "On the implementation of a robust and efficient finite element-based parallel solver for the compressible Navier-Stokes equations". *Computer Methods in Applied Mechanics and Engineering* 389 (2022), p. 114250.

[46] H. Guillard and C. Viozat. "On the behaviour of upwind schemes in the low Mach number limit". *Computers & fluids* 28.1 (1999), pp. 63-86.

[47] J. Haack, S. Jin, and J.-G. Liu. "An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations". *Communications in Computational Physics* 12.4 (2012), pp. 955-980.

[48] E. Hairer, S. Norsett, and G. Wanner. *Solving Ordinary Differential Equations I: Nonstiff Problems*. Jan. 1993.

[49] A. Harten, P.D. Lax, and B. van Leer. "On upstream differencing and Godunov-type schemes for hyperbolic conservation laws". *SIAM review* 25.1 (1983), pp. 35-61.

[50] G. Huang, Y. Xing, and T. Xiong. "High order well-balanced asymptotic preserving finite difference WENO schemes for the shallow water equations in all Froude numbers". *Journal of Computational Physics* 463 (2022), p. 111255.

[51] R.T. Jacobsen, R.B. Stewart, and M. Jahangiri. "Thermodynamic properties of nitrogen from the freezing line to 2000 K at pressures to 1000 MPa". *Journal of Physical and Chemical Reference Data* 15.2 (1986), pp. 735-909.

[52] J. Jung and V. Perrier. "A curl preserving finite volume scheme by space velocity enrichment. Application to the low Mach number accuracy problem". *Journal of Computational Physics* 515 (2024), p. 113252.
[53] J. Jung and V. Perrier. "Behavior of the discontinuous Galerkin method for compressible flows at low Mach number on triangles and tetrahedrons". *SIAM Journal on Scientific Computing* 46.1 (2024), A452–A482.

[54] G.E. Karniadakis and S. Sherwin. *Spectral hp-Element Methods for Computational Fluid Dynamics*. Oxford University Press, 2005.

[55] C.A. Kennedy and M.H. Carpenter. "Additive Runge-Kutta schemes for convection-diffusion-reaction equations". *Applied Numerical Mathematics* 44 (2003), pp. 139–181.

[56] C.A. Kennedy and M.H. Carpenter. "Higher-order additive Runge-Kutta schemes for ordinary differential equations". *Applied numerical mathematics* 136 (2019), pp. 183–205.

[57] S. Klainerman and A. Majda. "Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids". *Communications on pure and applied Mathematics* 34.4 (1981), pp. 481–524.

[58] R. Klein. "Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics I: One-dimensional flow". *Journal of Computational Physics* 121.2 (1995), pp. 213–237.

[59] M. Kronbichler, T. Heister, and W. Bangerth. "High accuracy mantle convection simulation through modern numerical methods". *Geophysical Journal International* 191.1 (2012), pp. 12–29.

[60] M. Kronbichler and K. Kormann. "Fast matrix-free evaluation of discontinuous Galerkin finite element operators". *ACM Transactions on Mathematical Software (TOMS)* 45.3 (2019), pp. 1–40.

[61] V. Kučera et al. "Asymptotic properties of a class of linearly implicit schemes for weakly compressible Euler equations". *Numerische Mathematik* 150 (2022), pp. 79–103.

[62] O. Le Métayer and R. Saurel. "The Noble-Abel Stiffened-Gas equation of state". *Physics of Fluids* 28 (Apr. 2016), p. 046102.

[63] R.J. LeVeque. *Finite volume methods for hyperbolic problems*. Vol. 31. Cambridge university press, 2002.

[64] S.G. Lias et al. "NIST chemistry webbook standard reference database number 69". *National Institute of Standards Technology* (Jan. 2010).

[65] C.D. Munz. "On Godunov-type schemes for Lagrangian gas dynamics". *SIAM Journal on Numerical Analysis* 31.1 (1994), pp. 17–42.

[66] C.D. Munz et al. "The extension of incompressible flow solvers to the weakly compressible regime". *Computers and Fluids* 32 (2003), pp. 173–196.

[67] S. Noelle et al. "A weakly asymptotic preserving low Mach number scheme for the Euler equations of gas dynamics". *SIAM Journal on Scientific Computing* 36.6 (2014), B989–B1024.

[68] S. Ocher and F. Solomon. "Upwind difference schemes for hyperbolic conservation laws". *Math. Comput* 38 (1982), pp. 339–374.
[69] G. Orlando. "A filtering monotonization approach for DG discretizations of hyperbolic problems". *Computers & Mathematics with Applications* 129 (2023), pp. 113-125.

[70] G. Orlando. "An implicit DG solver for incompressible two-phase flows with an artificial compressibility formulation". *International Journal for Numerical Methods in Fluids* 96 (2024), pp. 1932-1959.

[71] G. Orlando, P.F. Barbante, and L. Bonaventura. "An efficient IMEX-DG solver for the compressible Navier-Stokes equations for non-ideal gases". *Journal of Computational Physics* 471 (2022), p. 111653.

[72] G. Orlando, T. Benacchio, and L. Bonaventura. "An IMEX-DG solver for atmospheric dynamics simulations with adaptive mesh refinement". *Journal of Computational and Applied Mathematics* 427 (2023), p. 115124.

[73] G. Orlando, T. Benacchio, and L. Bonaventura. "Impact of curved elements for flows over orography with a Discontinuous Galerkin scheme". *Journal of Computational Physics* 519 (2024), p. 113445.

[74] G. Orlando and L. Bonaventura. "Asymptotic-preserving IMEX schemes for the Euler equations of non-ideal gases". *Journal of Computational Physics* 529 (2025), p. 113889.

[75] L. Pareschi and G. Russo. "Implicit-explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation". *Journal of Scientific computing* 25 (2005), pp. 129-155.

[76] S. Patankar. *Numerical heat transfer and fluid flow*. CRC press, 2018.

[77] Ding-Yu Peng and Donald B Robinson. "A new two-constant equation of state". *Industrial & Engineering Chemistry Fundamentals* 15.1 (1976), pp. 59-64.

[78] S. Reddy et al. "Schur complement IMPLICIT-EXPLICIT formulations for discontinuous Galerkin non-hydrostatic atmospheric models". *Journal of Computational Physics* 491 (2023), p. 112361.

[79] F. Rieper and G. Bader. "The influence of cell geometry on the accuracy of upwind schemes in the low Mach number regime". *Journal of Computational Physics* 228.8 (2009), pp. 2918-2933.

[80] V. Rusanov. "The calculation of the interaction of non-stationary shock waves and obstacles". *USSR Computational Mathematics and Mathematical Physics* (1962), pp. 304-320.

[81] S.I. Sandler. *Chemical, biochemical, and engineering thermodynamics*. John Wiley & Sons, 2017.

[82] C.-W. Shu. "High order weighted essentially nonoscillatory schemes for convection dominated problems". *SIAM review* 51.1 (2009), pp. 82-126.

[83] G.A. Sod. "A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws". *Journal of computational physics* 27.1 (1978), pp. 1-31.
[84] M. Tavelli and M. Dumbser. "A pressure-based semi-implicit space-time discontinuous Galerkin method on staggered unstructured meshes for the solution of the compressible Navier-Stokes equations at all Mach numbers". *Journal of Computational Physics* 341 (2017), pp. 341-376.

[85] A. Thomann et al. "An all speed second order IMEX relaxation scheme for the Euler equations". *Communications in Computational Physics* 28 (2019), pp. 591-620.

[86] E.F. Toro. *Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction*. Springer Science & Business Media, 2009.

[87] J. Vidal. *Thermodynamics: Applications to chemical engineering and petroleum industry*. Editions Technip, 2001.

[88] G. Wanner and E. Hairer. *Solving ordinary differential equations II*. Vol. 375. Springer Berlin Heidelberg New York, 1996.

[89] G.B. Whitham. *Linear and nonlinear waves*. Vol. 42. John Wiley & Sons, 2011.

[90] E. Zampa and E. Dumbser. "An asymptotic-preserving and exactly mass-conservative semi-implicit scheme for weakly compressible flows based on compatible finite elements". *Journal of Computational Physics* 521 (2025), p. 113551.

[91] J. Zeifang et al. "A Novel Full-Euler Low Mach Number IMEX Splitting". *Communications in Computational Physics* 27 (2019), pp. 292-320.