© G. Leidi et al. 2022

1 Introduction

The interplay between turbulent convection and shear is fundamental in understanding the role of small- and large-scale dynamo mechanisms in the generation of strong magnetic fields in stellar interiors. These processes can only be modeled self-consistently with multidimensional magnetohydrodynamic (MHD) simulations (Brun et al. 2004; Browning et al. 2006; Browning 2008; Brown et al. 2010; Ghizaru et al. 2010; Käpylä et al. 2012; Masada et al. 2013; Karak et al. 2015; Hotta et al. 2015; Yadav et al. 2016; Augustson et al. 2016; Brun & Browning 2017; Rempel 2018; Käpylä 2021). Nowadays, many codes used for astrophysical MHD rely on finite-volume discretization and Godunov-like methods to retain the conservative property of the MHD equations. This method is particularly suited for simulating flows in the transonic and supersonic regimes, which characterize many astrophysical systems. However, stars are objects in nearly magnetohydrostatic equilibrium (MHSE), and the flows arising from such stratifications have a very low sonic Mach number, typically ℳ_son = ∣V∣/a ≲ 10⁻² (Kupka & Muthsam 2017), where V is the flow speed and a is the adiabatic sound speed. It is well known that conventional finite-volume schemes are not designed to work in such regimes (Viallet et al. 2011; Miczek et al. 2015; Dumbser et al. 2019; Minoshima et al. 2020). First, the approximate Riemann solvers used in many astrophysical MHD codes, such as the Harten–Lax–van Leer (HLLE; Einfeldt et al. 1991), Roe (Cargo & Gallice 1997), and Harten–Lax–van Leer discontinuities (HLLD; Miyoshi & Kusano 2005) solvers, show excessive numerical dissipation when the typical Mach number of the flow is below 10⁻². Second, explicit time-steppers have to satisfy the Courant–Friedrichs–Lewy (CFL) stability criterion (Courant et al. 1928) so that the propagation of fast magnetosonic waves is resolved in time. This poses a severe limitation when simulating low-Mach-number flows. In this regime, the fast magnetosonic waves become parasitic since they transport very little energy and drastically reduce the time step. This makes convectional schemes exceedingly expensive for simulating the evolution of fluid motions and Alfvén waves, which are orders of magnitude slower than fast magnetosonic waves in deep layers of stars (Brun et al. 2005; Browning 2008; Käpylä 2011; Augustson et al. 2016). Lastly, standard Godunov-type schemes, by construction, cannot preserve stratifications in MHSE. This happens because hyperbolic fluxes and gravitational source terms are separately discretized and do not balance exactly in hydrostatic setups, which inevitably leads to the generation of spurious flows even in pure hydrodynamic simulations (Edelmann et al. 2021). This problem becomes even more critical in steep stratifications, where the spatial reconstruction from cell centers results in large jumps in the pressure at the cell interfaces, considerably accelerating the fluid along the gravity vector. Such numerical artifacts can dominate over the physical convective motions, leading to unreliable results.

Difficulties in modeling low-Mach-number MHD flows in stellar interiors are usually overcome by using alternative approaches based on a different formulation of the physical problem. One of them consists in artificially boosting the energy flux (or energy generation) to drive faster convective motions, such that the typical Mach number of the resulting flows falls above the low-Mach regime (ℳ_son ≳ 10⁻²), where explicit time-steppers can be used efficiently (Käpylä et al. 2011, 2012, 2013; Viviani et al. 2019; Käpylä 2021). However, this approach significantly enhances the relative fluctuations of thermodynamic quantities (Warnecke et al. 2016; Käpylä et al. 2020), alters the mixing at convective boundaries (Hotta 2017; Käpylä 2019), and modifies the spectrum of internal waves in radiative regions of stars (Rogers et al. 2013; Edelmann et al. 2019; Horst et al. 2020; Higl et al. 2021). Another approach consists in solving the set of MHD equations using the anelastic approximation (Glatzmaier 1984, 1985; Brun et al. 2004; Jones et al. 2009; Gastine & Wicht 2012; Smolarkiewicz & Charbonneau 2013; Featherstone & Hindman 2016), which filters out the fast magnetosonic waves, alleviating the overly strict constraint on the time step. However, such an approximation cannot model the excitation of the compressible pressure modes. Another way to overcome the constraint on the time step is to drastically reduce the speed of the fast magnetosonic waves (Rempel 2005; Hotta et al. 2015); again, at the cost of modifying the original set of MHD equations.

A numerical scheme that is capable of efficiently solving the fully compressible MHD equations at low sonic Mach numbers in strongly stratified setups is still missing. In this work we present a new method that aims to fill this gap. This can only be accomplished (i) by reducing the numerical dissipation, (ii) by overcoming the strict CFL condition, and (iii) by preserving the background stratification in MHSE over long timescales. For aspect (i), we use a low-Mach version of the five-wave HLLD solver (Minoshima & Miyoshi 2021), whose numerical dissipation is independent of the sonic Mach number of the modeled flow in subsonic regimes.

In order to deal with aspect (ii), time-implicit discretization techniques can be used. Most of the fully implicit (e.g., Aydemir & Barnes 1985; Charlton et al. 1990; Chacón 2008; Lütjens & Luciani 2010) and semi-implicit (e.g., Harned & Kerner 1985; Schnack et al. 1987; Lerbinger & Luciani 1991; Glasser et al. 1999; Jardin 2012; Fambri 2021) MHD schemes presented in the literature are designed to simulate magnetically confined plasmas and low-β environments, where β is defined as the ratio of the gas pressure to the magnetic pressure: β = p/p_B. In such strongly magnetized plasmas, both fast magnetosonic and Alfvén waves put a strong limit on the time step. On the contrary, plasmas in stellar interiors are characterized by high β values (Mestel 1999)¹, and the only source of stiffness is the generation of fast magnetosonic waves. Thus, Alfvén waves do not need to be treated implicitly, which greatly simplifies the numerical problem. Recently, Dumbser et al. (2019) developed a semi-implicit conservative method that treats only the fast magnetosonic waves implicitly; however, that scheme cannot be implemented easily within the framework of our hydrodynamic code. In this work, we construct an alternative time-marching scheme suitable for modeling high-β plasmas at low Mach numbers, based on the approach described by Fuchs et al. (2009), in which the induction equation is solved in a separate step and coupled to the rest of the system through Strang splitting (Strang 1968). As the high speed of the fast magnetosonic waves is mostly determined by the pressure flux in the momentum equation, we solve the subset containing the continuity, momentum, and energy equations implicitly, whereas the induction equation is integrated using an explicit time-stepper. For stability, the time step is now limited by the fastest fluid and Alfvén speeds on the grid; it is approximately 1/ℳ_son longer than that allowed by the CFL condition, which leads to a considerable speed-up when the Mach number of the flow is low. Since the update on the induction equation is performed in a separate step, the flux-Jacobian in the time-implicit part of the algorithm does not need to be evaluated with respect to the magnetic field components. This allows for more flexibility when choosing the method that evolves the magnetic field. In particular, we use a staggered formulation of constrained transport (CT-contact; Gardiner & Stone 2005) to keep ∇ · B = 0 to machine precision, at least for a specific discretization of the divergence of the magnetic field.

Finally, aspect (iii) is solved by using the deviation well-balancing method (Berberich et al. 2021; Edelmann et al. 2021), which allows the a priori known background stratification in MHSE to be preserved, dramatically reducing the magnitude of numerical errors and the strength of spurious flows². Recently, Canivete Cuissa & Teyssier (2022) performed fully compressible simulations of stellar magneto-convection at ℳ_son ~ 10⁻³ using a well-balancing technique similar to the deviation method. However, their scheme relied on explicit time-steppers, and it was used to simulate only 2.5 convective turnovers. Moreover, they did not cure the excessive dissipation of the HLLD solver at low Mach numbers.

These methods have been implemented in the SEVEN-LEAGUE HYDRO (SLH) code, which has already been used in the past to simulate convective boundary mixing, shear instabilities, and wave excitation in stellar interiors, even in regimes of low Mach numbers (Miczek 2013; Edelmann 2014; Miczek et al. 2015; Edelmann & Röpke 2016; Edelmann et al. 2017; Horst et al. 2020, 2021; Andrassy et al. 2022). We stress that the current MHD implementation in SLH is not suitable for modeling low-β plasmas. For simulating such regimes, a different method should be used instead, which is beyond the scope of this work.

In Sect. 2, we summarize the main properties of the fully compressible MHD equations with gravity. In Sects. 3 and 4, we provide details on the numerical algorithms implemented in SLH. In Sect. 5, several numerical experiments are run with the new MHD scheme in order to check its accuracy and efficiency in simulating flows at low Mach numbers, even in the presence of a steep stratification. Finally, in Sect. 6 we draw conclusions and summarize the fundamental aspects of the proposed algorithm.

2 Equations of compressible ideal MHD with gravity

The MHD scheme implemented in SLH is designed to solve the set of compressible ideal MHD equations with a (time-independent) gravitational source term³: $$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho V\right)=0,$$(1) $$\frac{\partial \left(\rho V\right)}{\partial t}+\nabla \cdot \left[\rho V\otimes V+\left(p+p_{\text{B}}\right)I-B\otimes B\right]=\rho g,$$(2) $$\frac{\partial \left(\rho E_{\varphi }\right)}{\partial t}+\nabla \cdot \left[\left(\rho E_{\varphi }+p+p_{\text{B}}\right)V-B\left(B\cdot V\right)\right]=0,$$(3) $$\frac{\partial B}{\partial t}+\nabla \cdot \left(V\otimes B-B\otimes V\right)=0,$$(4)

where ρ denotes the density, V = (V_x, V_y, V_z) the velocity field, B = (B_x, B_y, B_z) the magnetic field⁴, I the unit tensor, g = (g_x, g_y, g_z) the gravitational acceleration, p the gas pressure and p_B = ∣B∣²/2 the magnetic pressure. The total energy density ρE_ϕ is defined as $$\rho E_{\varphi }=\rho e_{\mathrm{int}}+\frac{1}{2}\rho {\left|V\right|}^2+\frac{1}{2}{\left|B\right|}^2+\rho e_{\varphi },$$(5)

where e_int and e_ϕ are the specific internal and gravitational energies⁵.

The system is closed by an equation of state (EoS), which provides the numerical value of the gas pressure. Several different definitions for the EoS can be used in SLH, including a simple ideal gas law, radiation pressure, and a tabulated EoS (Helmholtz EoS; Timmes & Swesty 2000) that allows the effects of electron degeneracy and Coulomb corrections to be included, which are often needed to properly describe the thermodynamic conditions found in stellar material.

2.1 Eigenstructure of the MHD system and definition of Mach numbers

The homogeneous MHD system (left-hand side of Eqs. (1)–(4)) reduced to one spatial dimension has seven nonzero eigenvalues⁶, $${\lambda }_{1,7}=V_x\mp c_{\text{f},x},{\lambda }_{2,6}=V_x\mp c_{\text{a},x},{\lambda }_{3,5}=V_x\mp c_{\text{s},x},{\lambda }_4=V_x,$$(6)

associated with different modes of propagation: left/right fast magnetosonic waves, left/right Alfvén waves, left/right slow magnetosonic waves and one entropy wave. c_s,x, c_a,x, and c_f,x are the slow magnetosonic, Alfvén, and fast magnetosonic speeds, $$c_{\text{f},x}={\left[\frac{1}{2}\left(a^2+\frac{{\left|B\right|}^2}{\rho }\pm \sqrt{{\left(a^2\frac{{\left|B\right|}^2}{\rho }\right)}^2-4a^2{c}_{\text{a},x}^2}\right)\right]}^{\frac{1}{2}},$$(7) $$c_{\text{a},x}=\left|B_x\right|/\sqrt{\rho },$$(8)

with the adiabatic sound speed a defined as $$a={\left(\frac{\partial p}{\partial \rho }+\frac{p}{{\rho }^2}\frac{\partial p}{\partial e_{\mathrm{int}}}\right)}^{\frac{1}{2}}.$$(9)

As illustrated in Fig. 1, the MHD wave pattern has a fixed ordering: $${\lambda }_1<{\lambda }_2<{\lambda }_3<{\lambda }_4<{\lambda }_5<{\lambda }_6<{\lambda }_7.$$(10)

Alfvén waves correspond to incompressible modes of propagation, as they only carry perturbations in the velocity and magnetic field components orthogonal to the wave vector. Effects of compressibility are due to the propagation of slow and fast magnetosonic waves, while the entropy wave is simple advection of fluid.

In contrast to pure hydrodynamics, the more complex structure of the MHD waves allows several Mach numbers to be defined, depending on the considered reference velocity. In addition to ℳ_son (see Sect. 1), the Alfvén Mach number is defined as $${ℳ}_{\text{Alf}}=\frac{\left|V\right|}{\sqrt{c_{\text{a},x}^2+c_{\text{a},y}^2+c_{\text{a},z}^2}},$$(11)

while the directional slow and fast magnetosonic Mach numbers are given by $${ℳ}_{\text{slow,fast,}x}=\frac{\left|V_x\right|}{c_{\text{s,f},x}}.$$(12)

Fig. 1 Wave structure of the MHD system.

2.2 Low-Mach limit of the MHD system

Magnetic fields amplified by dynamo mechanisms in deep convective layers of stars are likely to approach equipartition with respect to the kinetic energy content of the flow (Brandenburg & Subramanian 2005; Featherstone et al. 2009; Augustson et al. 2016; Hotta 2017; Käpylä 2019). To model such processes, the MHD system in Eqs. (1)–(4) must be solved in regimes of Mach numbers ℳ_fast,x ≲ ℳ_son ≤ 10⁻² and ℳ_Alf ~ 1. To infer the structure of the solution under such conditions, it is useful to consider the nondimensional form of the fully compressible MHD equations⁷, $$\frac{\partial \widehat{\rho }}{\partial \widehat{t}}+\widehat{\nabla }\cdot \left(\widehat{\rho }\widehat{V}\right)=0,$$(13) $$\frac{\partial \left(\widehat{\rho }\widehat{V}\right)}{\partial \widehat{t}}+\widehat{\nabla }\cdot [\widehat{\rho }\widehat{V}\otimes \widehat{V}+\left(\frac{\widehat{p}}{{\widehat{ℳ}}_{\text{son}}^2}+\frac{{\widehat{p}}_B}{{\widehat{ℳ}}_{\text{Alf}}^2}\right)I-\frac{\widehat{B}\otimes \widehat{B}}{{\widehat{ℳ}}_{\text{Alf}}^2}]=0,$$(14) $$\frac{\partial \left(\widehat{\rho }\widehat{E}\right)}{\partial \widehat{t}}+\widehat{\nabla }\cdot \left[\left(\widehat{\rho }\widehat{E}+\widehat{p}+{\widehat{p}}_B\frac{{\widehat{ℳ}}_{\text{son}}^2}{{\widehat{ℳ}}_{\text{Alf}}^2}\right)\widehat{V}-\widehat{B}\left(\widehat{B}\cdot \widehat{V}\right)\frac{{\widehat{ℳ}}_{\text{son}}^2}{{\widehat{ℳ}}_{\text{Alf}}^2}\right]=0,$$(15) $$\frac{\partial \widehat{B}}{\partial \widehat{t}}+\widehat{\nabla }\cdot \left(\widehat{V}\otimes \widehat{B}-\widehat{B}\otimes \widehat{V}\right)=0.$$(16)

Here, the different variables have been rescaled by some reference quantity representative of the physical system of interest: $x=\widehat{x}{x}_{\text{r}}$, $\rho =\widehat{\rho }{\rho }_{\text{r}}$, $V=\widehat{V}{V}_{\text{r}}$, $E=\widehat{E}{a}_{\text{r}}^2$, $p=\widehat{p}{\rho }_{\text{r}}{a}_{\text{r}}^2$, $B=\widehat{B}{B}_{\text{r}}$. ${\widehat{ℳ}}_{\text{son}}=\left|V_{\text{r}}\right|/a_{\text{r}}$ and ${\widehat{ℳ}}_{\text{Alf}}=\left|V_{\text{r}}\right|/\left(\left|B_{\text{r}}\right|/\sqrt{{\rho }_{\text{r}}}\right)$ are the characteristic sonic and Alfvén Mach numbers of the flow.

In the limit of ${\widehat{ℳ}}_{\text{son}}\to 0$, Eqs. (13)–(15) approach the incompressible regime (see Matthaeus & Brown 1988), in which the gas pressure is constant in space except for fluctuations ∝ ${\widehat{ℳ}}_{\text{son}}^2$. As this solution does not allow for compressible modes of propagation, only Alfvén and entropy waves can transport fluctuations in the state variables across the physical domain. In the regime we are interested in $\left({\widehat{ℳ}}_{\text{Alf}}~1\right)$, these waves travel at similar speeds. However, the incompressible limit is not the only solution to the compressible MHD equations at ${\widehat{ℳ}}_{\text{son}}\ll 1$. In fact, both slow and fast magnetosonic waves can be propagated with arbitrary small velocity fluctuations. Fast magnetosonic waves in particular travel at much higher speed than Alfvén waves and fluid motions in the low-Mach limit. If the plasma-β is high, the large fast magnetosonic speed c_f,x is mostly determined by the pressure flux $\widehat{p}/{\widehat{ℳ}}_{\text{son}}^2$ in Eq. (13). Since both slow incompressible flows and fast magnetosonic waves are permitted in the limit of ${\widehat{ℳ}}_{\text{son}}\to 0$, the system of fully compressible MHD equations is stiff.

2.3 Magnetohydrostatic solutions

Magnetohydrostatic stratifications are a special class of solutions to the MHD system with gravitational source terms (see Eqs. (1)–(4)), where all the time derivatives are zero and the velocity is zero everywhere. Under these conditions, the distribution of density, pressure and magnetic field is given by the magnetohydrostatic equation $$\nabla \cdot \left[\left(p+p_{\text{B}}\right)I-B\otimes B\right]=\rho g.$$(17)

Any set (ρ, p, B) that solves Eq. (17) is called a “magnetohydrostatic solution”. Equation (17) is undetermined, so a whole continuum of magnetohydrostatic solutions exists.

The stratification of stars is very well described by MHSE over a large fraction of their lifetime. Large deviations from MHSE are only expected in the late evolutionary stages of massive stars, in atmospheric layers and in stellar winds. Even though their stratification is continuously perturbed by a whole variety of physical processes over fast dynamical timescales, the amplitude of such perturbations remains small, and the overall structure of the star can be considered to be in MHSE. Significant changes to the stratification only happen over the much longer thermal and nuclear timescales (Kippenhahn et al. 2013).

2.4 The solenoidal constraint

The magnetic field satisfies the solenoidal constraint $$\nabla \cdot B=0.$$(18)

This constraint has its origin in Maxwell’s equation and states that physically no magnetic monopoles can exist. Solutions to Eqs. (1)–(4) automatically satisfy this condition at all times if the initial field obeys the constraint. This can easily be illustrated by rewriting Eq. (4) into the equivalent form $$\frac{\partial B}{\partial t}+\nabla \times ℰ=0,$$(19)

where ℰ = −V × B is the electromotive force. Applying the divergence to Eq. (19) results in $$\frac{\partial \left(\nabla \cdot B\right)}{\partial t}=0.$$(20)

3 Spatial discretization

The system of partial differential equations (PDEs) described in Sect. 2 takes the general conservative form $$\frac{\partial U}{\partial t}+\frac{\partial F\left(U\right)}{\partial x}+\frac{\partial G\left(U\right)}{\partial y}+\frac{\partial H\left(U\right)}{\partial z}=S\left(U\right),$$(21)

with the respective vector of conservative variables U, physical fluxes F, G, H and source term S. In SLH, Eq. (21) is solved numerically using the finite-volume method (LeVeque 2002; Toro 2009), which is briefly summarized in the next section.

3.1 Finite-volume discretization

In a first step, the physical system is mapped on a 3D Cartesian grid⁸ divided into N_x × N_y × N_z cells, whose spatial extent is given by [x_L, x_R] × [y_L,y_R] x [z_L, z_R]. Each cell in the computational domain is defined by the set of indices (i, j, k), and its volume is given by the product of the spatial resolution elements along each axis:Θ_i,j,k = ΔxΔyΔz. Any quantity located at the center of the cell refers to the same indices, while quantities located at the cell boundaries are denoted by sets of indices like (i + 1/2, j, k), which in this case defines the interface between cells (i, j, k) and (i + 1, j,k).

Integrating Eq. (21) over the cell volume leads to $$\begin{array}{lll}\frac{\partial {\widehat{U}}_{i,j,k}}{\partial t}\hfill & =\hfill & -\frac{1}{\text{Δ}x}\left({\widehat{F}}_{i+1/2,j,k}-{\widehat{F}}_{i-1/2,j,k}\right)\hfill \\ \hfill & \hfill & -\frac{1}{\text{Δ}y}\left({\widehat{G}}_{i,j+1/2,k}-{\widehat{G}}_{i,j-1/2,k}\right)\hfill \\ \hfill & \hfill & -\frac{1}{\text{Δ}z}\left({\widehat{H}}_{i,j,k+1/2}-{\widehat{H}}_{i,j,k-1/2}\right)\hfill \\ \hfill & \hfill & +{\widehat{S}}_{i,j,k},\hfill \end{array}$$(22)

where ${\widehat{U}}_{i,j,k}$ is the volume-averaged vector of conserved quantities $${\widehat{U}}_{i,j,k}=\frac{1}{{\text{Θ}}_{i,j,k}}{\displaystyle {\int }_{\text{Θ,}i,j,k}U\text{dΘ.}}$$(23)

The same procedure applies to Ŝ_i,j,k, while the surface-averaged fluxes are defined as $${\widehat{F}}_{i+1/2,j,k}=\frac{1}{A_{i+1/2,j,k}}{\displaystyle {\int }_{Ai+1/2,j,k}F\cdot \widehat{n}\text{d}A,}$$(24)

where A_i+1/2,j,k is the area of the interface (i + 1/2, j, k) and $\widehat{n}$ is the normal to the surface pointing outward from the cell.

The right-hand side of Eq. (22) can be discretized in space if suitable numerical values for ${\widehat{F}}_{i+1/2,j,k}$ and Ŝ_i,j,k are provided. For the latter, a typical choice consists in substituting the volume-averaged quantity with its value in the center of the cell, which is accurate to second order: $${\widehat{S}}_{i,j,k}\simeq S_{i,j,k}.$$(25)

The computation of numerical fluxes, in contrast, needs more care, and upwind techniques must be used to achieve stability. The resulting system of ordinary differential equations (ODEs) is then discretized in time using the methods described in Sect. 4.

3.2 Numerical flux function

In order to get a proper estimate of the fluxes ${\widehat{F}}_{i+1/2,j,k}$, we use the Godunov method (Godunov & Bohachevsky 1959). First, a pair of left and right states $U_{i+1/2,j,k}^L$, $U_{i+1/2,j,k}^R$ is reconstructed⁹ (through 1D sweeping) to the center of each cell boundary, starting from the cell-centered states ${\widehat{U}}_{i,j,k}$. These states define a 1D Riemann problem, which can then be solved (either exactly or approximately) to provide the value of a flux function $ℱ\left(U_{i+1/2,j,k}^L,U_{i+1/2,j,k}^R\right)$. The surface-averaged flux is then approximated (to second-order accuracy) as $${\widehat{F}}_{i+1/2,j,k}\simeq ℱ\left(U_{i+1/2,j,k}^L{U}_{i+1/2,j,k}^R\right).$$(26)

Many MHD Riemann solvers used nowadays are designed to work in supersonic regimes. In order to achieve numerical stability, such solvers need to add upwind numerical diffusion terms to the physical fluxes¹⁰, which smear out any discontinuity present in the flow on a timescale comparable to the cell crossing time of the shock. The choice of these terms depends on the specific approximate Riemann solver used. In particular, the diffusion term associated with the pressure flux (see Eq. (2)) usually takes the form (see, e.g., Einfeldt et al. 1991; Cargo & Gallice 1997; Miyoshi & Kusano 2005) $$D_{\text{p}}^x\propto -\overline{\rho }{\overline{c}}_{\text{f}}\left(V_x^R+V_x^L\right),$$(27)

where $\overline{\rho }$ and ${\overline{c}}_{\text{f}}$ are suitable averages of the density and the fast magnetosonic speed at the cell interface. However, in low-Mach regimes, discontinuities in the flow are only transported by the linearly degenerate entropy and Alfvén waves. These modes propagate with small speeds (see Sect. 2.2), and by the time they cross one cell in the computational grid they are strongly dissipated by the action of the numerical term in Eq. (27). This effect can also be explained by noticing that the pressure-diffusion coefficient scales as O(1/ℳ_fast,x), so that it overwhelms the physical flux proportional to O(1) at low sonic Mach numbers. To remove this excessive dissipation, we use the low-dissipation HLLD solver (LHLLD) of Minoshima & Miyoshi (2021). This is a variation of the original five-wave HLLD solver (see Fig. 2) of Miyoshi & Kusano (2005). LHLLD introduces a Mach-dependent parameter ϕ ∝ ℳ_fast,x in the intermediate state of the total pressure p_T = p + p_B: $$\begin{array}{lll}{p}_{\text{T}}^{*}\hfill & =\hfill & \frac{\left(S_{\text{R}}-V_x^{\text{R}}\right){\rho }_{\text{R}}{p}_{\text{T}}^{\text{L}}-\left(S_{\text{L}}-V_x^{\text{L}}\right){\rho }_{\text{L}}{p}_{\text{T}}^{\text{R}}}{\left(S_{\text{R}}-V_x^{\text{R}}\right){\rho }_{\text{R}}-\left(S_{\text{L}}-V_x^{\text{L}}\right){\rho }_{\text{L}}}\hfill \\ \hfill & \hfill & +\varphi \frac{{\rho }_{\text{L}}{\rho }_{\text{R}}\left(S_{\text{R}}-V_x^{\text{R}}\right)\left(S_{\text{L}}-V_x^{\text{L}}\right)\left(V_x^{\text{R}}-V_x^{\text{L}}\right)}{\left(S_{\text{R}}-V_x^{\text{R}}\right){\rho }_{\text{R}}-\left(S_{\text{L}}-V_x^{\text{L}}\right){\rho }_{\text{L}}}.\hfill \end{array}$$(28)

In this context, S_L and S_R are conservative estimates of the speeds λ_1,7. In SLH they are evaluated as $$\begin{array}{lll}{S}_{\text{L}}\hfill & =\hfill & \min \left(V_x^{\text{L}},V_x^{\text{R}}\right)-\max \left(c_{\text{f},x}^{\text{L}},c_{\text{f},x}^{\text{R}}\right),\hfill \\ S_{\text{R}}\hfill & =\hfill & \max \left(V_x^{\text{L}},V_x^{\text{R}}\right)+\max \left(c_{\text{f},x}^{\text{L}},c_{\text{f},x}^{\text{R}}\right).\hfill \end{array}$$(29)

The low-Mach fix ϕ is computed according to the following formulas: $$\begin{array}{lll}{c}_{\text{u}}^1\hfill & =\hfill & {\left[\frac{1}{2}\left(\frac{{\left|B_{\text{L}}\right|}^2}{{\rho }_{\text{L}}}+{\left|V_{\text{L}}\right|}^2+\sqrt{{\left(\frac{{\left|B_{\text{L}}\right|}^2}{{\rho }_{\text{L}}}+{\left|V_{\text{L}}\right|}^2\right)}^2-4\frac{{\left|V_{\text{L}}\right|}^2{B}_x^2}{{\rho }_{\text{L}}}}\right)\right]}^{\frac{1}{2}},\hfill \\ c_{\text{u}}^{\text{R}}\hfill & =\hfill & {\left[\frac{1}{2}\left(\frac{{\left|B_{\text{R}}\right|}^2}{{\rho }_{\text{R}}}+{\left|V_{\text{R}}\right|}^2+\sqrt{{\left(\frac{{\left|B_{\text{R}}\right|}^2}{{\rho }_{\text{R}}}+{\left|V_{\text{R}}\right|}^2\right)}^2-4\frac{{\left|V_{\text{R}}\right|}^2{B}_x^2}{{\rho }_{\text{R}}}}\right)\right]}^{\frac{1}{2}},\hfill \\ \chi \hfill & =\hfill & \max \left\{\frac{c_{\text{u}}^{\text{L}}}{c_{\text{f,}x}^{\text{L}}},\frac{c_{\text{u}}^{\text{R}}}{c_{\text{f,}x}^{\text{R}}}\right\},\hfill \\ \varphi \hfill & =\hfill & \chi \left(2-\chi \right).\hfill \end{array}$$(30)

Since the fast magnetosonic wave speeds $c_{\text{f},x}^{\text{L,R}}$ and consequently also S_LR scale as O(1/ℳ_fast,x), the second term in Eq. (28) would scale as O(1/ℳ_fast,x) if ϕ = 1, as in the original formulation of the HLLD solver. As previously described, this would lead to excessive numerical dissipation for small values of ℳ_fast,x. Instead, by computing ϕ according to Eq. (30), the dissipation term becomes independent of the fast magnetosonic Mach number, since ϕ∝ ℳ_fast,x. This modification does not affect the other properties of the HLLD solver, such as preserving positivity of density and internal energy (Miyoshi & Kusano 2005; Minoshima & Miyoshi 2021). We note that the combined diffusion coefficient in Eq. (28) has a residual scaling O(1/ℳ_Alf), which would still introduce too much dissipation in very sub-Alfvén regimes. However, these are far from our main astrophysical applications.

Fig. 2 Wave structure of HLLD-type solvers: only the fast magnetosonic (SL, SR), Alfvén $\left(S_{\text{L}}^{*},S_{\text{R}}^{*}\right)$, and entropy (SM) waves are considered in the computation of the states across the Riemann fan. The slow magnetosonic waves are discarded.

3.3 Well-balancing method

As already noted in Sects. 1 and 3.1, hyperbolic fluxes and gravitational source terms are discretized with different methods. As a consequence, Godunov-type schemes do not automatically preserve magnetohydrostatic solutions on a discrete grid exactly. Therefore, whenever a stratification needs to be enforced to be in MHSE on the computational grid, we use the deviation well-balancing method (Berberich et al. 2021; Edelmann et al. 2021). The main ingredient of this method is an a priori known target state Ũ that is a magnetohydrostatic solution to Eq. (21), $$\frac{\partial F\left(\tilde{U}\right)}{\partial x}+\frac{\partial G\left(\tilde{U}\right)}{\partial y}+\frac{\partial H\left(\tilde{U}\right)}{\partial z}=S\left(\tilde{U}\right),$$(31)

with Ṽ = 0. Subtracting Eq. (31) from the original balance law in Eq. (21) yields a system of PDEs for the deviations from the target solution ΔU = U - Ũ: $$\begin{array}{lll}\frac{\partial \left(\text{Δ}U\right)}{\partial t}\hfill & +\hfill & \left(\frac{\partial F\left(\tilde{U}+\text{Δ}U\right)}{\partial x}-\frac{\partial F\left(\tilde{U}\right)}{\partial x}\right)\hfill \\ \hfill & +\hfill & \left(\frac{\partial G\left(\tilde{U}+\text{Δ}U\right)}{\partial y}-\frac{\partial G\left(\tilde{U}\right)}{\partial y}\right)\hfill \\ \hfill & +\hfill & \left(\frac{\partial H\left(\tilde{U}+\text{Δ}U\right)}{\partial z}-\frac{\partial H\left(\tilde{U}\right)}{\partial z}\right)=S\left(\tilde{U}+\text{Δ}U\right)-S\left(\tilde{U}\right).\hfill \end{array}$$(32)

Now, to obtain a well-balanced method, Eq. (32) is discretized according to the finite-volume method described in Sect. 3.1, which leads to the semi-discrete form $$\begin{array}{lll}\frac{\partial {\left(\text{Δ}U\right)}_{i,j,k}}{\partial t}\hfill & =\hfill & -\frac{1}{\text{Δ}x}\left(F_{i+1/2,j,k}^{\text{dev}}-F_{i-1/2,j,k}^{\text{dev}}\right)\hfill \\ \hfill & \hfill & -\frac{1}{\text{Δ}y}\left(G_{i,j+1/2,k}^{\text{dev}}-G_{i,j-1/2,k}^{\text{dev}}\right)\hfill \\ \hfill & \hfill & -\frac{1}{\text{Δ}z}\left(H_{i,j,k+1/2}^{\text{dev}}-H_{i,j,k-1/2}^{\text{dev}}\right)\hfill \\ \hfill & \hfill & +S_{i,j,k}^{\text{dev}}.\hfill \end{array}$$(33)

In this formulation, the deviation fluxes and source terms are defined by $$F_{i+1/2,j,k}^{\text{dev}}={\widehat{F}}_{i+1/2,j,k}-F\left({\tilde{U}}_{i+1/2,j,k}\right),$$(34) $$G_{i,j+1/2,k}^{\text{dev}}={\widehat{G}}_{i,j+1/2,k}-G\left({\tilde{U}}_{i,j+1/2,k}\right),$$(35) $$H_{i,j,k+1/2}^{\text{dev}}={\widehat{H}}_{i,j,k+1/2}-H\left({\tilde{U}}_{i,j,k+1/2}\right),$$(36) $$S_{i,j,k}^{\text{dev}}={\widehat{S}}_{i,j,k}-S\left({\tilde{U}}_{i,j,k}\right),$$(37)

where ${\widehat{F}}_{i+1/2,j,k}$ is computed according to Eq. (26) in the states $$U_{i+1/2,j,k}^{\text{L,R}}={\tilde{U}}_{i+1/2,j,k}+\text{Δ}{U}_{i+1/2,j,k}^{\text{L,R}},$$(38)

while F (Ũ_i+1/2,j,k) corresponds to the physical fluxes in Eq. (21) evaluated in the target solution at the cell boundary. The deviations ΔU_i,j,k, rather than the states U_i,j,k, are reconstructed to the boundary of the cell¹¹. This guarantees that magnetohydrostatic solutions are preserved on the discrete grid, since in that case ΔU_i,j,k = 0, which leads to $${\widehat{F}}_{i+1/2,j,k}=ℱ\left({\tilde{U}}_{i+1/2,j,k,}{\tilde{U}}_{i+1/2,j,k}\right)=F\left({\tilde{U}}_{i+1/2,j,k}\right).$$(39)

Thus, the resulting method is well-balanced. Moreover, by removing the numerical errors arising from the magnetohydrostatic stratification, this method allows low-Mach flows to be simulated in stratified setups, which only cause small deviations from the MHSE state and would be completely dominated by spurious flows otherwise.

3.4 Constrained transport method

The divergence-free constraint described in Sect. 2.4 is not automatically satisfied if the induction equation is solved with Godunov-type schemes. As a result, magnetic monopoles are created locally at each time step and tend to accumulate, as they cannot be transported away by any of the MHD waves. If not properly treated, these artifacts can accelerate the flow along the magnetic field lines, generate wrong field topologies, and ultimately lead to severe stability problems (Brackbill & Barnes 1980).

Different strategies have been presented in the literature to cure this problem (for a review of these methods, see Tóth 2000). Among these, the eight-wave formulation (Powell 1997; Powell et al. 1999) modifies the MHD equations by including additional source terms that are proportional to ∇ · B. The modified system has an additional nonzero eigenvalue λ₈ = V_x, which transports jumps in the normal component (to the cell interface) of the magnetic field, so numerical monopoles are advected with the flow and do not accumulate over time.

Other solutions rely on divergence cleaning schemes (Dedner et al. 2002), where the divergence constraint is coupled to the MHD system using a generalized Lagrangian multiplier, ψ. This allows numerical monopoles to be transported with the maximum available speed on the grid and divergence errors to be damped at the same time.

One downside of both the eight-wave formulation and divergence cleaning is that they are not conservative and they cannot enforce any discretization of ∇ · B. to zero. Furthermore, these methods are most effective when open boundaries are used, so that the magnetic monopoles can leave the domain. However, this is rarely the case for simulations of stellar interiors, where impermeable boundaries are often used to avoid a significant mass loss from the system.

Constrained transport methods based on a staggered formulation, instead, conserve the magnetic flux through the boundaries of each cell and force one particular discretization of ∇ · B. to remain zero within round-off errors (Evans & Hawley 1988; Dai & Woodward 1998; Balsara & Spicer 1999; Tóth 2000; Londrillo & del Zanna 2004; Gardiner & Stone 2008; Mignone & Del Zanna 2021). Although a conservative scheme cannot guarantee that the discretized Lorentz force is orthogonal to the magnetic field lines in each cell of the computational grid (Tóth 2000), the magnitude of the parallel component of the force acting on the fluid is much smaller than in other methods.

The key point of staggered constrained transport methods is to compute the surface integral of Eq. (19) over cell boundaries using Stokes’s theorem, which leads to the finite-area equation¹² $$\begin{array}{lll}\frac{\partial {\widehat{B}}_{x,i+1/2,j,k}}{\partial t}\hfill & =\hfill & -\frac{1}{\text{Δ}y}\left({\widehat{ℰ}}_{z,i+1/2,j+1/2,k}-{\widehat{ℰ}}_{z,i+1/2,j-1/2,k}\right)\hfill \\ \hfill & \hfill & +\frac{1}{\text{Δ}z}\left({\widehat{ℰ}}_{y,i+1/2,j,k+1/2}-{\widehat{ℰ}}_{y,i+1/2,j,k-1/2}\right).\hfill \end{array}$$(40)

Here, ${\widehat{B}}_{x,i+1/2,j,k}$ is the surface-averaged magnetic field component normal to the cell boundary $${\widehat{B}}_{x,i+1/2,j,k}=\frac{1}{A_{i+1/2,j,k}}{\displaystyle {\int }_{A_{i+1/2,j,k}}B\cdot \widehat{n}\text{d}A,}$$(41)

while the line-averaged electromotive force is defined as $${\widehat{ℰ}}_{z,i+1/2,j+1/2,k}=\frac{1}{{\text{Δ}}_Z{}_{{}_{i+1/2,j+1/2,k}}}{\displaystyle {\int }_{\text{Δ}{Z}_{i+1/2,j+1/2,k}}{ℰ}_z{\text{d}}_z}.$$(42)

Analogous formulas can be derived for the other components of the magnetic field and the electromotive force.

In order to solve Eq. (40) numerically, a proper estimate for ${\widehat{B}}_{x,i+1/2,j,k}$ and ${\widehat{ℰ}}_{z,i+1/2,j+1/2,k}$ must be provided. For the former, we approximate the surface-averaged quantity with its value at the center of the cell boundary, $${\widehat{B}}_{x,i+1/2,j,k}\simeq B_{x,i+1/2,j,k},$$(43)

which is accurate to second order. In contrast to the standard finite-volume approach, here the magnetic field component normal to the interface is stored at cell boundaries, while the line-averaged electromotive force is evaluated at cell edges. Thus, the operation is performed on a staggered grid. Since the parallel magnetic field still needs to be reconstructed to compute the flux function, its value at cell-center locations is estimated as a simple arithmetic average between the neighboring cell interfaces: $$\begin{array}{lll}{\widehat{B}}_{x,i,j,k}\hfill & =\hfill & \frac{1}{2}\left({\widehat{B}}_{x,i-1/2,j,k}+{\widehat{B}}_{x,i+1/2,j,k}\right),\hfill \\ {\widehat{B}}_{y,i,j,k}\hfill & =\hfill & \frac{1}{2}\left({\widehat{B}}_{y,i,j-1/2,k}+{\widehat{B}}_{y,i,j+1/2,k}\right),\hfill \\ {\widehat{B}}_{z,i,j,k}\hfill & =\hfill & \frac{1}{2}\left({\widehat{B}}_{z,i,j,k-1/2}+{\widehat{B}}_{z,i,j,k+1/2}\right).\hfill \end{array}$$(44)

To compute the line-averaged electromotive force in Eq. (42), in SLH we use the CT-contact algorithm of Gardiner & Stone (2005). In this method, the electric field at cell edges is computed as a simple arithmetic average of the four neighboring face-centered electromotive force components, with the addition of a diffusion term that helps removing spurious oscillations when the magnetic field is advected. For instance, ${\widehat{ℰ}}_{z,i+1/2,j+1/2,k}$

is approximated to second-order accuracy by $$\begin{array}{lll}{\widehat{ℰ}}_{z,i+1/2,j+1/2,k}\hfill & \simeq \hfill & \frac{1}{4}({\overline{ℰ}}_{z,i+1/2,j,k}+{\overline{ℰ}}_{z,i+1/2,j+1,k}\hfill \\ \hfill & \hfill & +{\overline{ℰ}}_{z,i,j+1/2,k}+{\overline{ℰ}}_{z,i+1,j+1/2,k})\hfill \\ \hfill & \hfill & +\frac{\text{Δ}y}{8}\left\{{\left(\frac{\partial {ℰ}_z}{\partial y}\right)}_{i+1/2,j+1/4,k}-{\left(\frac{\partial {ℰ}_z}{\partial y}\right)}_{i+1/2,j+3/4,k}\right\}\hfill \\ \hfill & \hfill & +\frac{\text{Δ}x}{8}\left\{{\left(\frac{\partial {ℰ}_z}{\partial x}\right)}_{i+1/4,j+1/2,k}-{\left(\frac{\partial {ℰ}_z}{\partial x}\right)}_{i+3/4,j+1/2,k}\right\}\hfill \end{array}$$(45)

where ${\overline{ℰ}}_z$ can be computed from the solution to the Riemann problem in Eq. (26). The calculation for the x- and y-component is again analogous. The upwind diffusion term enters in the derivatives of the electromotive force in Eq. (45), which are obtained according to the sign s_i+1/2,j,,k of the entropy (contact) waves at the cell interfaces: $$\begin{array}{lll}{\left(\frac{\partial {ℰ}_z}{\partial y}\right)}_{i+1/2,j+1/4,k}\hfill & =\hfill & \frac{1+s_{i+1/2,j,k}}{2}\left(\frac{{\overline{ℰ}}_{z,i,j+1/2,k}-{ℰ}_{z,i,j,k}^{\text{cc}}}{\text{Δ}y/2}\right)\hfill \\ \hfill & \hfill & +\frac{1-s_{i+1/2,j,k}}{2}\left(\frac{{\overline{ℰ}}_{z,i+1,j+1/2,k}-{ℰ}_{z,i+1,j,k}^{\text{cc}}}{\text{Δ}y/2}\right).\hfill \end{array}$$(46)

Here ${ℰ}_{z,i,j,k}^{\text{cc}}={\left(-V_{i,j,k}\times B_{i,j,k}\right)}_z$ represents the z-component of the cell-centered electromotive force. The discretization of the line-averaged electromotive force leads to a semi-discrete form of Eq. (40) that can be integrated numerically in time. Any time-stepper that solves the resulting system of ODEs can keep the cell-volume average of ∇ · B, $$\begin{array}{lll}{\left(\nabla \cdot B\right)}_{i,j,k}\hfill & =\hfill & \frac{{\widehat{B}}_{x,i+1/2,j,k}-{\widehat{B}}_{x,i-1/2,j,k}}{\text{Δ}x}\hfill \\ \hfill & \hfill & +\frac{{\widehat{B}}_{y,i,j+1/2,k}-{\widehat{B}}_{y,i,j-1/2,k}}{\text{Δ}y}\hfill \\ \hfill & \hfill & +\frac{{\widehat{B}}_{z,i,j,k+1/2}-{\widehat{B}}_{z,i,j,k-1/2}}{\text{Δ}z},\hfill \end{array}$$(47)

within rounding errors.

4 Time integration algorithm

The CFL constraint in time-explicit marching schemes restricts the time step to the crossing time of the fastest wave resulting from the underlying PDEs over a grid cell. In low-Mach-number flows, the fast magnetosonic wave speeds become very large, so that the time step needs to be reduced accordingly. Thus, simulating the evolution of slow fluid motions and Alfvén waves becomes expensive. In these regimes, implicit methods, in which the time step is not limited by stability conditions but only by the desired accuracy, represent an attractive alternative. When using such methods, the time step should be restricted to the shortest advection and Alfvén crossing time over one grid cell. If the sonic Mach number is small enough, the possibility of larger chosen time steps then outweighs the disadvantage of higher computational costs for a single time step by using the implicit solver.

As outlined in Sect. 1, we split the induction equation (see Eq. (4)) from the continuity, momentum and energy equations (see Eqs. (1)–(3)), based on the approach described by Fuchs et al. (2009). This allows different spatial and temporal discretizations to be used depending on the problem at hand. In regimes of low Mach numbers and high β values, the stiffness is mostly generated by the pressure flux $\propto 1/{\widehat{ℳ}}_{\text{son}}^2$ in the momentum equation, while the nondimensional form of the induction equation does not depend on the Mach number of the flow (see Sect. 2.2). This suggests that implicit time discretization only needs to be applied to the subset of continuity, momentum and energy equations, whereas the induction equation can be solved with explicit time-steppers. These two updates can be combined to second-order accuracy with Strang splitting (Strang 1968): $$U^{n+1}=I^{\left(\frac{1}{2}\text{Δ}t\right)}{\mathscr{H}}^{\left(\text{Δ}t\right)}{I}^{\left(\frac{1}{2}\text{Δ}t\right)}{U}^n.$$(48)

Here, I represents a linear operator that updates only the magnetic field with an explicit marching scheme, while the nonlinear operator ℋ updates density, momentum and total energy (including source terms) using an implicit stepper. In each sub-step of Strang splitting, the discretization of the fluxes, source terms, and electromotive force is performed according to the methods described in Sect. 3.

In SLH, several implicit time-steppers can be used to solve the semi-discrete form of Eq. (22), such as first-order backward-Euler, higher-order ESDIRK schemes, and Crank-Nicolson. The resulting nonlinear system of equations is solved iteratively with a root-finding Raphson–Newton algorithm, which relies on the analytic formulation of the flux-Jacobian. Iterative linear solvers (such as BiCGSTAB(l), GMRES, and Multigrid) are used in combination with preconditioning techniques to solve each substep of the nonlinear solver¹³. In contrast, the semi-discrete form of the induction equation (see Eq. (40)) is solved with the time-explicit SSP-RK2 method of Shu & Osher (1988).

Numerical experiments performed with the proposed implicit-explicit Strang splitting (IESS) approach suggest that the maximum time step allowed for stability is approximately determined by $$\text{Δ}t=\underset{\text{Ω}=\left(i,j,k\right)}{min}\left\{\frac{\text{Δ}x}{\left|V_{x,\text{Ω}}\right|+c_{\text{a},x,\text{Ω}}},\frac{\text{Δ}y}{\left|V_{y,\text{Ω}}\right|+c_{\text{a},y,\text{Ω}}},\frac{\text{Δ}z}{\left|V_{z,\text{Ω}}\right|+c_{\text{a},z,\text{Ω}}}\right\},$$(49)

so that the propagation of fluid motions and Alfvén waves is well resolved in time. This time step is approximately 1/ℳ_son larger than that allowed by the conventional CFL condition if the plasma-β is high, which considerably reduces the computational effort when simulating low-Mach-number flows. The price one has to pay is that the propagation of fast magnetosonic waves is not well resolved in time. Another advantage of IESS is that it can easily be implemented within the framework of the SLH code, which already had fully implicit time integration capabilities to solve the compressible Euler equations.

A single step of the described time-marching scheme can be summarized in the following way. First, Δt is obtained from Eq. (49) given Bⁿ, ρⁿ, ρVⁿ, and $\rho E_{\varphi }^n$. If gravity is not present, eϕ does not appear in Eq. (5).

Second, SSP-RK2 and CT-contact are used to solve the induction equation over the first half of the time step, Δt/2. This results in an intermediate solution for the magnetic field, B^n+1/2.

Third, this intermediate solution, B^n+1/2, is used to solve the continuity, momentum, and energy equations over the full time step, Δt, with an implicit time-stepper. If gravity is present, then the well-balancing method described in Sect. 3.3 can be used. Any other source term is also considered in this step. This allows the solution for density, momentum, and energy to be obtained at the next step, ρⁿ⁺¹, ρVⁿ⁺¹, and $\rho E_{\varphi }^{n+1}$.

Fourth, B^n+1/2, ρⁿ⁺¹, ρVⁿ⁺¹, and $\rho E_{\varphi }^{n+1}$ are used to solve the induction equation over ΔT/2. This yields the magnetic field at the final step Bⁿ⁺¹.

The proposed MHD scheme is extremely modular, so different time-steppers, spatial reconstruction schemes and approximate Riemann solvers can be used in each sub-step of the algorithm, and well-balancing can be switched off if required. For instance, in addition to LHLLD, the original five-wave HLLD solver of Miyoshi & Kusano (2005) is also implemented in SLH. The performance and accuracy of both Riemann solvers are checked in some of the numerical experiments described in Sect. 5. Finally, in case slightly subsonic or transonic regimes need to be modeled with SLH, a fully un-split SSP-RK2 explicit time-stepper can be used.

5 Numerical tests

In order to assess the accuracy and performance of the newly implemented MHD algorithm, we have to rely on numerical experiments. Since the main purpose of the scheme is to be able to simulate MHD flows at low sonic Mach numbers in strong stratifications, we decide not to show the typical tests commonly run by other MHD codes. These usually include shock-tubes, supersonic vortices and magnetic blasts, which, however, are designed to test the shock-capturing capabilities of a numerical scheme. Instead, we ran a series of verification benchmarks that are more suited for testing the low-Mach properties of an MHD code.

As a first test, we solved the homogeneous MHD equations in three different cases (i, ii, and iii). In order to check the convergence and scaling of the methods for the whole MHD wave family, we performed a 1D linear analysis (i). The scaling was also checked against the advection of a stable MHD vortex in a wide range of Mach numbers (ii). Such a setup is particularly important as it resembles the typical vortex structures present in magneto-convection. The simulations were also run in fully-explicit mode using SSP-RK2, which allows the speed-up of IESS to be quantified as a function of the Mach number.

The ability of accurately evolving shear instabilities is fundamental in the context of simulations of turbulence as they generate additional vorticity, which leads to the cascade of energy. For this reason, we ran simulations of a magnetized Kelvin–Helmholtz instability (iii). We followed the growth and evolution of the instability in a resolution study from low-Mach to slightly subsonic regimes. A comparison between the HLLD and LHLLD solvers was performed to show the advantage of using low-dissipation fluxes over conventional methods in regimes of low Mach numbers.

Then we considered two setups in which gravity is present (iv, v). To check the entropy-conservation properties of the scheme based on the deviation well-balancing method, we modeled the rise of a parcel of fluid with higher entropy content than the (isentropic) background stratification, that is, a “hot bubble” (iv). By changing the magnitude of the entropy perturbation, we simulated different rise velocities of the bubble, down to Mach numbers of ℳ_son ~ 10⁻⁴. To quantify the magnitude of the numerical errors generated by an unbalanced stratification, we also simulated the rise of the bubble at ℳ_son ~ 10⁻² without well-balancing.

Finally, we simulated a fully 3D small-scale dynamo (SSD) amplification in a star-like environment at moderate grid resolutions (ν). By changing the rate at which energy is injected in the system, we simulated progressively slower flows down to ℳ_son ~ 10⁻³.

For all of the following tests, an ideal gas EoS was used with γ = 5/3 except when specified otherwise. Within the framework of the IESS time-marching scheme described in Sect. 4, ESDIRK23 was chosen to treat the implicit part of the algorithm. This guarantees second-order accuracy in time. The time step in Eq. (49) was reduced by 20% to get a more conservative stability criterion. Finally, unlimited linear reconstruction, which is second-order accurate in space, was applied to primitive variables. Overall, the proposed scheme is (globally) second-order accurate.

5.1 Linear analysis

In this test, we followed the propagation of linear modes for all the MHD waves (see Sect. 2) as a way to get quantitative estimates of diffusion errors and to check the scaling of the numerical scheme. The setup is based on Stone et al. (2008), which we modified by considering much larger values of the gas pressure to increase the fast magnetosonic speed relative to the Alfvén and entropy wave speeds. Such a stiff system is characteristic of low-Mach flows in high-β environments (see Sect. 2.2).

The homogeneous MHD equations were solved on a periodic ID Cartesian grid divided into N cells, with the spatial domain ranging from 0 to L = 1. For the ith wave, the solution at t = 0 was obtained by perturbing a uniform medium w₀ = (1,0,0,0,10³/γ, 1, $\sqrt{2}$, 1/2) with δw = AR_i(w₀) sin(2πx/L), where w is the vector of primitive variables (ρ, V_x, V_y, V_z, p, Β_x, B_y, Β_z), A = 10⁻⁴ is the amplitude of the perturbation, x is the spatial coordinate and R_i is the ith column of the right-eigenvector matrix (see Appendix A.3 in Stone et al. 2008). The chosen values for w₀ are such that c_f,x ≃ 40.85, c_a,x = 1 and c_s,x ≃ 0.99. For the entropy mode, we set V_x = 1, so the sonic Mach number of the wave is ℳ_son ≃ 0.032. The simulations were run for one crossing time defined as t_c = L/λ_i, where λ_i is the wave speed. The L₁ error was then computed for each primitive variable w_k as $$L_{1,i}\left(w_k\right)=\frac{1}{N}{\displaystyle \sum _j\left|w_{k,j}\left(t=t_{\text{c}}\right)-w_{k,j}\left(t=0\right)\right|,}$$(50)

and the global error associated with the ith wave was then computed as $${\delta }_i=\frac{1}{A}\sqrt{{\displaystyle \sum _k{\left(\frac{L_{1,i}\left(w_k\right)}{\xi \left(w_{0,k}\right)}\right)}^2}},$$(51)

with $$\xi \left(w_{0,k}\right)=\{\begin{array}{ll}{w}_{0,k},\hfill & \text{for\hspace{0.17em}}{w}_{0,k}>0,\hfill \\ 1,\hfill & \text{otherwise}.\hfill \end{array}$$(52)

In Eq. (50), j is the spatial index. The tolerance of the Raphson-Newton algorithm was set to 10⁻¹⁰, so that the errors computed using Eq. (51) were not dominated by the finite convergence of the nonlinear solver.

Figure 3 shows the change of the global error as a function of N (ranging from 10 to 100) for all the seven MHD waves. Leftward and rightward propagating waves have identical errors. Since the MHD scheme relies on second-order methods to treat both the spatial and the temporal parts, the scheme converges with second-order accuracy for all waves except the fast magnetosonic waves, which are characterized by much larger errors. This is expected since IESS allows the MHD equations to be integrated over much longer time steps than the CFL constraint (see Sect. 4). As a consequence, the propagation of fast magnetosonic waves is not properly resolved in time and discretization errors strongly deteriorate the numerical solution. This effect is further quantified in Fig. 4, where we show the global error associated with the left-going fast magnetosonic wave as a function of the grid resolution for different values of the gas pressure p, such that c_f,x/c_a,x = 2, 3, 6, 12, 27, 59, 129. Overall, the errors tend to decrease for smaller values of c_f,x/c_a,x, since the time step in Eq. (49) gets closer to the CFL time step and fast magnetosonic waves are progressively better resolved. Moreover, the simulations run with c_f,x/c_a,x ≤ 27 converge with second-order accuracy on the grids considered in this study.

Fig. 3 Global error as a function of the grid resolution for the seven MHD waves after one crossing time, tc. The dashed black line represents the second-order scaling.

Fig. 4 Global error as a function of the grid resolution for the left-going fast magnetosonic wave. Different colors are associated with different values of cf,x/ca,x. The dashed black line is the second-order scaling.

5.2 Balsara vortex

In the previous section, we demonstrated that the current MHD scheme is capable of simulating linear waves with the expected (second-order) scaling with respect to resolution on 1D grids. In order to check the scaling in 2D and to test the low-Mach capabilities of the scheme, we considered the MHD vortex first described by Balsara (2004). This is an exact stationary solution of the ideal 2D homogeneous MHD equations, in which the distribution of the centrifugal acceleration, magnetic tension, gas and magnetic pressure gradients is such that the vortex is stable. The spatial domain is (x, y) ∈ [−5, 5] × [−5, 5], and we used 64 × 64 grid cells with periodic boundaries in both directions. The initial conditions are given by $$\begin{array}{lll}\left(V_x,V_y\right)\hfill & =\hfill & \tilde{V}{e}^{\frac{1-r^2}{2}}\left(-y,x\right),\hfill \\ \left(B_x,B_y\right)\hfill & =\hfill & \tilde{B}{e}^{\frac{1-r^2}{2}}\left(-y,x\right),\hfill \\ \hfill p& =\hfill & 1+\left[\frac{{\tilde{B}}^2}{2}\left(1-r^2\right)-\frac{{\tilde{V}}^2}{2}\right]e^{1-r^2}\hfill \\ \hfill \rho & =\hfill & 1,\hfill \end{array}$$(53)

with r² = x² + y². Ṽ is the maximum rotational velocity of the vortex and $\tilde{B}$ sets the value of the maximum Alfvén speed on the grid. To make this problem numerically more challenging, the vortex is advected along the diagonal of the computational grid, with ∣V_adv∣ = Ṽ. The vortex is evolved for one advective crossing time $t_{\text{adv}}=10\sqrt{2}/\tilde{V}$, after which it returns to the initial position. In this time interval, the vortex rotates 2.25 times.

We ran the grid of models $$\begin{array}{lll}\left(\tilde{V}\right)\times \left({\beta }_{\text{K}}\right)\hfill & =\hfill & \left({10}^{-5},{10}^{-4},{10}^{-3},{10}^{-2},{10}^{-1}\right)\hfill \\ \hfill & \hfill & \times \left({10}^{-2},{10}^{-1},1,{10}^1,{10}^2\right),\hfill \end{array}$$(54)

with ${\beta }_{\text{K}}={\tilde{B}}^2/{\tilde{V}}^2$ being the ratio of the magnetic to the rotational kinetic energy, which is constant across the domain. Given this choice of parameters, the initial maximum Mach number ℳ_son ranges from 1.55 × 10⁻⁵ to 1.55 × 10⁻¹, so this parameter study covers both low Mach numbers and slightly subsonic regimes, in both weakly and strongly magnetized fluids.

Figure 5 shows the magnetic energy distribution after one advective crossing time t_adv. Numerical dissipation converts a fraction of kinetic and magnetic energy into internal energy, but the shape of the vortex is well preserved in all runs. The dissipation rate is virtually independent of ℳ_son. In contrast, dissipation of magnetic energy depends on the value of β_K. As already pointed out in Sect. 3.2, the pressure-diffusion coefficient in LHLLD has a residual scaling O(1/ℳ_Alf). A larger value of β_K corresponds to lower ℳ_Alf, which then increases the magnitude of the numerical dissipation. The velocity field is progressively more diffused out and becomes less efficient in sustaining the magnetic field through induction against numerical resistivity.

To check the convergence of the scheme in 2D, we ran a vortex with V = 10⁻³ (corresponding to max(ℳ_son)_t=0 = 1.55 × 10⁻³) and β_K = 1 at different resolutions¹⁴. At the end of the simulation, the L₁ error was computed for each primitive variable w_k as $$L_1\left(w_k\right)=\frac{1}{N^2}{\displaystyle \sum _{i,j}\left|w_{k,i,j}\left(t=t_{\text{adv}}\right)-w_{k,i,j}\left(t=0\right)\right|,}$$(55)

where i, j are the spatial indices. Figure 6 shows the convergence of the L₁ error for different grids from N = 32 up to N = 512 cells per dimension. Convergence is second order for all primitive variables.

To compare the amount of numerical dissipation introduced by a standard and a low-Mach MHD flux function, we reran this last set of simulations with HLLD. In Fig. 7, we show the final rotational kinetic energy distribution obtained with the two methods: $$E_{\text{R}}=\frac{1}{2}\rho \left[{\left(V_x-\tilde{V}/\sqrt{2}\right)}^2+{\left(V_y-\tilde{V}/\sqrt{2}\right)}^2\right].$$(56)

At low resolution, HLLD considerably stretches the vortex and a large fraction of kinetic energy is dissipated into internal energy. In contrast, simulations run with LHLLD show mild dissipation and dispersion errors are only visible at the lowest resolutions. All simulations converge with increasing resolution, but the kinetic energy conservation in the vortex simulated with HLLD is still two orders of magnitude worse than that obtained with LHLLD at the highest resolution considered in this study.

As explained in Sect. 4, one advantage of IESS is that the MHD equations can be integrated on time steps longer than that allowed by the CFL condition without sacrificing stability. However, a single step of the proposed scheme is much more expensive than a single step of a more standard time-explicit marching scheme, as a large nonlinear system has to be solved iteratively with a Raphson–Newton method. Because of these competing effects, we expect the IESS scheme to be more efficient than an explicit time-stepper below a certain Mach number. To determine this threshold, we ran sets of simulations with the parameters $$\begin{array}{lll}\left(\tilde{V}\right)\times \left({\beta }_{\text{K}}\right)\hfill & =\hfill & \left({10}^{-4},{10}^{-3},{10}^{-2},{10}^{-1}\right)\hfill \\ \hfill & \hfill & \times \left({10}^{-1},1,{10}^1\right),\hfill \end{array}$$(57)

using both IESS and the explicit SSP-RK2¹⁵ on 40 × 40 grid cells. Every other sub-step of the Godunov method (like the spatial reconstruction, the LHLLD flux function and constrained transport) remained unchanged, so the only difference was in the time discretization. At the end of each simulation, the ratio of the wall-clock times WCT_SSP–RK2/WCT_IESS was taken as a measure of the relative efficiency between the marching schemes¹⁶. The results are shown in Fig. 8. As expected, the speed-up of IESS increases as the Mach number of the vortex is decreased. The simulations with β_K = 10 are slower than the other cases, as the larger Alfvén speed considerably reduces the time step estimate in Eq. (49), while no significant difference is seen between β_K = 0.1 and β_K = 1.0. IESS overtakes SSP-RK2 at max(ℳ_son)_t=0 ≃ 4 × 10⁻² for β_K = (0.1, 1) and max(ℳ_son)_t=0 ≃ 2 × 10⁻² for β_K = 10. At max(ℳ_son)_t=0 = 10⁻³, IESS is 10–20 times faster than SSP-RK2. This justifies the implementation efforts of apartially implicit time discretization algorithm for modeling slow flows.

Fig. 5 Magnetic energy distribution of the Balsara vortex after one advective crossing time, tadv, normalized by the maximum magnetic energy at t = 0. The ratio of the magnetic to the (rotational) kinetic energy of the vortex is varied along the y-axis (in descending order), while the initial maximum rotational velocity, Ṽ, varies along the x-axis. The inset in each subplot shows the ratio of the final to the initial magnetic energy. The vortex run with Ṽ = 10−1 and βK = 102 (bottom-right corner) has a maximum Mach number max(ℳson)t=0 = 1.65 × 10−1. In that system, the gas pressure drops in the regions around the center of the vortex to balance the large magnetic and centrifugal forces, which ultimately decreases the1 sound speed where the velocity is maximum.

5.3 Magnetized Kelvin–Helmholtz instability

For the following test, we ran MHD simulations of a Kelvin–Helmholtz instability. This is the primary instability that arises when there is a velocity shear within a continuous fluid, and it is the main source of vorticity that leads to the energy cascade in 3D turbulent flows. An accurate representation of this process is therefore a fundamental requirement for any numerical scheme to be used for simulating magneto-convection. We considered a 2D domain with (x, y) ∈ [0, 2] × [−0.5, 0.5], mapped on a 2N × N grid. The horizontal velocity profile is given by $$V_x={ℳ}_x\left[1-2\eta \left(x\right)\right],$$(58)

with $$\eta \left(x\right)=\{\begin{array}{ll}\frac{1}{2}\left\{1+\sin \left[16\pi \left(y+0.25\right)\right]\right\},\hfill & \text{for\hspace{0.17em}}y>-\frac{9}{32}\text{and\hspace{0.17em}}y<-\frac{7}{32},\hfill \\ 1,\hfill & \text{for\hspace{0.17em}}y\ge -\frac{7}{32}\text{and\hspace{0.17em}}y\le \frac{7}{32},\hfill \\ \frac{1}{2}\left\{1-\sin \left[16\pi \left(y-0.25\right)\right]\right\},\hfill & \text{for\hspace{0.17em}}y>\frac{7}{32}\text{and\hspace{0.17em}}y<\frac{9}{32},\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}$$(59)

The parameter ℳ_x is the maximum sonic Mach number of the horizontal flow, and p = 1 and ρ = γ, so initially the adiabatic sound speed c is 1 everywhere. In this test, γ = 1.4. The magnetic field at t = 0 is uniform and horizontal (B_x = 0.1 ℳ_x), and the minimum Alfvén Mach number ℳ_Alf is 11.82 for all values of ℳ_x.

It is well known that magnetic fields aligned with the shear flow have a stabilizing effect because they exert a restoring force on the perturbed interface (Chandrasekhar 1961). With a too strong field, the instability may reach saturation when the flow is still essentially laminar or it may be suppressed completely. Instead, weak magnetic stresses do not considerably affect the initial growth of the instability, so the flow can develop the typical vortex structures present in the pure hydrodynamic case. This leads to a much more complex evolution in the nonlinear phase (Frank et al. 1996). For this setup, nearly laminar flows are expected only when min(ℳ_Alf)_t=0 ≲ 1.1, as shown in Fig. A.1.

The instability is started by adding a perturbation to the y-velocity component in the initial state, V_y = 0.1ℳ_x sin(2πx) (see Fig. 9). The initial conditions are periodic in both directions. The evolution of the Kelvin–Helmholtz instability was studied for a wide range of Mach numbers and grid resolutions: $$\begin{array}{lll}\left({ℳ}_x\right)\times \left(N\right)\hfill & =\hfill & \left({10}^{-4},{10}^{-3},{10}^{-2},{10}^{-1}\right)\hfill \\ \hfill & \hfill & \times \left(32,64,128,256,512,1024\right).\hfill \end{array}$$(60)

The final time reached by each simulation was set according to the initial amplitude of the shear flow (t_max = 4.8/ℳ_x). The chosen initial conditions are such that the interface across the shear flow is smooth and resolved, which leads to convergent results at least in the early stages of the evolution of the flow. As in the previous test, we compared the results obtained with both the HLLD and LHLLD solvers.

Figures 10 and 11 show the time evolution of the y-direction kinetic energy $E_{\text{K},y}={\displaystyle {\sum }_{ij}{\left(\rho V_y^2\right)}_{ij}/2}$ and the total magnetic energy E_M = Σ_ij |B_ij|²/2 for all the simulations considered in this study. As in the previous problem, i, j are the spatial indices. Because of stretching and wrapping of the field lines within the vortices, the magnetic energy slowly increases with time at the expense of the kinetic energy content of the flow. After the primary rolls reach the top and bottom boundaries (t/t_max ≃ 0.25), E_K,y saturates due to the periodicity of the grid and starts to decrease. The secondary vortices keep winding up the magnetic field lines until Lorentz forces start to feedback on the velocity field, breaking down these inner structures. The two original shear interfaces get closer to each other (see Fig. A.3) until a strong numerical reconnection event happens at t/t_max ≃ 0.45, which violently decouples the primary rolls and causes a secondary peak in E_K,y at t/t_max ≃ 0.5. After this time, other reconnection events break down the flow into smaller structures, and both the magnetic and the kinetic energy are slowly dissipated away by the action of numerical resistivity and viscosity.

Since in this case we solved the ideal MHD equations, there is no characteristic scale on which magnetic and kinetic energy are dissipated into heat, so numerical effects play a significant role on progressively smaller scales at higher resolution. Thus, the amplification and dissipation of magnetic energy hardly converge for the resolutions considered in this study. The initial growth of E_K,y, in contrast, is not much influenced by the initial weak field, and it is mostly determined by the strength of the shear flows and the width of the shear interface, which is resolved. As a consequence, E_K,y converges until the major numerical reconnection event affects the velocity field. As shown in Fig. 10, the HLLD solver requires more resolution to reach convergence as the setup is run at progressively lower sonic Mach numbers. Eventually, the Mach-dependent pressure-diffusion coefficient in Eq. (27) completely dominates the evolution of the flow and deteriorates the numerical solution. For this reason, at ℳ_x = 10⁻⁴ we were able to successfully run with HLLD only the 64 × 32 and 128 × 64 grids, while for higher resolutions the nonlinear solver failed to converge.

The effects of numerical dissipation are also shown in Fig. 12, where the distributions of the sonic Mach number obtained with HLLD and LHLLD are compared at fixed resolution (128 × 64 cells) for different values of ℳ_x at t/t_max = 1/6. While in moderately subsonic regimes the large-scale structures in the flow are qualitatively similar, for lower Mach numbers HLLD introduces progressively more dissipation and the instability is eventually halted. When LHLLD is used instead, the morphology of the flow seems to be independent of the Mach number.

Finally, we performed a quantitative convergence study by computing the L₁ error associated with E_K,y at t/t_max = 1/6. At this time, the first rolls have developed to considerable vertical wavelengths (see Fig. 12) so that the instability has already entered the nonlinear regime, and the flow is expected to converge as shown in Fig. 10. The L₁ error was computed against a reference solution, which was taken from the highest grid resolution runs considered in this test (N = 1024) using the LHLLD solver. All simulations (including the reference solutions) were down-sampled to a 64 × 32 grid, so that the errors could directly be computed for different resolutions. This analysis was repeated for different values of ℳ_x using both HLLD and LHLLD. The results are shown in Fig. 13. The errors are rescaled by ${ℳ}_x^2$. so that curves corresponding to different sonic Mach numbers lie on the same scale. Overall, the convergence is second-order with N for all simulations. LHLLD provides almost identical (rescaled) errors at given resolution in different regimes of Mach numbers. This is expected because the numerical dissipation introduced by this solver does not depend on ℳ_son, thanks to the low-Mach fix in Eq. (30). Instead, the errors computed for the HLLD runs show a clear dependence on the sonic Mach number, and the errors get larger for slower flows. In particular, at ℳ_x = 0.1, HLLD needs approximately 1.2 times the resolution of LHLLD to achieve the same accuracy, which justifies the use of HLLD in this regime of Mach numbers. Instead, when ℳ_x = 10⁻² and ℳ_x = 10⁻³, HLLD needs respectively twice or four times the resolution to be as accurate as the low-dissipation flux, which increases the amount of computing time by 8 or 64. Thus, the use of a low-Mach approximate Riemann solver becomes indispensable for providing accurate results in regimes of low sonic Mach numbers with moderate grid resolutions, which would be unfeasible with more standard solvers.

Fig. 6 Convergence of the L1 error in the Balsara vortex for each primitive variable as a function of resolution. For these simulations, Ṽ = 10−3 and βK = 1. The dashed black line is the second-order scaling.

Fig. 7 Distribution of the rotational kinetic energy (normalized by the maximum initial value) of the Balsara vortex after one advective time, tadv, at Ṽ = 10−3 and βK = 1. The top panels show the vortices obtained with the HLLD flux function as a function of resolution, while the plots in the bottom panels are obtained with LHLLD. The insets show the fraction of rotational kinetic energy that has been dissipated by the end of the simulation: ${\left(E_{\text{R},t=t_{\text{adv}}}\right)}_{\text{tot}}$/(ER,t=0)tot − 1.

Fig. 8 Ratio of the wall clock times obtained with SSP-RK2 and IESS as a function of initial maximum sonic Mach number of the Balsara vortex for different magnetic to (rotational) kinetic energy ratios. The dashed black line is drawn to represent the same relative efficiency.

5.4 Hot bubble

Flows in deep stellar convection zones are usually characterized by the presence of slow parcels of fluid that move in a stratification that is unstable against convection. In the absence of volume heating and cooling processes, these packets of fluid preserve their entropy content until they mix with the surroundings. Therefore, a numerical scheme designed to simulate such flows should have good entropy-conservation properties. However, entropy conservation is hard to achieve if the density, temperature and pressure stratifications span several orders of magnitude and if the flows are very slow, since their entropy content would only be slightly higher or lower than the adiabatic surroundings¹⁷. Under these conditions, discretization errors caused by an imperfect balance of the background MHSE stratification can dominate the dynamics and deteriorate the numerical solution. As shown in Sect. 3.3, the magnitude of such errors can be drastically reduced by using well-balancing techniques.

In this section, we check the entropy-conservation properties of the MHD scheme implemented in SLH by running simulations of the hot bubble setup described by Edelmann et al. (2021), where a bubble of higher entropy content with respect to the surroundings buoyantly rises in an adiabatic stratification. The physical domain is mapped on a 2D Cartesian grid (N_x = 2/3 × N_y), and the background stratification is in MHSE. Boundary conditions are periodic everywhere and the gravitational acceleration takes the form $$g_y\left(y\right)=g_0\sin \left(k_{y}y\right),$$(61)

where g₀ = −1.09904373 × 10⁵ cm s⁻², k_y = 2π/L_y, y is the vertical spatial coordinate, and L_y is the vertical extent of the grid. The value of g₀ is set such that the ratio of the maximum to the minimum gas pressure¹⁸ p(y) is 100, which corresponds to 4.6 pressure scale heights. The entropy profile inside the bubble is given by $$A=A_0\left\{1+{\left(\frac{\text{Δ}A}{A}\right)}_{t=0}\cos {\left(\frac{\pi }{2}\frac{r}{r_0}\right)}^2\right\},$$(62)

where A₀ is background entropy, r₀ is the radius of the bubble, r is the distance from the center of the bubble and (ΔA/A)_t=0 is the initial entropy perturbation. The density is $$\rho \left(y\right)={\left(\frac{p\left(y\right)}{A}\right)}^{1/\gamma },$$(63)

so that the (initial) buoyant acceleration of the bubble is proportional to the entropy perturbation, $$a_{\text{b}}=\frac{\text{Δ}\rho }{\rho }{g}_y\propto {\left(\frac{\text{Δ}A}{A}\right)}_{t=0}.$$(64)

We ran the models for the set of parameters $$\begin{array}{lll}{\left(\right(\frac{\Delta A}{A}\left)\right)}_{t=0}\times \left(\right(N_y\left)\right)\hfill & =\hfill & \left(\right({10}^{-7},{10}^{-5},{10}^{-3},{10}^{-1}\left)\right)\hfill \\ \hfill & \hfill & \times \left(\right(96,192,384,768\left)\right),\hfill \end{array}$$(65)

and we set the maximum time such that in each run the bubble raised approximately the same distance l. This allowed different regimes of sonic Mach numbers to be simulated, as the velocity, V, reached by the bubble over a length, l, scales as $$V\propto {\left(a_{\text{b}}l\right)}^{1/2}.$$(66)

This ultimately leads to the relation $${ℳ}_{\text{son}}\propto {\left(\frac{\text{Δ}A}{A}\right)}_{t=0}^{1/2}.$$(67)

A uniform horizontal magnetic field was added to the system, and its strength was rescaled depending on the entropy perturbation, $$B_x=B_0{\left(\frac{\text{Δ}A}{A}\right)}_{t=0}^{1/2}.$$(68)

This ensures that the relative magnitude of magnetic stresses compared to the ram pressure of the bubble remains the same for all simulations, and that the morphology of the flow is unaltered. B₀ = 47.3 was chosen such that the final Alfvén Mach number at the position of largest entropy is in the range ℳ_Alf ≃ 2–3 depending on the grid resolution. Thus, magnetic fields are dynamically important but not strong enough to suppress buoyancy.

In Fig. B.1, we show the final entropy excess for all the simulations run in the parameter study. The center of the bubble accelerates faster than other regions as it is the point with maximum entropy, and the acceleration profile across the bubble leads to the development of shear at its outer edges. As the bubble rises in the stratification, the magnetic field lines are stretched into thin tubes, which locally amplifies the magnetic energy (see Fig. B.2). The amount of amplification depends on the numerical resistivity and so on resolution. In contrast to the pure hydrodynamic case studied by Edelmann et al. (2021), here the presence of a magnetic field suppresses the formation of vortices at the sides of the bubble. Overall, the entropy content of the bubble is well preserved even on the coarsest grid, but some negative entropy fluctuations are present at the very top of the bubble. These negative fluctuations are numerical artifacts. In fact, the entropy fluctuations may locally increase as a fraction of magnetic and kinetic energy is dissipated into internal energy, but they cannot become negative physically. These artifacts do not depend on the entropy perturbation, and they are limited to a very narrow region in the spatial domain that tends to shrink as the resolution is increased. All models converge upon grid refinement.

According to Eq. (67), the sonic Mach number of the bubble is expected to scale as the square root of the initial entropy perturbation. Any deviation from this relation, which has been obtained on the basis of physical arguments, can be due to difficulties in modeling slow flows in a stratified setup and the build-up of significant numerical errors. In Fig. 14, we show this scaling for the coarsest grid resolution. All data points overlap with the theoretical curve, and the minimum Mach number ℳ_son achieved in this parameter study is 3.32 × 10⁻⁴ (see also Fig. B.3). The ratio of the rising velocity of the bubble to the Alfvén speed (in the point of maximum entropy) does not depend on the amplitude of the entropy perturbation. Since the initial magnetic field is proportional to ${\left(\text{Δ}A/A\right)}_{t=0}^{1/2}$, the amount of amplification due to induction only depends on the velocity of the bubble V and the timescale over which magnetic induction operates (∝ 1/V).

Finally, to quantify the strength of the spurious flows that are expected to arise if the stratification is left unbalanced, in Fig. 15 we show a comparison between simulations obtained with and without deviation well-balancing, where the vertical resolution N_y ranges from 96 to 768. For this comparison, we fixed ${\left(\text{Δ}A/A\right)}_{t=0}^{1/2}$ = 10⁻³ such that the maximum sonic Mach number of the bubble is approximately 3 × 10⁻². The unbalanced simulations develop large entropy fluctuations, both negative and positive, which strongly deteriorate the numerical solution. As the grid is refined, the simulations tend to converge, but wide regions of negative entropy fluctuations are still present even on the finest grid. Thus, this test demonstrates that well-balancing techniques are fundamental to correctly simulate the evolution of small entropy perturbations in steep isentropic stratifications and to reduce the effects of numerical errors when using moderately coarse grids.