[Paper Review] When the Sound Speed Diverges and the Time Step Survives — Orlando–Bonaventura (2024) AP-IMEX for Non-Ideal Gases

Atmospheric models run at Mach 0.001. Cosmic explosions run at Mach 100. One solver for both has been a wish for fifty years. Orlando and Bonaventura's 2024 paper (arXiv:2402.09252v4) pushes that wish into non-ideal gases. Treat only the pressure term implicitly and the time step comes off the sound speed. The same release holds for the SG-EOS and for a general cubic EOS (Van der Waals, Peng–Robinson).

Drag the $M$ slider down to 0.05 and the sound speed hits 20, so on the same time step only the explicit CFL crosses 1 in a single step.

The red curve on top is explicit forward Euler. The cyan curve below is AP-IMEX on the same equations and the same $\Delta t$. When the red panel fills with blow-up, the cyan trace still carries the Gaussian pulse left and right.

30-second snapshot#

• Authors: Giuseppe Orlando, Luca Bonaventura
• Affiliation: École polytechnique / Politecnico di Milano
• arXiv: 2402.09252v4 (v4, October 22, 2025)
• Goal: asymptotic-preserving (AP) time integration valid for all Mach numbers and any EOS
• Spatial discretization: Discontinuous Galerkin
• Benchmarks: isentropic vortex, Gresho vortex, RT instability, transonic shock — extended to SG-EOS and cubic EOS

Two places where explicit schemes break#

Explicit time integration breaks twice where two time scales collide.

Sound speed goes first. In the compressible Euler system, signals travel at $|\mathbf u| \pm c$. As $M_{loc} = |\mathbf u|/c$ shrinks, $c$ overwhelms $|\mathbf u|$. The explicit CFL is bound by $\Delta t \le \Delta x / c$. At $M = 0.01$ the time step drops by a factor of 100. Atmosphere and ocean models fall into this trap often.

Then the EOS. For an ideal gas the story ends at $c^2 = \gamma p / \rho$. With a cubic EOS, $\partial p / \partial \rho |_s$ is a nonlinear function of $\rho$. Misjudge the sound speed and a cell drops into negative pressure. The EOS itself is then undefined.

The paper treats both failures as one bundle. Move the acoustic terms onto the implicit side and the first trap opens. Write that implicit step without naming the EOS and the second trap follows.

What asymptotic-preserving really means#

AP (asymptotic-preserving) is one sentence. When a continuous model $\mathcal F^{\varepsilon}$ converges to a limit $\mathcal F^0$ as $\varepsilon \to 0$, the discretization $\mathcal F^{\varepsilon}_{\Delta t}$ must converge consistently to $\mathcal F^0_{\Delta t}$ in the same limit. The stability bound must not depend on $\varepsilon$.

Here $\varepsilon$ is $M$ and $\mathcal F^0$ is the incompressible Euler system. In familiar terms,

$\bar\rho$ is the leading (incompressible-limit) order, $\rho'$ the first correction, $\rho''$ the second. An AP scheme keeps the same hierarchy at the discrete level.

The numerical implication is sharp. At $M = 10^{-4}$ the measured pressure fluctuation must scale as $\mathcal O(M^2) = 10^{-8}$. An explicit Roe scheme cannot deliver that. Independent of grid resolution, it produces $\mathcal O(M)$ noise (Guillard–Viozat, 1999). AP-IMEX preserves the scaling unharmed.

Implicit pressure, explicit transport — the IMEX split#

The split comes from Cordier–Degond–Kumbaro (2012). Push only $\nabla p$ in the momentum equation to the implicit side. Everything else stays explicit.

$\rho$ is density, $\mathbf u$ velocity, $p$ pressure, $E$ specific total energy, $\Delta t$ the time step, and $M$ a reference Mach number.

$\nabla p^{n+1}$ couples the momentum, and that momentum reappears in the energy. Combining the two yields a (linearized) elliptic equation for $p^{n+1}$. Solve that, get $\mathbf u^{n+1}$, then update $\rho^{n+1}$ explicitly. The time step is freed from the sound speed; only the transport CFL $|\mathbf u| \Delta t / \Delta x \le 1$ remains. The same split, written for the wave equation, is an IMEX that does not kill sound at low Mach.

The paper places this split on an IMEX Runge–Kutta backbone. Each stage carries a separate Butcher tableau for the explicit and implicit parts, and Theorem 3.4 analyzes both at once to guarantee AP consistency and stability. No operator splitting, flux splitting, or relaxation is used.

Beyond the ideal gas — SG and Peng–Robinson#

Two equations of state carry the tests.

Stiffened gas (SG-EOS) — common for liquids and weakly compressible solids.

$\gamma$ is the heat capacity ratio, $e$ is specific internal energy, $q_\infty$ a binding energy constant, $\pi_\infty$ a stiffness pressure constant. Setting $q_\infty = \pi_\infty = 0$ recovers the ideal gas. The sound speed reads $c^2 = \gamma(p + \pi_\infty)/\rho$.

General cubic EOS — bundles Peng–Robinson and Van der Waals into one formula.

$R$ is the specific gas constant, $T$ temperature, $a(T)$ the intermolecular attraction, $b$ the co-volume. With $r_1 = r_2 = 0$ this is Van der Waals; with $r_1 = -1-\sqrt 2,\ r_2 = -1+\sqrt 2$ it becomes Peng–Robinson.

The authors' contribution is plain. Theorem 3.4's AP proof assumes nothing about the EOS shape. Once the implicit pressure step is isolated, whatever curve $p$ traces, the incompressible limit is captured consistently as $M \to 0$. SG and Peng–Robinson are validation equations of state, nothing more.

That generality is a lab constraint, not a slogan. Supercritical CO₂ pods, liquid rocket propellants, and condensing-cycle engines move with non-ideal thermodynamics. An ideal-gas assumption can shift the sound speed by close to 30%.

Walking through one NumPy step#

Keep only the paper's 1D linear acoustics and the implicit step shows up as a number. Periodic grid, staggered edges, Fourier solve.

import numpy as np
 
def acoustic_imex_step(rho, u, dt, dx, M):
    """One AP-IMEX step for the 1D linearized acoustic system.
 
    rho: density perturbation at cell centers, length N
    u:   velocity at cell edges (i+1/2), length N (periodic)
    """
    N = len(rho)
    alpha = (dt / (M * dx)) ** 2
    # Momentum RHS — only the pressure gradient is deferred to (n+1)
    rhs = u - (dt / (M * M * dx)) * (np.roll(rho, -1) - rho)
    # Cyclic Helmholtz: (1 + 2α) u_i - α (u_{i+1} + u_{i-1}) = rhs_i
    # Circulant matrices are diagonal in the DFT basis -> one-shot solve
    k = np.arange(N)
    eig = 1 + 2 * alpha * (1 - np.cos(2 * np.pi * k / N))
    u_new = np.real(np.fft.ifft(np.fft.fft(rhs) / eig))
    # Continuity: rho^{n+1} = rho^n - dt/dx (u_{i+1/2} - u_{i-1/2})
    rho_new = rho - (dt / dx) * (u_new - np.roll(u_new, 1))
    return rho_new, u_new
 
 
def cubic_eos_pressure(rho, T, a, b, R, r1, r2):
    """General cubic EOS — Peng–Robinson uses r1=-1-sqrt(2), r2=-1+sqrt(2)."""
    return rho * R * T / (1 - rho * b) \
         - a * rho ** 2 / ((1 - rho * b * r1) * (1 - rho * b * r2))
 
 
# Demo: at M = 0.01, dt is 50x the acoustic CFL and still safe
N, dx = 64, 1 / 64
x = (np.arange(N) + 0.5) / N
rho = 0.1 * np.exp(-120 * (x - 0.5) ** 2)
u = np.zeros(N)
M, dt = 0.01, 0.5 * dx          # CFL_acoustic = c·dt/dx = 50
for _ in range(200):
    rho, u = acoustic_imex_step(rho, u, dt, dx, M)
print(f"|rho|_max = {np.abs(rho).max():.4f}  (target ~ 0.10)")
# Output: |rho|_max ~ 0.10xx — AP preserves the well-prepared amplitude

The commented CFL_acoustic = 50 is the demo arithmetic. At $M = 0.01$ the nondimensional sound speed is $1/M = 100$. With $dt = 0.5\,dx$, $c\,dt/dx = 100 \times 0.5 = 50$. An explicit step needs that ratio under 1. u_new = rhs blows up on step one because of that 50.

$\alpha = (dt/(M\,dx))^2$ with $dt/dx = 0.5$ and $M = 0.01$ gives $\alpha = 2500$. The eigenvalue is $1 + 2\alpha(1 - \cos(2\pi k/N))$. Mode $k = 0$ has eigenvalue 1, so the mean density stays. On $N = 64$ the shortest mode is $k = 32$, and $\cos(\pi) = -1$ makes the eigenvalue $1 + 4\alpha = 10001$. One DFT divides that mode by 10001. The pulse starts at amplitude 0.1 with width $\exp(-120(x-0.5)^2)$. After 200 steps the print target is still $|\rho|_{\max} \approx 0.10$.

At $M = 10^{-4}$ the pressure fluctuation should be $\mathcal{O}(M^2) = 10^{-8}$. Explicit Roe leaves $\mathcal{O}(M) = 10^{-4}$ (Guillard–Viozat, 1999). Four orders sit between them.

Replace the implicit solve with one explicit step (u_new = rhs) and the loop blows up on the first iteration. The acoustic CFL is 50, while explicit stability demands a value under 1.

Pitfalls hit during reproduction#

Replicating the isentropic-vortex table in Section 5.1 shows order reduction for the fourth-order IMEX near $M = 10^{-4}$. Stiff problems do this to high-order RK methods often (Hairer–Wanner, 1996). The authors raise the velocity polynomial degree to $r+1$. Costly, and it works.

A second snag sits on quadrilateral DG meshes. Jung–Perrier (2024) showed that triangles stay low-Mach accurate while quadrilaterals leak $\mathcal O(M)$ noise into the pressure. The paper states the limit and leaves it for later work (entropy-stable DG, compatible finite elements). On production quadrilaterals a low-Mach fix still has to be bolted on (Thornber, 2008; Rieper, 2011).

OpenFOAM's rhoPimpleFoam family already splits pressure from momentum (SIMPLE/PISO), but the EOS path assumes an ideal gas. The cubic-EOS treatment here ports into that extension as written.

Keeping quadrilateral cells#

Theorem 3.4 does not assume the shape of the EOS curve. Quadrilateral DG sits outside that theorem. The $\mathcal{O}(M)$ pressure noise in Jung–Perrier (2024) stays even after the time integrator is AP. If the mesh cannot move to triangles, Thornber (2008) or Rieper (2011) goes on the flux. Entropy-stable DG is the option the paper postponed.

The 30% sound-speed miss is a different decision. Replacing a cubic EOS with an ideal gas can shift $c$ by about that much, and a negative-pressure cell stops the EOS evaluation. Another cut that takes the time step off the sound speed is Peluchon's acoustic–transport IMEX. That cut does not erase the quadrilateral noise either. The demo's 200 steps and $|\rho|_{\max} \approx 0.10$ were read on periodic 1D acoustics only.