Issue 
A&A
Volume 668, December 2022



Article Number  A143  
Number of page(s)  29  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202244665  
Published online  19 December 2022 
A finitevolume scheme for modeling compressible magnetohydrodynamic flows at low Mach numbers in stellar interiors
^{1}
Heidelberger Institut für Theoretische Studien,
SchlossWolfsbrunnenweg 35,
69118
Heidelberg, Germany
email: giovanni.leidi@hits.org
^{2}
Zentrum für Astronomie der Universität Heidelberg, Astronomisches RechenInstitut,
Mönchhofstr. 12–14,
69120
Heidelberg, Germany
^{3}
Department of Mathematics, Würzburg University,
EmilFischerStr. 40,
97074
Würzburg, Germany
^{4}
Computer, Computational and Statistical Sciences (CCS) Division and Center for Theoretical Astrophysics (CTA), Los Alamos National Laboratory,
Los Alamos, NM
87545, USA
^{5}
Zentrum für Astronomie der Universität Heidelberg, Institut für Theoretische Astrophysik,
Philosophenweg 12,
69120
Heidelberg, Germany
Received:
2
August
2022
Accepted:
18
September
2022
Fully compressible magnetohydrodynamic (MHD) simulations are a fundamental tool for investigating the role of dynamo amplification in the generation of magnetic fields in deep convective layers of stars. The flows that arise in such environments are characterized by low (sonic) Mach numbers (ℳ_{son} ≲ 10^{−2}). In these regimes, conventional MHD codes typically show excessive dissipation and tend to be inefficient as the Courant–Friedrichs–Lewy (CFL) constraint on the time step becomes too strict. In this work we present a new method for efficiently simulating MHD flows at low Mach numbers in a spacedependent gravitational potential while still retaining all effects of compressibility. The proposed scheme is implemented in the finitevolume SEVENLEAGUE HYDRO (SLH) code, and it makes use of a lowMach version of the fivewave Harten–Lax–van Leer discontinuities (HLLD) solver to reduce numerical dissipation, an implicit–explicit time discretization technique based on Strang splitting to overcome the overly strict CFL constraint, and a wellbalancing method that dramatically reduces the magnitude of spatial discretization errors in strongly stratified setups. The solenoidal constraint on the magnetic field is enforced by using a constrained transport method on a staggered grid. We carry out five verification tests, including the simulation of a smallscale dynamo in a starlike environment at ℳ_{son} ~ 10^{−3}. We demonstrate that the proposed scheme can be used to accurately simulate compressible MHD flows in regimes of low Mach numbers and strongly stratified setups even with moderately coarse grids.
Key words: magnetohydrodynamics (MHD) / methods: numerical
© G. Leidi et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1 Introduction
The interplay between turbulent convection and shear is fundamental in understanding the role of small and largescale dynamo mechanisms in the generation of strong magnetic fields in stellar interiors. These processes can only be modeled selfconsistently 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 finitevolume discretization and Godunovlike 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^{−2} (Kupka & Muthsam 2017), where V is the flow speed and a is the adiabatic sound speed. It is well known that conventional finitevolume 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^{−2}. Second, explicit timesteppers 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 lowMachnumber 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 Godunovtype 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 lowMachnumber 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 lowMach regime (ℳ_{son} ≳ 10^{−2}), where explicit timesteppers 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 lowMach version of the fivewave 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), timeimplicit 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 semiimplicit (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)^{1}, 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 semiimplicit 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 timemarching 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 timestepper. 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 speedup when the Mach number of the flow is low. Since the update on the induction equation is performed in a separate step, the fluxJacobian in the timeimplicit 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 (CTcontact; 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 wellbalancing 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^{2}. Recently, Canivete Cuissa & Teyssier (2022) performed fully compressible simulations of stellar magnetoconvection at ℳ_{son} ~ 10^{−3} using a wellbalancing technique similar to the deviation method. However, their scheme relied on explicit timesteppers, 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 SEVENLEAGUE 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 (timeindependent) gravitational source term^{3}: (1) (2) (3) (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^{4}, I the unit tensor, g = (g_{x}, g_{y}, g_{z}) the gravitational acceleration, p the gas pressure and p_{B} = ∣B∣^{2}/2 the magnetic pressure. The total energy density ρE_{ϕ} is defined as (5)
where e_{int} and e_{ϕ} are the specific internal and gravitational energies^{5}.
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 (lefthand side of Eqs. (1)–(4)) reduced to one spatial dimension has seven nonzero eigenvalues^{6}, (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, (7) (8)
with the adiabatic sound speed a defined as (9)
As illustrated in Fig. 1, the MHD wave pattern has a fixed ordering: (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 (11)
while the directional slow and fast magnetosonic Mach numbers are given by (12)
Fig. 1 Wave structure of the MHD system. 
2.2 LowMach 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^{−2} 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^{7}, (13) (14) (15) (16)
Here, the different variables have been rescaled by some reference quantity representative of the physical system of interest: , , , , , . and are the characteristic sonic and Alfvén Mach numbers of the flow.
In the limit of , Eqs. (13)–(15) approach the incompressible regime (see Matthaeus & Brown 1988), in which the gas pressure is constant in space except for fluctuations ∝ . 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 , these waves travel at similar speeds. However, the incompressible limit is not the only solution to the compressible MHD equations at . 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 lowMach limit. If the plasmaβ is high, the large fast magnetosonic speed c_{f,x} is mostly determined by the pressure flux in Eq. (13). Since both slow incompressible flows and fast magnetosonic waves are permitted in the limit of , 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 (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 (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 (19)
where ℰ = −V × B is the electromotive force. Applying the divergence to Eq. (19) results in (20)
3 Spatial discretization
The system of partial differential equations (PDEs) described in Sect. 2 takes the general conservative form (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 finitevolume method (LeVeque 2002; Toro 2009), which is briefly summarized in the next section.
3.1 Finitevolume discretization
In a first step, the physical system is mapped on a 3D Cartesian grid^{8} 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 (22)
where is the volumeaveraged vector of conserved quantities (23)
The same procedure applies to Ŝ_{i,j,k}, while the surfaceaveraged fluxes are defined as (24)
where A_{i+1/2,j,k} is the area of the interface (i + 1/2, j, k) and is the normal to the surface pointing outward from the cell.
The righthand side of Eq. (22) can be discretized in space if suitable numerical values for and Ŝ_{i,j,k} are provided. For the latter, a typical choice consists in substituting the volumeaveraged quantity with its value in the center of the cell, which is accurate to second order: (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 , we use the Godunov method (Godunov & Bohachevsky 1959). First, a pair of left and right states , is reconstructed^{9} (through 1D sweeping) to the center of each cell boundary, starting from the cellcentered states . These states define a 1D Riemann problem, which can then be solved (either exactly or approximately) to provide the value of a flux function . The surfaceaveraged flux is then approximated (to secondorder accuracy) as (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^{10}, 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) (27)
where and are suitable averages of the density and the fast magnetosonic speed at the cell interface. However, in lowMach 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 pressurediffusion 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 lowdissipation HLLD solver (LHLLD) of Minoshima & Miyoshi (2021). This is a variation of the original fivewave HLLD solver (see Fig. 2) of Miyoshi & Kusano (2005). LHLLD introduces a Machdependent parameter ϕ ∝ ℳ_{fast,x} in the intermediate state of the total pressure p_{T} = p + p_{B}: (28)
In this context, S_{L} and S_{R} are conservative estimates of the speeds λ_{1,7}. In SLH they are evaluated as (29)
The lowMach fix ϕ is computed according to the following formulas: (30)
Since the fast magnetosonic wave speeds 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 subAlfvén regimes. However, these are far from our main astrophysical applications.
Fig. 2 Wave structure of HLLDtype solvers: only the fast magnetosonic (S_{L}, S_{R}), Alfvén , and entropy (S_{M}) waves are considered in the computation of the states across the Riemann fan. The slow magnetosonic waves are discarded. 
3.3 Wellbalancing method
As already noted in Sects. 1 and 3.1, hyperbolic fluxes and gravitational source terms are discretized with different methods. As a consequence, Godunovtype 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 wellbalancing 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), (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  Ũ: (32)
Now, to obtain a wellbalanced method, Eq. (32) is discretized according to the finitevolume method described in Sect. 3.1, which leads to the semidiscrete form (33)
In this formulation, the deviation fluxes and source terms are defined by (34) (35) (36) (37)
where is computed according to Eq. (26) in the states (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^{11}. This guarantees that magnetohydrostatic solutions are preserved on the discrete grid, since in that case ΔU_{i,j,k} = 0, which leads to (39)
Thus, the resulting method is wellbalanced. Moreover, by removing the numerical errors arising from the magnetohydrostatic stratification, this method allows lowMach 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 divergencefree constraint described in Sect. 2.4 is not automatically satisfied if the induction equation is solved with Godunovtype 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 eightwave 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 λ_{8} = 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 eightwave 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 roundoff 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 finitearea equation^{12} (40)
Here, is the surfaceaveraged magnetic field component normal to the cell boundary (41)
while the lineaveraged electromotive force is defined as (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 and must be provided. For the former, we approximate the surfaceaveraged quantity with its value at the center of the cell boundary, (43)
which is accurate to second order. In contrast to the standard finitevolume approach, here the magnetic field component normal to the interface is stored at cell boundaries, while the lineaveraged 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 cellcenter locations is estimated as a simple arithmetic average between the neighboring cell interfaces: (44)
To compute the lineaveraged electromotive force in Eq. (42), in SLH we use the CTcontact 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 facecentered electromotive force components, with the addition of a diffusion term that helps removing spurious oscillations when the magnetic field is advected. For instance,
is approximated to secondorder accuracy by (45)
where can be computed from the solution to the Riemann problem in Eq. (26). The calculation for the x and ycomponent 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: (46)
Here represents the zcomponent of the cellcentered electromotive force. The discretization of the lineaveraged electromotive force leads to a semidiscrete form of Eq. (40) that can be integrated numerically in time. Any timestepper that solves the resulting system of ODEs can keep the cellvolume average of ∇ · B, (47)
within rounding errors.
4 Time integration algorithm
The CFL constraint in timeexplicit marching schemes restricts the time step to the crossing time of the fastest wave resulting from the underlying PDEs over a grid cell. In lowMachnumber 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 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 timesteppers. These two updates can be combined to secondorder accuracy with Strang splitting (Strang 1968): (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 substep 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 timesteppers can be used to solve the semidiscrete form of Eq. (22), such as firstorder backwardEuler, higherorder ESDIRK schemes, and CrankNicolson. The resulting nonlinear system of equations is solved iteratively with a rootfinding Raphson–Newton algorithm, which relies on the analytic formulation of the fluxJacobian. 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^{13}. In contrast, the semidiscrete form of the induction equation (see Eq. (40)) is solved with the timeexplicit SSPRK2 method of Shu & Osher (1988).
Numerical experiments performed with the proposed implicitexplicit Strang splitting (IESS) approach suggest that the maximum time step allowed for stability is approximately determined by (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 lowMachnumber 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 timemarching scheme can be summarized in the following way. First, Δt is obtained from Eq. (49) given B^{n}, ρ^{n}, ρV^{n}, and . If gravity is not present, eϕ does not appear in Eq. (5).
Second, SSPRK2 and CTcontact 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 timestepper. If gravity is present, then the wellbalancing 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, ρ^{n+1}, ρV^{n+1}, and .
Fourth, B^{n+1/2}, ρ^{n+1}, ρV^{n+1}, and are used to solve the induction equation over ΔT/2. This yields the magnetic field at the final step B^{n+1}.
The proposed MHD scheme is extremely modular, so different timesteppers, spatial reconstruction schemes and approximate Riemann solvers can be used in each substep of the algorithm, and wellbalancing can be switched off if required. For instance, in addition to LHLLD, the original fivewave 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 unsplit SSPRK2 explicit timestepper 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 shocktubes, supersonic vortices and magnetic blasts, which, however, are designed to test the shockcapturing capabilities of a numerical scheme. Instead, we ran a series of verification benchmarks that are more suited for testing the lowMach 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 magnetoconvection. The simulations were also run in fullyexplicit mode using SSPRK2, which allows the speedup 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 lowMach to slightly subsonic regimes. A comparison between the HLLD and LHLLD solvers was performed to show the advantage of using lowdissipation 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 entropyconservation properties of the scheme based on the deviation wellbalancing 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^{−4}. To quantify the magnitude of the numerical errors generated by an unbalanced stratification, we also simulated the rise of the bubble at ℳ_{son} ~ 10^{−2} without wellbalancing.
Finally, we simulated a fully 3D smallscale dynamo (SSD) amplification in a starlike 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^{−3}.
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 timemarching scheme described in Sect. 4, ESDIRK23 was chosen to treat the implicit part of the algorithm. This guarantees secondorder 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 secondorder accurate in space, was applied to primitive variables. Overall, the proposed scheme is (globally) secondorder 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 lowMach 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_{0} = (1,0,0,0,10^{3}/γ, 1, , 1/2) with δw = AR_{i}(w_{0}) 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^{−4} is the amplitude of the perturbation, x is the spatial coordinate and R_{i} is the ith column of the righteigenvector matrix (see Appendix A.3 in Stone et al. 2008). The chosen values for w_{0} 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_{1} error was then computed for each primitive variable w_{k} as (50)
and the global error associated with the ith wave was then computed as (51)
In Eq. (50), j is the spatial index. The tolerance of the RaphsonNewton algorithm was set to 10^{−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 secondorder methods to treat both the spatial and the temporal parts, the scheme converges with secondorder 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 leftgoing 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 secondorder 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, t_{c}. The dashed black line represents the secondorder scaling. 
Fig. 4 Global error as a function of the grid resolution for the leftgoing fast magnetosonic wave. Different colors are associated with different values of c_{f,x}/c_{a,x}. The dashed black line is the secondorder 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 (secondorder) scaling with respect to resolution on 1D grids. In order to check the scaling in 2D and to test the lowMach 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 (53)
with r^{2} = x^{2} + y^{2}. Ṽ is the maximum rotational velocity of the vortex and 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 , after which it returns to the initial position. In this time interval, the vortex rotates 2.25 times.
We ran the grid of models (54)
with 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^{−5} to 1.55 × 10^{−1}, 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 pressurediffusion 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^{−3} (corresponding to max(ℳ_{son})_{t=0} = 1.55 × 10^{−3}) and β_{K} = 1 at different resolutions^{14}. At the end of the simulation, the L_{1} error was computed for each primitive variable w_{k} as (55)
where i, j are the spatial indices. Figure 6 shows the convergence of the L_{1} 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 lowMach 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: (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 timeexplicit 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 timestepper below a certain Mach number. To determine this threshold, we ran sets of simulations with the parameters (57)
using both IESS and the explicit SSPRK2^{15} on 40 × 40 grid cells. Every other substep 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 wallclock times WCT_{SSP–RK2}/WCT_{IESS} was taken as a measure of the relative efficiency between the marching schemes^{16}. The results are shown in Fig. 8. As expected, the speedup 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 SSPRK2 at max(ℳ_{son})_{t=0} ≃ 4 × 10^{−2} for β_{K} = (0.1, 1) and max(ℳ_{son})_{t=0} ≃ 2 × 10^{−2} for β_{K} = 10. At max(ℳ_{son})_{t=0} = 10^{−3}, IESS is 10–20 times faster than SSPRK2. 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, t_{adv}, 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 yaxis (in descending order), while the initial maximum rotational velocity, Ṽ, varies along the xaxis. The inset in each subplot shows the ratio of the final to the initial magnetic energy. The vortex run with Ṽ = 10^{−1} and β_{K} = 10^{2} (bottomright 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 magnetoconvection. 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 (58)
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 yvelocity 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: (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 ydirection kinetic energy and the total magnetic energy E_{M} = Σ_{ij} B_{ij}^{2}/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 Machdependent pressurediffusion coefficient in Eq. (27) completely dominates the evolution of the flow and deteriorates the numerical solution. For this reason, at ℳ_{x} = 10^{−4} 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 largescale 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_{1} 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_{1} 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 downsampled 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 . so that curves corresponding to different sonic Mach numbers lie on the same scale. Overall, the convergence is secondorder 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 lowMach 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^{−2} and ℳ_{x} = 10^{−3}, HLLD needs respectively twice or four times the resolution to be as accurate as the lowdissipation flux, which increases the amount of computing time by 8 or 64. Thus, the use of a lowMach 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 L_{1} 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 secondorder scaling. 
Fig. 7 Distribution of the rotational kinetic energy (normalized by the maximum initial value) of the Balsara vortex after one advective time, t_{adv}, 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: /(E_{R,t=0})_{tot} − 1. 
Fig. 8 Ratio of the wall clock times obtained with SSPRK2 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 entropyconservation 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^{17}. 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 wellbalancing techniques.
In this section, we check the entropyconservation 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 (61)
where g_{0} = −1.09904373 × 10^{5} cm s^{−2}, 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_{0} is set such that the ratio of the maximum to the minimum gas pressure^{18} p(y) is 100, which corresponds to 4.6 pressure scale heights. The entropy profile inside the bubble is given by (62)
where A_{0} is background entropy, r_{0} 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 (63)
so that the (initial) buoyant acceleration of the bubble is proportional to the entropy perturbation, (64)
We ran the models for the set of parameters (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 (66)
This ultimately leads to the relation (67)
A uniform horizontal magnetic field was added to the system, and its strength was rescaled depending on the entropy perturbation, (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_{0} = 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 buildup 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^{−4} (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 , 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 wellbalancing, where the vertical resolution N_{y} ranges from 96 to 768. For this comparison, we fixed = 10^{−3} such that the maximum sonic Mach number of the bubble is approximately 3 × 10^{−2}. 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 wellbalancing 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.
Fig. 9 Initial setups of V_{x} and V_{y} (here rescaled by ℳ_{x}) used for simulating the growth of the Kelvin–Helmholtz instability. 
Fig. 10 Time evolution of the ydirection kinetic energy rescaled by its initial value in the magnetized Kelvin–Helmholtz instability test problem. Each panel corresponds to a different initial Mach number, ℳ_{x}. Different colors are used for different grid resolutions (the 64 × 32 and 128 × 64 grids cells have been left out for clarity). Results obtained with the HLLD solver are represented by dotdashed lines, while solid lines are used for LHLLD. The solid black line in each panel is the reference solution. As explained in the text, the nonlinear solver does not converge when using HLLD at ℳ_{x} = 10^{−4} for N > 64. 
Fig. 12 Distribution of the sonic Mach number in the KelvinHelmholtz instability test at t/t_{max} =1/6 obtained with the HLLD (top panels) and the LHLLD (bottom panels) solvers on a 128 × 64 grid for different values of ℳ_{x}. All panels are rescaled by the corresponding value of M_{x}. 
Fig. 13 Convergence with resolution N of the L_{1} error associated with E_{K,y} rescaled by in the simulations of the Kelvin–Helmholtz instability. Different colors are used for different initial sonic Mach numbers ℳ_{x} using the LHLLD (solid lines) and HLLD (dotdashed lines) solvers. The dashed black line is the secondorder scaling. 
Fig. 14 Maximum sonic Mach number, minimum and maximum entropy fluctuations and Alfvén speed of the hot bubble as a function of the initial entropy perturbation obtained on a 64 × 96 grid. The black lines represent the physical scalings. 
5.5 Smallscale dynamo in a stratified setup
The previous tests demonstrated that the proposed MHD implementation can accurately simulate slow flows even in strongly stratified setups. As this numerical method will mostly be applied to simulate stellar interiors, it seems natural to test the scheme for dynamo amplification, which is the main cause for the generation of strong magnetic fields in a wide mass range of stars (see Brun & Browning 2017, and references therein). In this section we focus on simulations of SSDs, where the magnetic energy is amplified on scales comparable to or smaller than the scales at which turbulence is forced (Meneguzzi et al. 1981; Schekochihin et al. 2004, 2007; Brandenburg & Subramanian 2005; Iskakov et al. 2007), in contrast to largescale dynamos where most of the magnetic energy is at scales larger than the forcing scale (Brun et al. 2004; Käpylä et al. 2008; Charbonneau 2013; Augustson et al. 2016). Even though the efficiency of the dynamo amplification depends on many physical parameters, including the magnetic Prandtl number P_{m} = ν/η (Schekochihin et al. 2004, 2007; Pietarila Graham et al. 2010; Brandenburg 2011, 2014), here we do not perform a parameter study for P_{m}, since in the current MHD scheme the viscosity (ν) and resistivity (η) coefficients are not fixed, but are intrinsic to the underlying numerical methods, so they are not easy to constrain. Instead, we aim to check if it is possible to excite an SSD using the SLH code at low sonic Mach numbers.
We built the initial conditions based on the work of Andrassy et al. (2022), who performed a pure hydrodynamic study of a 3D convection zone with a stable layer on top of it, where the convective flows had a typical maximum Mach number ℳ_{son} ~ 0.1. The stratification of that model resembled oxygen shell burning in a massive star, even though some simplifications were adopted. Among these, an ideal gas EoS was used and effects of neutrino cooling were ignored. Here, we modified that setup even further by retaining only the convection zone and by decreasing the rate of energy injection to test our method in the lowMachnumber regime. Removing the stable layer greatly simplifies the problem at low Mach numbers, as the propagation of internal gravity waves does not need to be resolved. Since the wavelength of these modes becomes shorter for progressively slower convective flows (Sutherland 2010; Edelmann et al. 2021), high grid resolutions would be necessary to capture this process accurately at low Mach numbers, which would make the simulations very expensive.
For our experiment, we used N_{x} × N_{y} × N_{z} = 2N × N × 2N grid cells and the spatial domain (normalized by a characteristic length L = 4 × 10^{8} cm) is −1 ≤ x ≤ 1, 1 ≤ y ≤ 2, −1 ≤ z ≤ 1. Periodic boundaries were used in the horizontal directions, while reflecting boundaries were used in the vertical direction. The initial stratification is adiabatic and in MHSE, and it is given by the polytropic relation (69)
The stratification spans 2.2 pressure scale heights. The gravitational acceleration takes the form (70)
where g_{0} = 1.414870, and (71)
Moreover, in contrast to Andrassy et al. (2022), we are not interested in studying convective boundary mixing, so we only used a single species with mean molecular weight µ = 1.
Convection is driven by a constant in time heat source placed close to the bottom boundary, (72)
where = 3.795720 × 10^{−4} erg cm^{−3} s^{−1} and b is a nondimensional factor that allows the strength of the convective flows to be controlled through the scaling (see, e.g., Kippenhahn et al. 2013; Andrassy et al. 2020; Horst et al. 2021; Edelmann et al. 2021; Käpylä 2021) (73)
We considered the following grid of models: (74)
We ran the simulations at nominal luminosity (b = 1) with SSPRK2, where the convective flows are characterized by a maximum ℳ_{son} of 0.1. In this regime of sonic Mach numbers, IESS is less efficient than explicit timesteppers (see Fig. 8). In contrast, the cases with b = 10^{−3} and b = 10^{−6} were run with IESS, since according to Eq. (73), the maximum ℳ_{son} is 10^{−2} and 10^{−3}, respectively.
In order to initiate an SSD, a weak seed magnetic field was added to the system^{19}, (75)
where the dependence on b is such that the timescale on which the magnetic energy reaches saturation (expressed in units of convective turnovers) does not depend on the value of b. The convective turnover timescale was estimated as (76)
is averaged over time, and ,, represent the standard deviation of each velocity component computed over the whole domain. All simulations were run until t/τ_{conv} = 70, with τ_{conv} ≃ 5000, 500, 50 s for b = 10^{−6}, 10^{−3}, 1, respectively. The initial MHSE density stratification was perturbed to initiate convection^{20}. This perturbation leads to the development of buoyant structures that rise in the stratification until they hit the top boundary, after which they quickly become turbulent (see Fig. C.2). Convection fully develops after one convective turnover.
Figure 16 shows the time evolution of the kinetic and the magnetic energy for all the simulations run in the grid of models considered in this test. In the kinematic phase, the magnetic field is irrelevant to the dynamics, and it is amplified exponentially by the action of a dynamo process, with most of the magnetic energy distributed close to the resistive scale. As visible in the horizontal cuts shown in Fig. 17, the magnetic field distribution is characterized by smallscale structures with mixed polarity, while the velocity field is distributed on slightly larger scales, which suggests that P_{m} ≳ 1. The growth rate γ_{M} = d(lnE_{M})/dτ_{conv} increases with resolution (see also Fig. C.4), which is compatible with SSD amplification. In particular, we find that γ_{M} ∝ Δx^{−0.8}. The dependence of the growth rate on Δx is weaker than γ_{M} ∝ Δx^{−4/3}, which is typically observed in simulations of SSDs in solar and stellar convection zones (Pietarila Graham et al. 2010; Rempel 2014; Hotta et al. 2015; Riva & Steiner 2022; Canivete Cuissa & Teyssier 2022), and steeper than γ_{M} ∝ Δx^{−2/3}, which is predicted by the Kazantsev dynamo theory (Kazantsev 1968; Brandenburg & Subramanian 2005).
When the magnetic field becomes strong enough, the Lorentz force starts to influence the evolution of the turbulent flows, damping the velocity on the small scales (see the bottom panels in Fig. 17). A statistically steady state configuration is then reached where the magnetic energy achieves subequipartition values E_{M}/E_{K} ≃ 0.1−0.2. In all simulations, an SSD is successfully excited, and the mean value of the amplified magnetic energy increases with resolution (see Table 1). In fact, the size of the resistive scale is smaller on finer grids, which in turn increases the maximal stretching rate of the field lines and the SSD becomes more efficient. In contrast, no systematic difference is observed in the magnetic to the kinetic energy ratio in simulations run with the same resolution but different values of b. This is due to the fact that, thanks to the use of the LHLLD solver, the size of the viscous and resistive scales does not depend on the sonic Mach number of the flow, which is mostly determined by the chosen value of b. At given resolution, the SSD amplifies the magnetic field on the same spatial (resistive) scales, so the evolution of E_{M}/E_{K} becomes virtually independent of ℳ_{son} (and so of b) if the time is rescaled by the convective turnover τ_{conv}, except for statistical fluctuations caused by the chaotic nature of the turbulent flows.
Figure 18 shows the root mean square sonic Mach number averaged over 20 convective turnovers in the saturated regime as a function of b. As noted in Edelmann et al. (2021), numerical errors introduced by an unbalanced stratification can cause deviations from the scaling in Eq. (73). In this case, the use of deviation wellbalancing and LHLLD allows a good agreement between the computed Mach numbers and the scaling law to be reached, within three standard deviations. This proves that the convective flows are correctly simulated and are not dominated by numerical errors. The smallest (ℳ_{son})rms achieved in these runs is approximately 3.7 × 10^{−4}, which is close to what typically found in simulations of coreconvective stars (Augustson et al. 2016; Edelmann et al. 2019; Horst et al. 2020; Higl et al. 2021).
Finally, Fig. 19 shows the kinetic and magnetic energy spectra (taken in the midplane of the box) in the saturated stage for different values of b and grid resolutions. Both spectra have been rescaled by b^{2/3} to take into account the different energy contents of the flows achieved with different values of b. The kinetic energy spectra converge to the k^{−5/3} Kolmogorov law (Kolmogorov 1941) upon grid refinement, and the dissipation range shifts toward progressively larger wave numbers k. The magnetic energy distributions peak in the inertial range, as expected in SSD simulations, and on the large scales they show a shallower dependence on k than the Kazantsev isotropic dynamo theory, k^{3/2} (Kazantsev 1968). This can be explained by the fact that, in this setup, turbulence is not isotropic, and largescale anisotropic convective flows are present because of the steep stratification and the use of closed vertical boundaries (see Fig. C.2). Magnetic and kinetic energy achieve equipartition at the bottom of the inertial range, except for the simulations run on 64 × 32 × 64 grid cells, in which equipartition is reached only in the dissipation range. The maximum magnetic to kinetic energy ratio is around 3 in the dissipation range for the highest resolution considered in this study. Again, since the numerical diffusion of the MHD scheme is Machindependent, the shape and amplitude of the rescaled spectra do not depend on b on any of the three grids. Thus, this test indicates that SLH is capable of correctly simulating fully compressible magnetoconvection and SSDs in regimes of low sonic Mach numbers even with moderate grid resolutions.
Fig. 15 Final distribution of the entropy fluctuations of the hot bubble for (ΔA/A)_{t=0} = 10^{−3} at different grid resolutions. The entropy fluctuations are rescaled by (ΔA/A)_{t=0}. The top row shows the results obtained with deviation wellbalancing, whereas no wellbalancing method was used in the simulations shown in the bottom row. The insets show the minimum and maximum values of the entropy fluctuation in each panel. 
Fig. 16 Time evolution of the magnetic energy (sky blue) and kinetic energy (vermillion) in the simulations of the SSD for different grid resolutions (dotted: 64 × 32 × 64, dotdashed: 128 × 64 × 128, solid: 192 × 96 × 192). Each panel shows the results obtained with a specific value of b (from left to right: b = 10^{−6}, 10^{−3}, 1). The time is expressed in units of the convective turnover, while the magnetic and kinetic energy curves are rescaled by b^{2/3} to take into account the different energy contents of the flows. 
Fig. 17 Horizontal slices in the y = 1.5 plane taken in the kinematic (upper plots) and saturated (lower plots) regimes of the SSD with b = 10^{−6} on the 192 × 96 × 192 grid. The panels on the left show the vertical sonic Mach number ℳ_{son,y} = V_{y/a} multiplied by 10^{3}, while the plots on the right show the vertical magnetic field rescaled by the root mean square value across the plane. 
Time averages of E_{M}/b^{2/3} [10^{45} × erg] over the time interval 20 < t/τ_{conv} < 70 for the different resolutions, N, and boost factors, b, considered in this study.
Fig. 18 Root mean square of the sonic Mach number as a function of the driving luminosity b. The data points are averages computed in the time interval 20 < t/τ_{com} < 40. The error bars represent three standard deviations over the time series, while the dotted black line is the b^{1/3} scaling. 
Fig. 19 Kinetic (solid lines) and magnetic (dashed lines) energy spectra as a function of the wave number k obtained in an horizontal slice of the convective box at y = 1.5. All spectra are averaged over the time interval 20 < t/τ_{conv} < 40 and are rescaled by b^{2/3}. The three panels show the results for different numbers of grid cells (from top to bottom: 64 × 32 × 64, 128 × 64 × 128, 192 × 96 × 192). The dotted black lines are the Kolmogorov (k^{−5/3}) and Kazantsev (k^{3/2}) scalings. 
6 Summary and conclusions
In this work, we have presented a new finitevolume scheme for solving the fully compressible MHD equations with gravity in regimes of low sonic Mach numbers and highβ environments, which is suitable for simulating magnetoconvection and dynamo processes in deep layers of stars. This method relies on a lowdissipation MHD Riemann solver (LHLLD; Minoshima & Miyoshi 2021) to avoid the excessive numerical dissipation typical of highresolution, shockcapturing solvers as ℳ_{son} → 0.
The strict CFL condition on the time step is overcome by using an implicitexplicit time discretization algorithm, for which the induction equation is integrated using an explicit timestepper, while the rest of the MHD system is integrated implicitly. The solutions to the two subsets of equations are coupled through Strang splitting following the prescription of Fuchs et al. (2009). The combined marching scheme has a less restrictive condition on the time step, which is limited only by the fastest fluid and Alfvén speeds instead of the fast magnetosonic speed, leading to a considerable speedup in regimes of low sonic Mach numbers.
Whenever required, a magnetohydrostatic solution can be enforced on the discrete grid with the deviation wellbalancing method (Berberich et al. 2021; Edelmann et al. 2021). This technique leads to better entropyconservation properties of the numerical scheme, even in cases where the pressure and density stratifications span several orders of magnitude across the computational domain. Finally, the ∇ · Β = 0 condition is enforced using the CTcontact method (Gardiner & Stone 2005). This new scheme is implemented in the SLH code, and it has been tested in five numerical experiments.
First, we checked the global convergence of the methods by following the propagation of linear modes for all the MHD waves on a ID grid. This test proves that the scheme is globally secondorder accurate.
For the second test, we ran simulations of an MHD vortex advected along the diagonal of a square 2D grid. The characteristic velocities involved in the problem were varied such that the maximum Mach number, ℳ_{son}, ranged from 1.55 × 10^{−5} to 1.55 × 10^{−1}. This experiment showed that the MHD scheme also scales with secondorder accuracy on 2D grids and that the numerical dissipation is independent of ℳ_{son}, even though it becomes larger for lower ℳ_{Alf} However, we observed a considerable dissipation only when the magnetic energy of the vortex was 100 times its rotational kinetic energy. This regime is far from our main astrophysical applications. The dissipation of kinetic energy for both the LHLLD and the standard HLLD solvers has been quantified for different resolutions at a maximum Mach number ℳ_{son} = 1.55 × 10^{−3}. Even though all the results converged upon grid refinement, conservation of rotational kinetic energy was two orders of magnitude worse when using the HLLD flux instead of LHLLD. We also quantified the efficiency of IESS over a standard SSPRK2 as a function of ℳ_{son}. When the maximum sonic Mach number of the flow is below (2–5) × 10^{−2}, IESS becomes considerably more efficient than explicit timesteppers.
In the third experiment we considered the growth of a KelvinHelmholtz instability under the effects of a weak magnetic field parallel to the shear flow, which is known to generate more complex vortex structures than the case with a strong magnetic field (Frank et al. 1996). The secondorder convergence of the ydirection kinetic energy was checked using both the LHLLD and HLLD solvers for different shear velocities, such that the corresponding sonic Mach number ranged from 10^{−4} to 10^{−1}. Again, we observed that the amount of dissipation was virtually independent of ℳ_{son} for LHLLD, while the numerical solution obtained with HLLD was progressively more degraded as the Mach number was further decreased. This test showed that HLLD needed twice or four times the resolution to be as accurate as LHLLD when ℳ_{son} ~ 10^{−3} and ℳ_{son} ~ 10^{−2}, while it only needed 20% more resolution at ℳ_{son} ~ 10^{−1}.
In our fourth test we simulated the rise of a hot bubble in an adiabatic stratification in MHSE. The initial magnetic field was horizontal and uniform. Different entropy perturbations were considered to test the capabilities of the MHD scheme in modeling slow flows in steep stratifications (in this case, the vertical domain spanned 4.6 pressure scale heights). Overall, the entropy content of the bubble was very well preserved and all results converged upon grid refinement. Thanks to the deviation method, we were able to successfully simulate the rise of the bubble for entropy perturbations as low as 10^{−7}, leading to typical sonic Mach numbers of 3 × 10^{−4}. A relation between the rise velocity of the bubble and the entropy perturbation has been obtained on the basis of physical considerations. We show that the results obtained with this MHD scheme could satisfy that relation even on coarse grids, which suggests that discretization errors arising from the background stratification did not play any significant role in the evolution of the bubble. For comparison, we also ran the same setup at an intermediate entropy perturbation without wellbalancing. The unbalanced states led to the generation of large pressure jumps at the cell interfaces, which launched strong waves that degraded the numerical solution at low resolution. Even when the magnitude of these errors was progressively reduced at higher resolution, they were still significant on the finest grid.
Lastly, we ran a fully convective box on a 3D Cartesian grid (with 2N × N × 2N grid cells) to simulate an SSD. The initial stratification was in MHSE and it resembled the thermodynamic conditions found in oxygen shell burning in a massive star (Andrassy et al. 2022). Convection was driven by a heat source placed at the bottom of the box, and a weak seed magnetic field was planted in the system to initiate the dynamo. By changing the rate of energy injection, we were able to study different velocity regimes. In particular, we simulated three different cases, with ℳ_{son} ~ 10^{−3}, 10^{−2}, and 10^{−1}. We only considered moderate grid resolutions (N = 32, 64, 96). In the kinematic phase, the magnetic field energy was exponentially amplified on the smallest scales of the turbulent flow, with a higher growth rate in finer grids, which is consistent with SSD amplification. When the Lorentz force started to affect the evolution of the fluid, the saturated nonlinear phase began. Because of the use of a lowMach Riemann solver, the amount of magnetic energy amplified (compared to the kinetic energy content of the flow) did not depend on the sonic Mach number of convection and achieved subequipartition values for the resolutions considered in this study (E_{M}/E_{K} ≃ 0.1ℳ0.2).
Overall, the results obtained in these tests demonstrate that the numerical methods implemented in SLH can accurately and efficiently tackle a variety of MHD processes that act in stellar interiors, in regimes that are inaccessible to conventional methods.
Acknowledgements
The work of G.L. and F.K.R. is supported by the German Research Foundation (DFG) through the grant RO 3676/31. The work of C.B. and C.K. is supported by DFG through the grant KL 566/221. P.V.F.E. was supported by the US Department of Energy through the Los Alamos National Laboratory (LANL). LANL is operated by Triad National Security, LLC, for the National Nuclear Security Administration of the US Department of Energy (Contract No. 89233218CNA000001). G.L., R.A., J.H., G.W., and F.K.R. acknowledge support by the Klaus Tschira Foundation. This work has been assigned a document release number LAUR2227864.
Appendix A Magnetized KelvinHelmholtz instability
This appendix explores the effects of the grid resolution and strength of the initial magnetic field on the evolution of the Kelvin–Helmholtz instability shown in Sect. 5.3.
Fig. A.1 Distribution of the sonic Mach number in the simulations of the KelvinHelmholtz instability at t/t_{max} = 1/6 for different values of the initial magnetic field, Β_{x}, computed as , with α = (0.1, 0.2, 0.5, 0.6, 0.7, 0.9, 0.9, 1.0, 1.2). These simulations were performed on a 512 × 256 grid with M_{x} = 10^{−3}. The title in each panel is the corresponding minimum Alfvén Mach number of the flow at t = 0. For a strong enough initial magnetic field, the magnetic stresses prevent the growth of the instability. 
Fig. A.2 Distribution of Β_{x} in the simulations of the Kelvin–Helmholtz instability at t/t_{max} = 1/6 for different grid resolutions, starting from the initial conditions described in Sect. 5.3 with ℳ_{x} = 10^{−3}. On grids with N ≤ 128, numerical discretization errors generate gridscale vorticity, which leads to the growth of secondary instabilities in the regions between the primary rolls. This effect does not appear in better converged simulations. 
Fig. A.3 Time evolution of B_{x} in the simulations of the Kelvin–Helmholtz instability starting from the initial conditions described in Sect. 5.3 with ℳ_{x} = 10^{−3}. The grid is 2048 × 1024. 
Appendix B Hot bubble
Here, we extend the study described in Sect. 5.4. In particular, we show the dependence of the entropy fluctuations, ℳ_{son}, and p_{B} on the magnitude of the initial entropy perturbation (ΔA/A)_{t=0}.
Fig. B.1 Final distribution of the entropy fluctuations, ΔA/A, of the hot bubble for different values of (ΔA/A)_{t=0} and grid resolutions. Each panel is rescaled by the corresponding value of (ΔA/A)_{t=0}. The insets provide the minimum and maximum values of the entropy fluctuations in each plot. 
Fig. B.2 Final distribution of ∣B∣/B_{x,t=0} for different values of (ΔA/A)_{t=0} and grid resolutions in the simulations of the hot bubble. The insets show the maximum ratio in each panel. The amount of numerical resistivity decreases upon grid refinement, which leads to the generation of narrower stripes with stronger magnetic fields. 
Fig. B.3 Final distribution of the sonic Mach number of the hot bubble for different values of (ΔA/A)_{t=0} and different grid resolutions. Each panel is rescaled by the corresponding value of . The insets show the maximum sonic Mach number. An entropy perturbation of (ΔA/A)_{t=0} = 0.1 drives flows that are far from the lowMachnumber regime. In this case, effects of compressibility caused by the high ram pressure of the bubble are large enough to cause a 6–7% deviation from the theoretical scaling discussed in Sect. 5.4. 
Appendix C Smallscale dynamo
In this section we extend the analysis of the SSD test described in Sect. 5.5. In particular, we show ID vertical averages of the velocity and magnetic field distributions, the time evolution of the numerical divergence of the magnetic field, vertical cuts of the sonic Mach number distribution and the time evolution of the total magnetic energy in the kinematic phase.
Fig. C.1 Vertical profiles of V_{y} (left) and Β_{ν} (right) in the simulations of the SSD, averaged over 20 < t/τ_{conv} < 40. Different colors represent different grid resolutions, while different line styles are used for representing different values of b. The vertical velocity is rescaled by b^{2/3} to remove the dependence of the energy injection rate, while Β_{y} is rescaled by the corresponding equipartition value, . 
Fig. C.2 Vertical cuts of the sonic Mach number taken at z = 0 in the simulations of the SSD run with b = 10^{−6} on the 192 × 96 × 192 grid. Each panel is taken at a different moment in time, which is given in the title. As the SSD approaches the saturated regime, the smallscale structures in the velocity field are damped by the Lorentz force, and vertical motions mainly happen in the form of largescale flows. 
Fig. C.3 Time evolution of the maximum (solid line) and mean (dotdashed) relative divergence of the magnetic field in the simulations of the SSD. Since the induction equation is solved using a staggered constrained transport method, the update on the magnetic field keeps the divergence in Eq. (47) within roundoff error. Although these errors accumulate in time, by the end of the simulation magnetic monopoles are still irrelevant to the dynamics of the SSD. 
Fig. C.4 Time evolution of the magnetic energy in the simulations of the SSD up to t/τ_{conv} = 20. Each line style represents a specific value of the boost factor, b, while different colors are used for different numbers of grid cells. 
References
 Andrassy, R., Herwig, F., Woodward, P., & Ritter, C. 2020, MNRAS, 491, 972 [NASA ADS] [CrossRef] [Google Scholar]
 Andrassy, R., Higl, J., Mao, H., et al. 2022, A&A, 659, A193 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Augustson, K. C., Brun, A. S., & Toomre, J. 2016, ApJ, 829, 92 [Google Scholar]
 Aydemir, A. Y., & Barnes, D. C. 1985, J. Comput. Phys., 59, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Balsara, D. S. 2004, ApJS, 151, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Balsara, D. S., & Spicer, D. S. 1999, J. Comput. Phys., 149, 270 [NASA ADS] [CrossRef] [Google Scholar]
 Berberich, J. P., Chandrashekar, P., & Klingenberg, C. 2021, Comput. Fluids, 219, 104858 [Google Scholar]
 Brackbill, J. U., & Barnes, D. C. 1980, J. Comput. Phys., 35, 426 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A. 2011, ApJ, 741, 92 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A. 2014, ApJ, 791, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A., & Subramanian, K. 2005, Phys. Rep., 417, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, B. P., Browning, M. K., Brun, A. S., Miesch, M. S., & Toomre, J. 2010, ApJ, 711, 424 [Google Scholar]
 Browning, M. K. 2008, ApJ, 676, 1262 [Google Scholar]
 Browning, M. K., Miesch, M. S., Brun, A. S., & Toomre, J. 2006, ApJ, 648, L157 [Google Scholar]
 Brun, A. S., & Browning, M. 2017, Living Rev. Sol. Phys., 14, 4 [Google Scholar]
 Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073 [Google Scholar]
 Brun, A. S., Browning, M. K., & Toomre, J. 2005, ApJ, 629, 461 [Google Scholar]
 Canivete Cuissa, J. R., & Teyssier, R. 2022, A&Amp;A 664, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cargo, P., & Gallice, G. 1997, J. Comput. Phys., 136, 446 [NASA ADS] [CrossRef] [Google Scholar]
 Chacón, L. 2008, Phys. Plasmas, 15, 056103 [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability (Oxford University Press) [Google Scholar]
 Charbonneau, P. 2013, Solar and Stellar Dynamos (SpringerVerlag) [CrossRef] [Google Scholar]
 Charlton, L. A., Holmes, J. A., Lynch, V. E., Carreras, B. A., & Hender, T. C. 1990, J. Comput. Phys., 86, 270 [NASA ADS] [CrossRef] [Google Scholar]
 Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174 [NASA ADS] [CrossRef] [Google Scholar]
 Courant, R., Friedrichs, K., & Lewy, H. 1928, Math. Ann., 100, 32 [Google Scholar]
 Dai, W., & Woodward, P. R. 1998, ApJ, 494, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Dedner, A., Kemm, F., Kröner, D., et al. 2002, J. Comput. Phys., 175, 645 [Google Scholar]
 Dumbser, M., Balsara, D. S., Tavelli, M., & Fambri, F. 2019, Int. J. Numer. Methods Fluids, 89, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Edelmann, P. V. F. 2014, Dissertation, Technische Universität München, Germany [Google Scholar]
 Edelmann, P. V. F., & Röpke, F. K. 2016, in JUQUEEN Extreme Scaling Workshop 2016, eds. D. Brömmel, W. Frings, & B. J. N. Wylie (Jülich Supercomputing Centre) [Google Scholar]
 Edelmann, P. V. F., Röpke, F. K., Hirschi, R., Georgy, C., & Jones, S. 2017, A&Amp;A, 604, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Edelmann, P. V. F., Ratnasingam, R. P., Pedersen, M. G., et al. 2019, ApJ, 876, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Edelmann, P. V. F., Horst, L., Berberich, J. P., et al. 2021, A&Amp;A, 652, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Einfeldt, B., Roe, P. L., Munz, C. D., & Sjogreen, B. 1991, J. Comput. Phys., 92, 273 [NASA ADS] [CrossRef] [Google Scholar]
 Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659 [Google Scholar]
 Fambri, F. 2021, Int. J. Numer. Methods Fluids, 93, 3447 [NASA ADS] [CrossRef] [Google Scholar]
 Featherstone, N. A., & Hindman, B. 2016, ApJ, 818, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Featherstone, N. A., Browning, M. K., Brun, A. S., & Toomre, J. 2009, ApJ, 705, 1000 [Google Scholar]
 Felipe, T., Khomenko, E., & Collados, M. 2010, ApJ, 719, 357 [Google Scholar]
 Frank, A., Jones, T. W., Ryu, D., & Gaalaas, J. B. 1996, ApJ, 460, 777 [Google Scholar]
 Fuchs, F. G., Mishra, S., & Risebro, N. H. 2009, J. Comput. Phys., 228, 641 [NASA ADS] [CrossRef] [Google Scholar]
 Gardiner, T. A., & Stone, J. M. 2005, J. Comput. Phys., 205, 509 [NASA ADS] [CrossRef] [Google Scholar]
 Gardiner, T. A., & Stone, J. M. 2008, J. Comput. Phys., 227, 4123 [NASA ADS] [CrossRef] [Google Scholar]
 Gastine, T., & Wicht, J. 2012, Icarus, 219, 428 [Google Scholar]
 Ghizaru, M., Charbonneau, P., & Smolarkiewicz, P. K. 2010, ApJ, 715, L133 [NASA ADS] [CrossRef] [Google Scholar]
 Glasser, A. H., Sovinec, C. R., Nebel, R. A., et al. 1999, Plasma Phys. Controlled Fusion, 41, A747 [Google Scholar]
 Glatzmaier, G. A. 1984, J. Comput. Phys., 55, 461 [NASA ADS] [CrossRef] [Google Scholar]
 Glatzmaier, G. A. 1985, ApJ, 291, 300 [NASA ADS] [CrossRef] [Google Scholar]
 Godunov, S. K., & Bohachevsky, I. 1959, Matematiceskij Sbornik, 47, 271 [Google Scholar]
 Harned, D. S., & Kerner, W. 1985, J. Comput. Phys., 60, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Higl, J., Müller, E., & Weiss, A. 2021, A&Amp;A, 646, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Horst, L., Edelmann, P. V. F., Andrássy, R., et al. 2020, A&Amp;A, 641, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Horst, L., Hirschi, R., Edelmann, P. V. F., Andrassy, R., & Roepke, F. K. 2021, A&Amp;A, 653, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hotta, H. 2017, ApJ, 843, 52 [Google Scholar]
 Hotta, H., Rempel, M., & Yokoyama, T. 2015, ApJ, 803, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Iskakov, A. B., Schekochihin, A. A., Cowley, S. C., McWilliams, J. C., & Proctor, M. R. E. 2007, Phys. Rev. Lett., 98, 208501 [NASA ADS] [CrossRef] [Google Scholar]
 Jardin, S. C. 2012, J. Comput. Phys., 231, 822 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, C. A., Kuzanyan, K. M., & Mitchell, R. H. 2009, J. Fluid Mech., 634, 291 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J. 2011, Astron. Nachr., 332, 43 [CrossRef] [Google Scholar]
 Käpylä, P. J. 2019, A&Amp;A, 631, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J. 2021, A&Amp;A, 651, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J., Korpi, M. J., & Brandenburg, A. 2008, A&Amp;A, 491, 353 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2011, Astron. Nachr., 332, 883 [CrossRef] [Google Scholar]
 Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2012, ApJ, 755, L22 [CrossRef] [Google Scholar]
 Käpylä, P. J., Mantere, M. J., Cole, E., Warnecke, J., & Brandenburg, A. 2013, ApJ, 778, 41 [CrossRef] [Google Scholar]
 Käpylä, P. J., Gent, F. A., Olspert, N., Käpylä, M. J., & Brandenburg, A. 2020, Geophys. Astrophys. Fluid Dyn., 114, 8 [CrossRef] [Google Scholar]
 Karak, B. B., Käpylä, P. J., Käpylä, M. J., et al. 2015, A&Amp;A, 576, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kazantsev, A. P. 1968, Sov. J. Exp. Theor. Phys., 26, 1031 [NASA ADS] [Google Scholar]
 Khomenko, E., & Collados, M. 2006, ApJ, 653, 739 [Google Scholar]
 Kippenhahn, R., Weigert, A., & Weiss, A. 2013, Stellar Structure and Evolution (Springer) [Google Scholar]
 Kolmogorov, A. 1941, Akad. Nauk SSSR Dokl., 30, 301 [Google Scholar]
 Kupka, F., & Muthsam, H. J. 2017, Living Rev. Comput. Astrophys., 3, 1 [Google Scholar]
 Lerbinger, K., & Luciani, J. F. 1991, J. Comput. Phys., 97, 444 [NASA ADS] [CrossRef] [Google Scholar]
 LeVeque, R. J. 2002, Finite Volume Methods for Hyperbolic Problems, Cambridge Texts in Applied Mathematics (Cambridge University Press) [Google Scholar]
 Londrillo, P., & del Zanna, L. 2004, J. Comput. Phys., 195, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Lütjens, H., & Luciani, J.F. 2010, J. Comput. Phys., 229, 8130 [CrossRef] [Google Scholar]
 Masada, Y., Yamada, K., & Kageyama, A. 2013, ApJ, 778, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Matthaeus, W. H., & Brown, M. R. 1988, Phys. Fluids, 31, 3634 [NASA ADS] [CrossRef] [Google Scholar]
 Meneguzzi, M., Frisch, U., & Pouquet, A. 1981, Phys. Rev. Lett., 47, 1060 [NASA ADS] [CrossRef] [Google Scholar]
 Mestel, L. 1999, Stellar Magnetism (Oxford University Press) [Google Scholar]
 Miczek, F. 2013, Dissertation, Technische Universität München, Germany [Google Scholar]
 Miczek, F., Röpke, F. K., & Edelmann, P. V. F. 2015, A&Amp;A, 576, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mignone, A., & Del Zanna, L. 2021, J. Comput. Phys., 424, 109748 [NASA ADS] [CrossRef] [Google Scholar]
 Minoshima, T., & Miyoshi, T. 2021, J. Comput. Phys., 446, 110639 [CrossRef] [Google Scholar]
 Minoshima, T., Kitamura, K., & Miyoshi, T. 2020, ApJS, 248, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Miyoshi, T., & Kusano, K. 2005, J. Comput. Phys., 208, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Müller, B. 2020, Living Rev. Comput. Astrophys., 6, 3 [CrossRef] [Google Scholar]
 Pietarila Graham, J., Cameron, R., & Schüssler, M. 2010, ApJ, 714, 1606 [Google Scholar]
 Powell, K. G. 1997, in Upwind and HighResolution Schemes, eds. M.Y. Hussaini, B. van Leer, J. Van Rosendale (Springer), 570 [CrossRef] [Google Scholar]
 Powell, K. G., Roe, P. L., Linde, T. J., Gombosi, T. I., & De Zeeuw, D. L. 1999, J. Comput. Phys., 154, 284 [NASA ADS] [CrossRef] [Google Scholar]
 Rempel, M. 2005, ApJ, 622, 1320 [NASA ADS] [CrossRef] [Google Scholar]
 Rempel, M. 2014, ApJ, 789, 132 [Google Scholar]
 Rempel, M. 2018, ApJ, 859, 161 [Google Scholar]
 Riva, F., & Steiner, O. 2022, A&Amp;A, 660, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rogers, T. M., Lin, D. N. C., McElwaine, J. N., & Lau, H. H. B. 2013, ApJ, 772, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Schekochihin, A. A., Cowley, S. C., Taylor, S. F., Maron, J. L., & McWilliams, J. C. 2004, ApJ, 612, 276 [NASA ADS] [CrossRef] [Google Scholar]
 Schekochihin, A. A., Iskakov, A. B., Cowley, S. C., et al. 2007, N. J. Phys., 9, 300 [NASA ADS] [CrossRef] [Google Scholar]
 Schnack, D., Barnes, D., Mikic, Z., Harned, D. S., & Caramana, E. 1987, J. Comput. Phys., 70, 330 [NASA ADS] [CrossRef] [Google Scholar]
 Seta, A., & Federrath, C. 2020, MNRAS, 499, 2076 [NASA ADS] [CrossRef] [Google Scholar]
 Shu, C.W., & Osher, S. 1988, J. Comput. Phys., 77, 439 [NASA ADS] [CrossRef] [Google Scholar]
 Smolarkiewicz, P. K., & Charbonneau, P. 2013, J. Comput. Phys., 236, 608 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Strang, G. 1968, SIAM J. Numer. Anal., 5, 506 [CrossRef] [MathSciNet] [Google Scholar]
 Sutherland, B. R. 2010, Internal Gravity Waves (Cambridge University Press) [Google Scholar]
 Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501 [NASA ADS] [CrossRef] [Google Scholar]
 Toro, E. F. 2009, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction (Springer) [CrossRef] [Google Scholar]
 Tóth, G. 2000, J. Comput. Phys., 161, 605 [Google Scholar]
 Viallet, M., Baraffe, I., & Walder, R. 2011, A&Amp;A, 531, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Viviani, M., Käpylä, M. J., Warnecke, J., Käpylä, P. J., & Rheinhardt, M. 2019, ApJ, 886, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&Amp;A, 429, 335 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Warnecke, J., Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2016, A&Amp;A, 596, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yadav, R. K., Christensen, U. R., Wolk, S. J., & Poppenhaeger, K. 2016, ApJ, 833, L28 [Google Scholar]
A similar approach in which the states are split into a background component and deviations is described in Vögler et al. (2005), Khomenko & Collados (2006), Felipe et al. (2010), and Hotta et al. (2015).
If the gravitational potential is time independent, solving the energy equation for ρEϕ instead of ρE allows the ρg · V source term to be removed. This leads to more accurate results and better entropy and energyconservation properties in simulations of gas dynamics with gravity (Müller 2020; Edelmann et al. 2021).
Several spatial reconstruction routines are implemented in SLH, from simple constant extrapolation to the Piecewise Parabolic Method of Colella & Woodward (1984). These reconstruction schemes can be applied to both conservative and primitive variables.
For more details on the implementation of implicit time stepping in SLH, see Miczek (2013) and Miczek et al. (2015).
We took these values as representative of the typical conditions found in stellar convection zones close to equipartition regimes (Augustson et al. 2016).
More details on how to compute the pressure profile can be found in Edelmann et al. (2021).
As observed by Seta & Federrath (2020), the evolution of the dynamo in the nonlinear regime does not depend on the form of the seed field, as long as its magnitude is weak enough to not affect the development of convection.
Details on how to compute the density perturbation can be found in Andrassy et al. (2022).
All Tables
Time averages of E_{M}/b^{2/3} [10^{45} × erg] over the time interval 20 < t/τ_{conv} < 70 for the different resolutions, N, and boost factors, b, considered in this study.
All Figures
Fig. 1 Wave structure of the MHD system. 

In the text 
Fig. 2 Wave structure of HLLDtype solvers: only the fast magnetosonic (S_{L}, S_{R}), Alfvén , and entropy (S_{M}) waves are considered in the computation of the states across the Riemann fan. The slow magnetosonic waves are discarded. 

In the text 
Fig. 3 Global error as a function of the grid resolution for the seven MHD waves after one crossing time, t_{c}. The dashed black line represents the secondorder scaling. 

In the text 
Fig. 4 Global error as a function of the grid resolution for the leftgoing fast magnetosonic wave. Different colors are associated with different values of c_{f,x}/c_{a,x}. The dashed black line is the secondorder scaling. 

In the text 
Fig. 5 Magnetic energy distribution of the Balsara vortex after one advective crossing time, t_{adv}, 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 yaxis (in descending order), while the initial maximum rotational velocity, Ṽ, varies along the xaxis. The inset in each subplot shows the ratio of the final to the initial magnetic energy. The vortex run with Ṽ = 10^{−1} and β_{K} = 10^{2} (bottomright 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. 

In the text 
Fig. 6 Convergence of the L_{1} 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 secondorder scaling. 

In the text 
Fig. 7 Distribution of the rotational kinetic energy (normalized by the maximum initial value) of the Balsara vortex after one advective time, t_{adv}, 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: /(E_{R,t=0})_{tot} − 1. 

In the text 
Fig. 8 Ratio of the wall clock times obtained with SSPRK2 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. 

In the text 
Fig. 9 Initial setups of V_{x} and V_{y} (here rescaled by ℳ_{x}) used for simulating the growth of the Kelvin–Helmholtz instability. 

In the text 
Fig. 10 Time evolution of the ydirection kinetic energy rescaled by its initial value in the magnetized Kelvin–Helmholtz instability test problem. Each panel corresponds to a different initial Mach number, ℳ_{x}. Different colors are used for different grid resolutions (the 64 × 32 and 128 × 64 grids cells have been left out for clarity). Results obtained with the HLLD solver are represented by dotdashed lines, while solid lines are used for LHLLD. The solid black line in each panel is the reference solution. As explained in the text, the nonlinear solver does not converge when using HLLD at ℳ_{x} = 10^{−4} for N > 64. 

In the text 
Fig. 11 Same as Fig. 10 but showing the total magnetic energy divided by its initial value. 

In the text 
Fig. 12 Distribution of the sonic Mach number in the KelvinHelmholtz instability test at t/t_{max} =1/6 obtained with the HLLD (top panels) and the LHLLD (bottom panels) solvers on a 128 × 64 grid for different values of ℳ_{x}. All panels are rescaled by the corresponding value of M_{x}. 

In the text 
Fig. 13 Convergence with resolution N of the L_{1} error associated with E_{K,y} rescaled by in the simulations of the Kelvin–Helmholtz instability. Different colors are used for different initial sonic Mach numbers ℳ_{x} using the LHLLD (solid lines) and HLLD (dotdashed lines) solvers. The dashed black line is the secondorder scaling. 

In the text 
Fig. 14 Maximum sonic Mach number, minimum and maximum entropy fluctuations and Alfvén speed of the hot bubble as a function of the initial entropy perturbation obtained on a 64 × 96 grid. The black lines represent the physical scalings. 

In the text 
Fig. 15 Final distribution of the entropy fluctuations of the hot bubble for (ΔA/A)_{t=0} = 10^{−3} at different grid resolutions. The entropy fluctuations are rescaled by (ΔA/A)_{t=0}. The top row shows the results obtained with deviation wellbalancing, whereas no wellbalancing method was used in the simulations shown in the bottom row. The insets show the minimum and maximum values of the entropy fluctuation in each panel. 

In the text 
Fig. 16 Time evolution of the magnetic energy (sky blue) and kinetic energy (vermillion) in the simulations of the SSD for different grid resolutions (dotted: 64 × 32 × 64, dotdashed: 128 × 64 × 128, solid: 192 × 96 × 192). Each panel shows the results obtained with a specific value of b (from left to right: b = 10^{−6}, 10^{−3}, 1). The time is expressed in units of the convective turnover, while the magnetic and kinetic energy curves are rescaled by b^{2/3} to take into account the different energy contents of the flows. 

In the text 
Fig. 17 Horizontal slices in the y = 1.5 plane taken in the kinematic (upper plots) and saturated (lower plots) regimes of the SSD with b = 10^{−6} on the 192 × 96 × 192 grid. The panels on the left show the vertical sonic Mach number ℳ_{son,y} = V_{y/a} multiplied by 10^{3}, while the plots on the right show the vertical magnetic field rescaled by the root mean square value across the plane. 

In the text 
Fig. 18 Root mean square of the sonic Mach number as a function of the driving luminosity b. The data points are averages computed in the time interval 20 < t/τ_{com} < 40. The error bars represent three standard deviations over the time series, while the dotted black line is the b^{1/3} scaling. 

In the text 
Fig. 19 Kinetic (solid lines) and magnetic (dashed lines) energy spectra as a function of the wave number k obtained in an horizontal slice of the convective box at y = 1.5. All spectra are averaged over the time interval 20 < t/τ_{conv} < 40 and are rescaled by b^{2/3}. The three panels show the results for different numbers of grid cells (from top to bottom: 64 × 32 × 64, 128 × 64 × 128, 192 × 96 × 192). The dotted black lines are the Kolmogorov (k^{−5/3}) and Kazantsev (k^{3/2}) scalings. 

In the text 
Fig. A.1 Distribution of the sonic Mach number in the simulations of the KelvinHelmholtz instability at t/t_{max} = 1/6 for different values of the initial magnetic field, Β_{x}, computed as , with α = (0.1, 0.2, 0.5, 0.6, 0.7, 0.9, 0.9, 1.0, 1.2). These simulations were performed on a 512 × 256 grid with M_{x} = 10^{−3}. The title in each panel is the corresponding minimum Alfvén Mach number of the flow at t = 0. For a strong enough initial magnetic field, the magnetic stresses prevent the growth of the instability. 

In the text 
Fig. A.2 Distribution of Β_{x} in the simulations of the Kelvin–Helmholtz instability at t/t_{max} = 1/6 for different grid resolutions, starting from the initial conditions described in Sect. 5.3 with ℳ_{x} = 10^{−3}. On grids with N ≤ 128, numerical discretization errors generate gridscale vorticity, which leads to the growth of secondary instabilities in the regions between the primary rolls. This effect does not appear in better converged simulations. 

In the text 
Fig. A.3 Time evolution of B_{x} in the simulations of the Kelvin–Helmholtz instability starting from the initial conditions described in Sect. 5.3 with ℳ_{x} = 10^{−3}. The grid is 2048 × 1024. 

In the text 
Fig. B.1 Final distribution of the entropy fluctuations, ΔA/A, of the hot bubble for different values of (ΔA/A)_{t=0} and grid resolutions. Each panel is rescaled by the corresponding value of (ΔA/A)_{t=0}. The insets provide the minimum and maximum values of the entropy fluctuations in each plot. 

In the text 
Fig. B.2 Final distribution of ∣B∣/B_{x,t=0} for different values of (ΔA/A)_{t=0} and grid resolutions in the simulations of the hot bubble. The insets show the maximum ratio in each panel. The amount of numerical resistivity decreases upon grid refinement, which leads to the generation of narrower stripes with stronger magnetic fields. 

In the text 
Fig. B.3 Final distribution of the sonic Mach number of the hot bubble for different values of (ΔA/A)_{t=0} and different grid resolutions. Each panel is rescaled by the corresponding value of . The insets show the maximum sonic Mach number. An entropy perturbation of (ΔA/A)_{t=0} = 0.1 drives flows that are far from the lowMachnumber regime. In this case, effects of compressibility caused by the high ram pressure of the bubble are large enough to cause a 6–7% deviation from the theoretical scaling discussed in Sect. 5.4. 

In the text 
Fig. C.1 Vertical profiles of V_{y} (left) and Β_{ν} (right) in the simulations of the SSD, averaged over 20 < t/τ_{conv} < 40. Different colors represent different grid resolutions, while different line styles are used for representing different values of b. The vertical velocity is rescaled by b^{2/3} to remove the dependence of the energy injection rate, while Β_{y} is rescaled by the corresponding equipartition value, . 

In the text 
Fig. C.2 Vertical cuts of the sonic Mach number taken at z = 0 in the simulations of the SSD run with b = 10^{−6} on the 192 × 96 × 192 grid. Each panel is taken at a different moment in time, which is given in the title. As the SSD approaches the saturated regime, the smallscale structures in the velocity field are damped by the Lorentz force, and vertical motions mainly happen in the form of largescale flows. 

In the text 
Fig. C.3 Time evolution of the maximum (solid line) and mean (dotdashed) relative divergence of the magnetic field in the simulations of the SSD. Since the induction equation is solved using a staggered constrained transport method, the update on the magnetic field keeps the divergence in Eq. (47) within roundoff error. Although these errors accumulate in time, by the end of the simulation magnetic monopoles are still irrelevant to the dynamics of the SSD. 

In the text 
Fig. C.4 Time evolution of the magnetic energy in the simulations of the SSD up to t/τ_{conv} = 20. Each line style represents a specific value of the boost factor, b, while different colors are used for different numbers of grid cells. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.