Quasi all-speed schemes for magnetohydrodynamic flows in a wide range of Mach numbers

We present novel numerical schemes for ideal magnetohydrodynamic (MHD) simulations aimed at enhancing stability against numerical shock instability and improving the accuracy of low-speed flows in multidimensions. Stringent benchmark tests confirm that our scheme is more robust against numerical shock instability and is more accurate for low-speed, nearly incompressible flows than conventional shock-capturing schemes. Our scheme is a promising tool for tackling MHD systems, including both high and low Mach number flows.


Introduction
A magnetohydrodynamic (MHD) simulation is an indispensable tool for investigating macroscopic dynamics of laboratory, space, and astrophysical plasmas.In the field of compressible MHD simulations, a family of upwind, shock-capturing schemes have been developed based on the solution to the Riemann problem in one-dimensional hyperbolic conservation laws.These schemes enable us to effectively resolve shocks and discontinuities in supersonic flows.Numerous studies have been dedicated to the development of numerical schemes to ensure the robustness and accuracy of MHD simulations.Notably, the Roe [1,2] and the Harten-Lax-van Leer discontinuities (HLLD) [3] approximate Riemann solvers are extensively implemented in modern MHD simulation codes.
In practical MHD simulations, however, conventional shock-capturing schemes may encounter numerical difficulties.For high Mach number flows, these schemes tend to be vulnerable to numerical shock instability when a multidimensional shock closely aligns with the grid spacing [4].Conversely, in the case of low Mach number flows, these schemes lead to a degradation in solution accuracy as the Mach number decreases.In the field of computational aerodynamics, a family of the advection upstream splitting method (AUSM) [5,6] has been extensively employed due to its high flexibility, and has been extended to the "all-speed" regime to permit reliable simulations in a wide range of Mach numbers, including nearly incompressible flows and high Mach number shocks [7,8,9].
In this paper, we present novel numerical schemes designed for MHD simulations including both low and high Mach number flows.The first scheme is the multistate low-dissipation advection upstream splitting method (MLAU) [10], an MHD extension of the all-speed AUSM scheme.The second scheme is the low-dissipation HLLD scheme (LHLLD) [11], which incorporates techniques from the MLAU scheme, allowing users of the HLLD scheme to easily enjoy its capabilities.Both schemes provide accurate solutions for nearly incompressible flows that remain super-Alfvénic, a feature we refer to as quasi all-speed capability.

Numerical schemes
We consider one-dimensional ideal MHD equations written in the following conservative form: where U and F are the state vector of conservative variables and the corresponding flux vector, respectively, and ρ, u = (u, v, w), B = (B x , B y , B z ), and e are the mass density, velocity, magnetic field, and total energy density.The gas pressure P is determined from the equation of state for the ideal gas, where γ is the specific heat ratio.The solenoidal condition of the magnetic field gives B x = constant in one dimension.Equation ( 1) is discretized on computational cells into a finite volume form as follows: where F i±1/2 is the numerical flux at the interfaces of a cell [x i − ∆x/2, x i + ∆x/2].The quality of numerical solutions is significantly influenced by the assessment of the numerical flux.

Multistate low-dissipation advection upstream splitting method (MLAU)
Following the methodology of the AUSM-family schemes, we split the flux (equation ( 3)) into three parts, namely, advection, pressure, and magnetic tension parts, such that where, N = (0, 1, 0, 0, 0, 0, 0) T , and is the total pressure (excluding the contribution of B x ).Subsequently, the numerical flux at the interface is where the subscripts L and R refer to the left and right states at the interface.The mass flux ṁ and the pressure flux Pt are built based on the all-speed AUSM scheme [7], which utilizes a weighted average of the left and right states with respect the Mach number of the fast magnetosonic speed c f .The magnetic tension flux T is built to be consistent with the HLLD scheme.
Liou argued that the scheme becomes vulnerable to numerical shock instability in multidimensions when the mass flux ṁ includes the term proportional to the pressure difference (P i+1/2,R − P i+1/2,L ) and aligns with the shock surface; although, this term is effective for stabilizing one-dimensional shocks [12].To prevent the instability while maintaining the robustness of one-dimensional shocks, we introduce a shock-detecting factor θ and multiply it by the pressure difference term (see equation (20) in the next subsection), which quickly approaches zero when the shock surface is nearly parallel to the z − x or x − y plane (i.e., the shock normal direction is orthogonal to x).
The pressure flux Pt is expressed as a combination of the central average, pressure difference, and the velocity difference of the left and right states.While the pressure difference term scales with the fluid velocity and acts as upwinding, the velocity difference term scales with the fast magnetosonic speed and serves as diffusion in one-dimensional compressible flows.However, it causes excessive diffusion for low-speed flows in multidimensions.To correct the scaling of the pressure flux, we introduce a modified fast magnetosonic speed c u , where c 2 a = |B| 2 /ρ and c 2 ax = B 2 x /ρ are Alfvén speeds.We then multiples c u /c f by the velocity difference term.This correction appropriately scales down the pressure flux to |u| = c a , and as a result, it is expected to provide accurate solutions for super-Alfvénic flows, particularly in high-beta and low-speed plasma.For more comprehensive information about the scheme design, please refer to [10,11].

Low-dissipation HLLD approximate Riemann solver (LHLLD)
The novel techniques used in the MLAU scheme are implemented to the HLLD scheme to enhance stability against numerical shock instability and improve the accuracy of low-speed flows.The normal velocity S M and the total pressure P * t in the Riemann fan are modified to incorporate the shock detection in equation ( 15) and the pressure correction using equation ( 19), such that where The factor ϕ multiplied by the third term in equation ( 21) should satisfy the condition ∝ c u /c f (c u /c f ≪ 1) to improve the accuracy of low Mach number flows and equal 1 (c u /c f ≥ 1) to reduce to the original scheme at high Mach numbers.We use For more comprehensive information about the scheme design, please refer to [11].

Numerical experiments
We conduct stringent numerical experiments to assess the capability of the present schemes.Since the numerical results obtained using the MLAU and LHLLD schemes are nearly indistinguishable, we present the LHLLD scheme for clarity.Our numerical code incorporates the second-order MUSCL scheme with a minmod limiter for interpolating physical variables, the third-order strong stability preserving Runge-Kutta method for time integration, and the central upwind constrained transport method to preserve the solenoidal condition of the magnetic field [13].We employ a CFL number of 0.4.The first problem is the two-dimensional Kelvin-Helmholtz instability (KHI) in nearly incompressible flows.The initial condition has a velocity shear u = (V 0 /2) tanh((y − 10λ)/λ)e x , uniform density and pressure ρ = ρ 0 , P = P 0 , and a uniform magnetic field B = B 0 (cos θ, 0, sin θ), with ρ 0 = V 0 = B 0 = λ = 1.0 and θ = 71.565• .We use γ = 2.0.To initiate the instability, we impose a small perturbation into the y-component of the velocity along the shear layer with a wavelength of 20λ, corresponding to the fastest growing mode.The computational domain covers the range from 0 ≤ x < 20λ and 0 ≤ y < 20λ, and it is resolved by N × N cells, where N = 64, 128, 256, and 1024 for reference.The boundary condition is periodic and symmetric in the x and y directions.
Figure 1(a)-(c) shows the stream lines at t = 119 with P 0 = 500 (corresponding Mach number of 0.0158), obtained using the HLLD and LHLLD schemes at N = 256, and the HLLD scheme at N = 1024 as a reference.The in-plane magnetic field is amplified by stretching, and the tension force distorts the flow within the vortex in 5 < x < 15 and 7 < y < 13.It is evident that the LHLLD scheme outperforms the HLLD scheme in resolving the flow pattern, exhibiting good agreement with the reference solution that has 4 2 times higher grid resolution.Figure 1(d) shows the growth rate at N = 64, 128, 256 with P 0 = 500.The linear growth shown in a dashed line, obtained by solving the linearized MHD equations (equation (68) in [10]), is mostly independent of the pressure in nearly incompressible flows.The LHLLD scheme shown in red lines converges to a linear growth rate of γ = 0.051V 0 /λ for all runs, while the HLLD scheme shown in blue lines requires N = 256 to achieve convergence.Figure 1(e) shows the growth rate at N = 64 with P 0 = 500, 5000, 50000, and 500000.Thanks to the pressure correction (equation ( 21)), the solutions obtained using the LHLLD scheme (red lines) are mostly independent of the pressure, while the HLLD scheme (blue lines) gets worse as the pressure increases.Consequently, we can conclude that the LHLLD scheme possesses the capability to provide accurate solutions for nearly incompressible flows.
The second problem is the two-dimensional Richtmyer-Meshkov instability (RMI) occurring in hypersonic flows with a very weak magnetic field.The initial condition is as follows: (ρ, u, B, P ) = (1.0,0, −1.0, 0, 0.000034641, 0, 0, 0.00006) (y > 0) (3.9988, 0, −0.250075, 0, 0.000138523, 0, 0, 0.75) (y < 0) (23) which satisfies the Rankine-Hugoniot condition for perpendicular MHD shocks.The upstream Mach number is 100 and the plasma beta is β = 2P/|B| 2 = 10 5 .We use γ = 5/3.A corrugated contact discontinuity is imposed in the upstream region, y cd = Y 0 + ψ 0 cos(2πx/λ), where Y 0 = 1.0, ψ 0 = 0.1 is the corrugation amplitude, and λ = 1.0 is the wavelength.We shift a frame moving at v = −0.625,which corresponds to the interface velocity after the collision with the incident shock, ensuring that the structure of the RMI remains at approximately y = 0 throughout the simulation run.The computational domain covers the range from 0 ≤ x < λ and −40λ ≤ y < 40λ, and it is resolved by N × 80N cells, where N = 128.The boundary condition is periodic in the x direction and is fixed to the initial state in the y direction.
Figure 2 shows the density profile at t = 25 obtained using the Roe, HLLD, and LHLLD schemes.The surface at y = 4 corresponds to the transmitted shock, and the distinctive mushroom-shaped structure in 0 < y < 1 corresponds to the nonlinear evolution of the RMI.In panel (a), the solution obtained using the Roe scheme is severely distorted by numerical shock  instability.The instability generates multiple unphysical voids in the shocked region, disrupts the structure of the region of interest (0 < y < 1), and causes a loss of solution symmetry.In panel (b), the HLLD scheme significantly mitigates the catastrophic structures observed in the Roe scheme, but it still exhibits grid-scale oscillation around the shock front.In contrast, the LHLLD scheme in panel (c) succeeds in completely eliminating unphysical structures by utilizing the shock detection, resulting in a smooth solution.Furthermore, the solution obtained using the LHLLD scheme exhibits faster roll-up flows in the region 0 < y < 1 compared to the HLLD scheme, which is attributed to the pressure correction for subsonic flows in the Mach number range of 0.01 − 0.1.
The last problem is the three-dimensional MHD turbulence induced by the magnetorotational instability (MRI).Building upon our previous research [14], we solve isothermal MHD equations  including the Coriolis force within the framework of a local shearing box approximation [15].We apply the orbital advection method to update the velocity fluctuation u ′ = u − u 0 , where u 0 = −qΩxe y corresponds to the background shear flow (with q = 1.5 and Ω = 1 representing angular frequency) [16] .Since the velocity fluctuation typically remains subsonic, the solution may be affected by the pressure correction in the present scheme.The initial condition is uniform, (ρ, u ′ , B, P ) = (1, 0, B 0 e z , P 0 ), with an incompressible uniform random noise imposed to u ′ to initiate the instability.The ambient magnetic field B 0 is fixed to 0.025, and the gas pressure P 0 is set to 3.125, 312.5, 3125, yielding initial plasma beta values of 10 4 , 10 6 , 10 7 .The computational domain covers the range from −1 < x < 1, −2 < y < 2, and −0.5 < z < 0.5, and it is resolved by 64 × 256 × 32 cells.The wavelength of the fastest growing mode of the MRI, λ FGM = 2πB 0 / √ ρΩ = 0.16, is resolved by approximately five grid points.The boundary condition is periodic in the y and z directions, and shearing periodic in the x direction.Figure 3 shows the time profiles of the volume-averaged stress w = ρu ′ v ′ − B x B y /4π obtained using the HLLD and LHLLD schemes.Notably, these results exhibit discrepancies that become more pronounced as the gas pressure (plasma beta) increases.The time-and volume-averaged quantities over the time interval t = 25 − 50 are summarized in Table 1.For various pressure values, the stress values remain approximately constant around 0.065 with the LHLLD scheme, while they increase from 0.07 to 0.18 with the HLLD scheme.Considering that the flow remains subsonic within the Mach number range of 0.005 − 0.2, we can conclude that the pressureindependent solution provided by the LHLLD scheme is indicative of incompressible turbulence.This conclusion is supported by examining the time-and volume-averaged pressure perturbation ∆P = P − P 0 , presented in the fifth column of Table 1.The LHLLD scheme satisfies the asymptotic behavior in the incompressible limit ∆P/P 0 = O(M 2 ) = O(P −1 0 ).The discrepancy in the solutions between the HLLD and LHLLD schemes can be attributed to the dependence of the saturation level of the MRI-induced turbulence on the magnetic Prandtl number [17].Specifically, ideal MHD simulations of the MRI-induced turbulence are sensitive to the numerical magnetic Prandtl number, which is determined by the scheme design and is not necessarily specified [14].The numerical magnetic Prandtl number is defined as the ratio of the numerical viscosity ν n to resistivity η n , both of which arise from the numerical diffusion terms involved in the momentum and induction equations.For subsonic flows, the numerical viscosity is evaluated from the pressure in the Riemann fan (equation ( 21)), resulting in ν n,HLLD = c f ∆x and ν n,LHLLD = c u ∆x.As shown in Figure 4, the velocity profile reveals that small-scale structures become less pronounced in the HLLD scheme as the gas pressure increases, while the structure in the LHLLD scheme remains similar regardless of the pressure.Meanwhile, the numerical resistivity is identical in both schemes, and is evaluated from an explicit expression of the magnetic tension term, yielding η n,HLLD = η n,LHLLD = max (|u|, c ax ) ∆x [10].Consequently, it is anticipated that the numerical magnetic Prandtl number in the HLLD scheme increases with increasing pressure, while it remains approximately constant in the LHLLD scheme.
As a proxy for the numerical diffusion coefficients, we measure characteristic wavelengths defined as follows: where f represents the volume average of the variable f .The time-averaged characteristic wavelengths of u ′ and B y are presented in the sixth and seventh columns of Table 1.The wavelength of u ′ in the HLLD scheme tends to increase as the gas pressure increases, while other quantities remain approximately constant.This suggests that the numerical viscosity and the numerical magnetic Prandtl number in the HLLD scheme are likely to increase with increasing pressure, causing the increase in the saturation level [14]; in contrast, these values remain constant in the LHLLD scheme, resulting in the pressure-independent solution.

Conclusion
We presented novel numerical schemes designed for MHD flows in a wide range of Mach numbers.The schemes, the multistate low-dissipation advection upstream splitting method (MLAU) [10] and the low-dissipation HLLD approximate Riemann solver (LHLLD) [11], enhance stability against numerical shock instability and improve the accuracy of low-speed flows in multidimensions by incorporating shock detection and pressure correction.Stringent benchmark tests confirm the performance of the scheme, outperforming conventional shock-capturing schemes.The scheme provides accurate solutions for nearly incompressible flows that remain super-Alfvénic, by virtue of the pressure correction that scales down appropriately to |u| = c a .We refer to this capability as quasi all-speed for super-Alfvénic flows.Importantly, this feature ensures that the numerical magnetic Prandtl number remains approximately constant, a critical condition for the numerical study of MHD turbulent phenomena.With this novel ability, the scheme becomes a promising tool to tackle MHD systems, including both high and low Mach number flows.The source code written in C++ programming language is available on GitHub1 .

Figure 1 .
Figure 1.Numerical simulation of the Kelvin-Helmholtz instability: (a-c) Stream lines obtained using the HLLD scheme at N = 256, the LHLLD scheme at N = 256, and the HLLD scheme at N = 1024.(d-e) Time profiles of the fastest growing mode of the y-component of the velocity, demonstrating the dependence on grid resolution and the initial pressure.

Figure 2 .
Figure 2. Density profiles in the Richtmyer-Meshkov instability at t = 25, obtained using (a) the Roe scheme, (b) the HLLD scheme, and (c) the LHLLD scheme.The arrows represent the rotational velocity component.

7 Figure 3 .
Figure 3.Time profiles of the volume-averaged stress in the MRI-induced turbulence with plasma beta values of 10 4 (black lines), 10 6 (blue lines), and 10 7 (red lines), obtained using (a) the HLLD scheme and (b) the LHLLD scheme.