Open Access
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/0004-6361/202244665
Published online 19 December 2022

© G. Leidi et al. 2022

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

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

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

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

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

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

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

2 Equations of compressible ideal MHD with gravity

The MHD scheme implemented in SLH is designed to solve the set of compressible ideal MHD equations with a (time-independent) gravitational source term3: (1) (2) (3) (4)

where ρ denotes the density, V = (Vx, Vy, Vz) the velocity field, B = (Bx, By, Bz) the magnetic field4, I the unit tensor, g = (gx, gy, gz) the gravitational acceleration, p the gas pressure and pB = ∣B2/2 the magnetic pressure. The total energy density ρEϕ is defined as (5)

where eint and eϕ are the specific internal and gravitational energies5.

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

2.1 Eigenstructure of the MHD system and definition of Mach numbers

The homogeneous MHD system (left-hand side of Eqs. (1)(4)) reduced to one spatial dimension has seven nonzero eigenvalues6, (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. cs,x, ca,x, and cf,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)

thumbnail Fig. 1

Wave structure of the MHD system.

2.2 Low-Mach limit of the MHD system

Magnetic fields amplified by dynamo mechanisms in deep convective layers of stars are likely to approach equipartition with respect to the kinetic energy content of the flow (Brandenburg & Subramanian 2005; Featherstone et al. 2009; Augustson et al. 2016; Hotta 2017; Käpylä 2019). To model such processes, the MHD system in Eqs. (1)(4) must be solved in regimes of Mach numbers ℳfast,x ≲ ℳson ≤ 10−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 equations7, (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 low-Mach limit. If the plasma-β is high, the large fast magnetosonic speed cf,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 finite-volume method (LeVeque 2002; Toro 2009), which is briefly summarized in the next section.

3.1 Finite-volume discretization

In a first step, the physical system is mapped on a 3D Cartesian grid8 divided into Nx × Ny × Nz cells, whose spatial extent is given by [xL, xR] × [yL,yR] x [zL, zR]. 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 volume-averaged vector of conserved quantities (23)

The same procedure applies to Ŝi,j,k, while the surface-averaged fluxes are defined as (24)

where Ai+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 right-hand 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 volume-averaged 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 reconstructed9 (through 1D sweeping) to the center of each cell boundary, starting from the cell-centered 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 surface-averaged flux is then approximated (to second-order 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 fluxes10, 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 low-Mach regimes, discontinuities in the flow are only transported by the linearly degenerate entropy and Alfvén waves. These modes propagate with small speeds (see Sect. 2.2), and by the time they cross one cell in the computational grid they are strongly dissipated by the action of the numerical term in Eq. (27). This effect can also be explained by noticing that the pressure-diffusion coefficient scales as O(1/ℳfast,x), so that it overwhelms the physical flux proportional to O(1) at low sonic Mach numbers. To remove this excessive dissipation, we use the low-dissipation HLLD solver (LHLLD) of Minoshima & Miyoshi (2021). This is a variation of the original five-wave HLLD solver (see Fig. 2) of Miyoshi & Kusano (2005). LHLLD introduces a Mach-dependent parameter ϕ ∝ ℳfast,x in the intermediate state of the total pressure pT = p + pB: (28)

In this context, SL and SR are conservative estimates of the speeds λ1,7. In SLH they are evaluated as (29)

The low-Mach fix ϕ is computed according to the following formulas: (30)

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

thumbnail Fig. 2

Wave structure of HLLD-type solvers: only the fast magnetosonic (SL, SR), Alfvén , and entropy (SM) waves are considered in the computation of the states across the Riemann fan. The slow magnetosonic waves are discarded.

3.3 Well-balancing method

As already noted in Sects. 1 and 3.1, hyperbolic fluxes and gravitational source terms are discretized with different methods. As a consequence, Godunov-type schemes do not automatically preserve magnetohydrostatic solutions on a discrete grid exactly. Therefore, whenever a stratification needs to be enforced to be in MHSE on the computational grid, we use the deviation well-balancing method (Berberich et al. 2021; Edelmann et al. 2021). The main ingredient of this method is an a priori known target state Ũ that is a magnetohydrostatic solution to Eq. (21), (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 well-balanced method, Eq. (32) is discretized according to the finite-volume method described in Sect. 3.1, which leads to the semi-discrete form (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 ΔUi,j,k, rather than the states Ui,j,k, are reconstructed to the boundary of the cell11. This guarantees that magnetohydrostatic solutions are preserved on the discrete grid, since in that case ΔUi,j,k = 0, which leads to (39)

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

3.4 Constrained transport method

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

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

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

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

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

The key point of staggered constrained transport methods is to compute the surface integral of Eq. (19) over cell boundaries using Stokes’s theorem, which leads to the finite-area equation12 (40)

Here, is the surface-averaged magnetic field component normal to the cell boundary (41)

while the line-averaged 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 surface-averaged quantity with its value at the center of the cell boundary, (43)

which is accurate to second order. In contrast to the standard finite-volume approach, here the magnetic field component normal to the interface is stored at cell boundaries, while the line-averaged electromotive force is evaluated at cell edges. Thus, the operation is performed on a staggered grid. Since the parallel magnetic field still needs to be reconstructed to compute the flux function, its value at cell-center locations is estimated as a simple arithmetic average between the neighboring cell interfaces: (44)

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

is approximated to second-order accuracy by (45)

where can be computed from the solution to the Riemann problem in Eq. (26). The calculation for the x- and y-component is again analogous. The upwind diffusion term enters in the derivatives of the electromotive force in Eq. (45), which are obtained according to the sign si+1/2,j,,k of the entropy (contact) waves at the cell interfaces: (46)

Here represents the z-component of the cell-centered electromotive force. The discretization of the line-averaged electromotive force leads to a semi-discrete form of Eq. (40) that can be integrated numerically in time. Any time-stepper that solves the resulting system of ODEs can keep the cell-volume average of ∇ · B, (47)

within rounding errors.

4 Time integration algorithm

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

As outlined in Sect. 1, we split the induction equation (see Eq. (4)) from the continuity, momentum and energy equations (see Eqs. (1)(3)), based on the approach described by Fuchs et al. (2009). This allows different spatial and temporal discretizations to be used depending on the problem at hand. In regimes of low Mach numbers and high β values, the stiffness is mostly generated by the pressure flux in the momentum equation, while the nondimensional form of the induction equation does not depend on the Mach number of the flow (see Sect. 2.2). This suggests that implicit time discretization only needs to be applied to the subset of continuity, momentum and energy equations, whereas the induction equation can be solved with explicit time-steppers. These two updates can be combined to second-order accuracy with Strang splitting (Strang 1968): (48)

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

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

Numerical experiments performed with the proposed implicit-explicit Strang splitting (IESS) approach suggest that the maximum time step allowed for stability is approximately determined by (49)

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

A single step of the described time-marching scheme can be summarized in the following way. First, Δt is obtained from Eq. (49) given Bn, ρn, ρVn, and . If gravity is not present, does not appear in Eq. (5).

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

Third, this intermediate solution, Bn+1/2, is used to solve the continuity, momentum, and energy equations over the full time step, Δt, with an implicit time-stepper. If gravity is present, then the well-balancing method described in Sect. 3.3 can be used. Any other source term is also considered in this step. This allows the solution for density, momentum, and energy to be obtained at the next step, ρn+1, ρVn+1, and .

Fourth, Bn+1/2, ρn+1, ρVn+1, and are used to solve the induction equation over ΔT/2. This yields the magnetic field at the final step Bn+1.

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

5 Numerical tests

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

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

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

Then we considered two setups in which gravity is present (iv, v). To check the entropy-conservation properties of the scheme based on the deviation well-balancing method, we modeled the rise of a parcel of fluid with higher entropy content than the (isentropic) background stratification, that is, a “hot bubble” (iv). By changing the magnitude of the entropy perturbation, we simulated different rise velocities of the bubble, down to Mach numbers of ℳson ~ 10−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 well-balancing.

Finally, we simulated a fully 3D small-scale dynamo (SSD) amplification in a star-like environment at moderate grid resolutions (ν). By changing the rate at which energy is injected in the system, we simulated progressively slower flows down to ℳson ~ 10−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 time-marching scheme described in Sect. 4, ESDIRK23 was chosen to treat the implicit part of the algorithm. This guarantees second-order accuracy in time. The time step in Eq. (49) was reduced by 20% to get a more conservative stability criterion. Finally, unlimited linear reconstruction, which is second-order accurate in space, was applied to primitive variables. Overall, the proposed scheme is (globally) second-order accurate.

5.1 Linear analysis

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

The homogeneous MHD equations were solved on a periodic ID Cartesian grid divided into N cells, with the spatial domain ranging from 0 to L = 1. For the ith wave, the solution at t = 0 was obtained by perturbing a uniform medium w0 = (1,0,0,0,103/γ, 1, , 1/2) with δw = ARi(w0) sin(2πx/L), where w is the vector of primitive variables (ρ, Vx, Vy, Vz, p, Βx, By, Βz), A = 10−4 is the amplitude of the perturbation, x is the spatial coordinate and Ri is the ith column of the right-eigenvector matrix (see Appendix A.3 in Stone et al. 2008). The chosen values for w0 are such that cf,x ≃ 40.85, ca,x = 1 and cs,x ≃ 0.99. For the entropy mode, we set Vx = 1, so the sonic Mach number of the wave is ℳson ≃ 0.032. The simulations were run for one crossing time defined as tc = L/λi, where λi is the wave speed. The L1 error was then computed for each primitive variable wk as (50)

and the global error associated with the ith wave was then computed as (51)

with (52)

In Eq. (50), j is the spatial index. The tolerance of the Raphson-Newton 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 second-order methods to treat both the spatial and the temporal parts, the scheme converges with second-order accuracy for all waves except the fast magnetosonic waves, which are characterized by much larger errors. This is expected since IESS allows the MHD equations to be integrated over much longer time steps than the CFL constraint (see Sect. 4). As a consequence, the propagation of fast magnetosonic waves is not properly resolved in time and discretization errors strongly deteriorate the numerical solution. This effect is further quantified in Fig. 4, where we show the global error associated with the left-going fast magnetosonic wave as a function of the grid resolution for different values of the gas pressure p, such that cf,x/ca,x = 2, 3, 6, 12, 27, 59, 129. Overall, the errors tend to decrease for smaller values of cf,x/ca,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 cf,x/ca,x ≤ 27 converge with second-order accuracy on the grids considered in this study.

thumbnail Fig. 3

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

thumbnail Fig. 4

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

5.2 Balsara vortex

In the previous section, we demonstrated that the current MHD scheme is capable of simulating linear waves with the expected (second-order) scaling with respect to resolution on 1D grids. In order to check the scaling in 2D and to test the low-Mach capabilities of the scheme, we considered the MHD vortex first described by Balsara (2004). This is an exact stationary solution of the ideal 2D homogeneous MHD equations, in which the distribution of the centrifugal acceleration, magnetic tension, gas and magnetic pressure gradients is such that the vortex is stable. The spatial domain is (x, y) ∈ [−5, 5] × [−5, 5], and we used 64 × 64 grid cells with periodic boundaries in both directions. The initial conditions are given by (53)

with r2 = x2 + y2. 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 ∣Vadv∣ = . 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 tadv. Numerical dissipation converts a fraction of kinetic and magnetic energy into internal energy, but the shape of the vortex is well preserved in all runs. The dissipation rate is virtually independent of ℳson. In contrast, dissipation of magnetic energy depends on the value of βK. As already pointed out in Sect. 3.2, the pressure-diffusion coefficient in LHLLD has a residual scaling O(1/ℳAlf). A larger value of βK corresponds to lower ℳAlf, which then increases the magnitude of the numerical dissipation. The velocity field is progressively more diffused out and becomes less efficient in sustaining the magnetic field through induction against numerical resistivity.

To check the convergence of the scheme in 2D, we ran a vortex with V = 10−3 (corresponding to max(ℳson)t=0 = 1.55 × 10−3) and βK = 1 at different resolutions14. At the end of the simulation, the L1 error was computed for each primitive variable wk as (55)

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

To compare the amount of numerical dissipation introduced by a standard and a low-Mach MHD flux function, we reran this last set of simulations with HLLD. In Fig. 7, we show the final rotational kinetic energy distribution obtained with the two methods: (56)

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

As explained in Sect. 4, one advantage of IESS is that the MHD equations can be integrated on time steps longer than that allowed by the CFL condition without sacrificing stability. However, a single step of the proposed scheme is much more expensive than a single step of a more standard time-explicit marching scheme, as a large nonlinear system has to be solved iteratively with a Raphson–Newton method. Because of these competing effects, we expect the IESS scheme to be more efficient than an explicit time-stepper below a certain Mach number. To determine this threshold, we ran sets of simulations with the parameters (57)

using both IESS and the explicit SSP-RK215 on 40 × 40 grid cells. Every other sub-step of the Godunov method (like the spatial reconstruction, the LHLLD flux function and constrained transport) remained unchanged, so the only difference was in the time discretization. At the end of each simulation, the ratio of the wall-clock times WCTSSP–RK2/WCTIESS was taken as a measure of the relative efficiency between the marching schemes16. The results are shown in Fig. 8. As expected, the speed-up of IESS increases as the Mach number of the vortex is decreased. The simulations with βK = 10 are slower than the other cases, as the larger Alfvén speed considerably reduces the time step estimate in Eq. (49), while no significant difference is seen between βK = 0.1 and βK = 1.0. IESS overtakes SSP-RK2 at max(ℳson)t=0 ≃ 4 × 10−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 SSP-RK2. This justifies the implementation efforts of apartially implicit time discretization algorithm for modeling slow flows.

thumbnail Fig. 5

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

5.3 Magnetized Kelvin–Helmholtz instability

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

with (59)

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

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

The instability is started by adding a perturbation to the y-velocity component in the initial state, Vy = 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 (tmax = 4.8/ℳx). The chosen initial conditions are such that the interface across the shear flow is smooth and resolved, which leads to convergent results at least in the early stages of the evolution of the flow. As in the previous test, we compared the results obtained with both the HLLD and LHLLD solvers.

Figures 10 and 11 show the time evolution of the y-direction kinetic energy and the total magnetic energy EM = Σij |Bij|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/tmax ≃ 0.25), EK,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/tmax ≃ 0.45, which violently decouples the primary rolls and causes a secondary peak in EK,y at t/tmax ≃ 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 EK,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, EK,y converges until the major numerical reconnection event affects the velocity field. As shown in Fig. 10, the HLLD solver requires more resolution to reach convergence as the setup is run at progressively lower sonic Mach numbers. Eventually, the Mach-dependent pressure-diffusion coefficient in Eq. (27) completely dominates the evolution of the flow and deteriorates the numerical solution. For this reason, at ℳx = 10−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/tmax = 1/6. While in moderately subsonic regimes the large-scale structures in the flow are qualitatively similar, for lower Mach numbers HLLD introduces progressively more dissipation and the instability is eventually halted. When LHLLD is used instead, the morphology of the flow seems to be independent of the Mach number.

Finally, we performed a quantitative convergence study by computing the L1 error associated with EK,y at t/tmax = 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 L1 error was computed against a reference solution, which was taken from the highest grid resolution runs considered in this test (N = 1024) using the LHLLD solver. All simulations (including the reference solutions) were down-sampled to a 64 × 32 grid, so that the errors could directly be computed for different resolutions. This analysis was repeated for different values of ℳx using both HLLD and LHLLD. The results are shown in Fig. 13. The errors are rescaled by . so that curves corresponding to different sonic Mach numbers lie on the same scale. Overall, the convergence is second-order with N for all simulations. LHLLD provides almost identical (rescaled) errors at given resolution in different regimes of Mach numbers. This is expected because the numerical dissipation introduced by this solver does not depend on ℳson, thanks to the low-Mach fix in Eq. (30). Instead, the errors computed for the HLLD runs show a clear dependence on the sonic Mach number, and the errors get larger for slower flows. In particular, at ℳx = 0.1, HLLD needs approximately 1.2 times the resolution of LHLLD to achieve the same accuracy, which justifies the use of HLLD in this regime of Mach numbers. Instead, when ℳx = 10−2 and ℳx = 10−3, HLLD needs respectively twice or four times the resolution to be as accurate as the low-dissipation flux, which increases the amount of computing time by 8 or 64. Thus, the use of a low-Mach approximate Riemann solver becomes indispensable for providing accurate results in regimes of low sonic Mach numbers with moderate grid resolutions, which would be unfeasible with more standard solvers.

thumbnail Fig. 6

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

thumbnail Fig. 7

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

thumbnail Fig. 8

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

5.4 Hot bubble

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

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

where g0 = −1.09904373 × 105 cm s−2, ky = 2π/Ly, y is the vertical spatial coordinate, and Ly is the vertical extent of the grid. The value of g0 is set such that the ratio of the maximum to the minimum gas pressure18 p(y) is 100, which corresponds to 4.6 pressure scale heights. The entropy profile inside the bubble is given by (62)

where A0 is background entropy, r0 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. B0 = 47.3 was chosen such that the final Alfvén Mach number at the position of largest entropy is in the range ℳAlf ≃ 2–3 depending on the grid resolution. Thus, magnetic fields are dynamically important but not strong enough to suppress buoyancy.

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

According to Eq. (67), the sonic Mach number of the bubble is expected to scale as the square root of the initial entropy perturbation. Any deviation from this relation, which has been obtained on the basis of physical arguments, can be due to difficulties in modeling slow flows in a stratified setup and the build-up of significant numerical errors. In Fig. 14, we show this scaling for the coarsest grid resolution. All data points overlap with the theoretical curve, and the minimum Mach number ℳson achieved in this parameter study is 3.32 × 10−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 well-balancing, where the vertical resolution Ny 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 well-balancing techniques are fundamental to correctly simulate the evolution of small entropy perturbations in steep isentropic stratifications and to reduce the effects of numerical errors when using moderately coarse grids.

thumbnail Fig. 9

Initial setups of Vx and Vy (here rescaled by ℳx) used for simulating the growth of the Kelvin–Helmholtz instability.

thumbnail Fig. 10

Time evolution of the y-direction 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 dot-dashed 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.

thumbnail Fig. 11

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

thumbnail Fig. 12

Distribution of the sonic Mach number in the Kelvin-Helmholtz instability test at t/tmax =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 Mx.

thumbnail Fig. 13

Convergence with resolution N of the L1 error associated with EK,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 (dot-dashed lines) solvers. The dashed black line is the second-order scaling.

thumbnail 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 Small-scale 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 large-scale 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 Pm = ν/η (Schekochihin et al. 2004, 2007; Pietarila Graham et al. 2010; Brandenburg 2011, 2014), here we do not perform a parameter study for Pm, 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 low-Machnumber 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 Nx × Ny × Nz = 2N × N × 2N grid cells and the spatial domain (normalized by a characteristic length L = 4 × 108 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 g0 = 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 SSP-RK2, 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 time-steppers (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 system19, (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)

where (77)

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 convection20. 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 small-scale structures with mixed polarity, while the velocity field is distributed on slightly larger scales, which suggests that Pm ≳ 1. The growth rate γM = d(lnEM)/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 sub-equipartition values EM/EK ≃ 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 EM/EK 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 well-balancing 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 core-convective 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 b2/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, k3/2 (Kazantsev 1968). This can be explained by the fact that, in this setup, turbulence is not isotropic, and large-scale 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 Mach-independent, 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 magneto-convection and SSDs in regimes of low sonic Mach numbers even with moderate grid resolutions.

thumbnail 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 well-balancing, whereas no well-balancing 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.

thumbnail 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, dot-dashed: 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 b2/3 to take into account the different energy contents of the flows.

thumbnail 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 = Vy/a multiplied by 103, while the plots on the right show the vertical magnetic field rescaled by the root mean square value across the plane.

Table 1

Time averages of EM/b2/3 [1045 × erg] over the time interval 20 < t/τconv < 70 for the different resolutions, N, and boost factors, b, considered in this study.

thumbnail 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 b1/3 scaling.

thumbnail 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 b2/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 (k3/2) scalings.

6 Summary and conclusions

In this work, we have presented a new finite-volume 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 magneto-convection and dynamo processes in deep layers of stars. This method relies on a low-dissipation MHD Riemann solver (LHLLD; Minoshima & Miyoshi 2021) to avoid the excessive numerical dissipation typical of high-resolution, shock-capturing solvers as ℳson → 0.

The strict CFL condition on the time step is overcome by using an implicit-explicit time discretization algorithm, for which the induction equation is integrated using an explicit time-stepper, 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 speed-up in regimes of low sonic Mach numbers.

Whenever required, a magnetohydrostatic solution can be enforced on the discrete grid with the deviation well-balancing method (Berberich et al. 2021; Edelmann et al. 2021). This technique leads to better entropy-conservation 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 CT-contact 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 second-order 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 second-order 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 SSP-RK2 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 time-steppers.

In the third experiment we considered the growth of a Kelvin-Helmholtz 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 second-order convergence of the y-direction 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 well-balancing. 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 low-Mach 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 sub-equipartition values for the resolutions considered in this study (EM/EK ≃ 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/3-1. The work of C.B. and C.K. is supported by DFG through the grant KL 566/22-1. 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 LA-UR-22-27864.

Appendix A Magnetized Kelvin-Helmholtz 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.

thumbnail Fig. A.1

Distribution of the sonic Mach number in the simulations of the Kelvin-Helmholtz instability at t/tmax = 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 Mx = 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.

thumbnail Fig. A.2

Distribution of Βx in the simulations of the Kelvin–Helmholtz instability at t/tmax = 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 grid-scale 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.

thumbnail Fig. A.3

Time evolution of Bx 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 pB on the magnitude of the initial entropy perturbation (ΔA/A)t=0.

thumbnail 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.

thumbnail Fig. B.2

Final distribution of ∣B/Bx,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.

thumbnail 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 low-Mach-number 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 Small-scale 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.

thumbnail Fig. C.1

Vertical profiles of Vy (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 b2/3 to remove the dependence of the energy injection rate, while Βy is rescaled by the corresponding equipartition value, .

thumbnail 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 small-scale structures in the velocity field are damped by the Lorentz force, and vertical motions mainly happen in the form of large-scale flows.

thumbnail Fig. C.3

Time evolution of the maximum (solid line) and mean (dot-dashed) 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 round-off error. Although these errors accumulate in time, by the end of the simulation magnetic monopoles are still irrelevant to the dynamics of the SSD.

thumbnail 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

  1. Andrassy, R., Herwig, F., Woodward, P., & Ritter, C. 2020, MNRAS, 491, 972 [NASA ADS] [CrossRef] [Google Scholar]
  2. Andrassy, R., Higl, J., Mao, H., et al. 2022, A&A, 659, A193 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Augustson, K. C., Brun, A. S., & Toomre, J. 2016, ApJ, 829, 92 [Google Scholar]
  4. Aydemir, A. Y., & Barnes, D. C. 1985, J. Comput. Phys., 59, 108 [NASA ADS] [CrossRef] [Google Scholar]
  5. Balsara, D. S. 2004, ApJS, 151, 149 [NASA ADS] [CrossRef] [Google Scholar]
  6. Balsara, D. S., & Spicer, D. S. 1999, J. Comput. Phys., 149, 270 [NASA ADS] [CrossRef] [Google Scholar]
  7. Berberich, J. P., Chandrashekar, P., & Klingenberg, C. 2021, Comput. Fluids, 219, 104858 [Google Scholar]
  8. Brackbill, J. U., & Barnes, D. C. 1980, J. Comput. Phys., 35, 426 [NASA ADS] [CrossRef] [Google Scholar]
  9. Brandenburg, A. 2011, ApJ, 741, 92 [NASA ADS] [CrossRef] [Google Scholar]
  10. Brandenburg, A. 2014, ApJ, 791, 12 [NASA ADS] [CrossRef] [Google Scholar]
  11. Brandenburg, A., & Subramanian, K. 2005, Phys. Rep., 417, 1 [NASA ADS] [CrossRef] [Google Scholar]
  12. Brown, B. P., Browning, M. K., Brun, A. S., Miesch, M. S., & Toomre, J. 2010, ApJ, 711, 424 [Google Scholar]
  13. Browning, M. K. 2008, ApJ, 676, 1262 [Google Scholar]
  14. Browning, M. K., Miesch, M. S., Brun, A. S., & Toomre, J. 2006, ApJ, 648, L157 [Google Scholar]
  15. Brun, A. S., & Browning, M. 2017, Living Rev. Sol. Phys., 14, 4 [Google Scholar]
  16. Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073 [Google Scholar]
  17. Brun, A. S., Browning, M. K., & Toomre, J. 2005, ApJ, 629, 461 [Google Scholar]
  18. Canivete Cuissa, J. R., & Teyssier, R. 2022, A&Amp;A 664, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Cargo, P., & Gallice, G. 1997, J. Comput. Phys., 136, 446 [NASA ADS] [CrossRef] [Google Scholar]
  20. Chacón, L. 2008, Phys. Plasmas, 15, 056103 [CrossRef] [Google Scholar]
  21. Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability (Oxford University Press) [Google Scholar]
  22. Charbonneau, P. 2013, Solar and Stellar Dynamos (Springer-Verlag) [CrossRef] [Google Scholar]
  23. 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]
  24. Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174 [NASA ADS] [CrossRef] [Google Scholar]
  25. Courant, R., Friedrichs, K., & Lewy, H. 1928, Math. Ann., 100, 32 [Google Scholar]
  26. Dai, W., & Woodward, P. R. 1998, ApJ, 494, 317 [NASA ADS] [CrossRef] [Google Scholar]
  27. Dedner, A., Kemm, F., Kröner, D., et al. 2002, J. Comput. Phys., 175, 645 [Google Scholar]
  28. Dumbser, M., Balsara, D. S., Tavelli, M., & Fambri, F. 2019, Int. J. Numer. Methods Fluids, 89, 16 [NASA ADS] [CrossRef] [Google Scholar]
  29. Edelmann, P. V. F. 2014, Dissertation, Technische Universität München, Germany [Google Scholar]
  30. 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]
  31. 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]
  32. Edelmann, P. V. F., Ratnasingam, R. P., Pedersen, M. G., et al. 2019, ApJ, 876, 4 [NASA ADS] [CrossRef] [Google Scholar]
  33. Edelmann, P. V. F., Horst, L., Berberich, J. P., et al. 2021, A&Amp;A, 652, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Einfeldt, B., Roe, P. L., Munz, C. D., & Sjogreen, B. 1991, J. Comput. Phys., 92, 273 [NASA ADS] [CrossRef] [Google Scholar]
  35. Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659 [Google Scholar]
  36. Fambri, F. 2021, Int. J. Numer. Methods Fluids, 93, 3447 [NASA ADS] [CrossRef] [Google Scholar]
  37. Featherstone, N. A., & Hindman, B. 2016, ApJ, 818, 32 [NASA ADS] [CrossRef] [Google Scholar]
  38. Featherstone, N. A., Browning, M. K., Brun, A. S., & Toomre, J. 2009, ApJ, 705, 1000 [Google Scholar]
  39. Felipe, T., Khomenko, E., & Collados, M. 2010, ApJ, 719, 357 [Google Scholar]
  40. Frank, A., Jones, T. W., Ryu, D., & Gaalaas, J. B. 1996, ApJ, 460, 777 [Google Scholar]
  41. Fuchs, F. G., Mishra, S., & Risebro, N. H. 2009, J. Comput. Phys., 228, 641 [NASA ADS] [CrossRef] [Google Scholar]
  42. Gardiner, T. A., & Stone, J. M. 2005, J. Comput. Phys., 205, 509 [NASA ADS] [CrossRef] [Google Scholar]
  43. Gardiner, T. A., & Stone, J. M. 2008, J. Comput. Phys., 227, 4123 [NASA ADS] [CrossRef] [Google Scholar]
  44. Gastine, T., & Wicht, J. 2012, Icarus, 219, 428 [Google Scholar]
  45. Ghizaru, M., Charbonneau, P., & Smolarkiewicz, P. K. 2010, ApJ, 715, L133 [NASA ADS] [CrossRef] [Google Scholar]
  46. Glasser, A. H., Sovinec, C. R., Nebel, R. A., et al. 1999, Plasma Phys. Controlled Fusion, 41, A747 [Google Scholar]
  47. Glatzmaier, G. A. 1984, J. Comput. Phys., 55, 461 [NASA ADS] [CrossRef] [Google Scholar]
  48. Glatzmaier, G. A. 1985, ApJ, 291, 300 [NASA ADS] [CrossRef] [Google Scholar]
  49. Godunov, S. K., & Bohachevsky, I. 1959, Matematiceskij Sbornik, 47, 271 [Google Scholar]
  50. Harned, D. S., & Kerner, W. 1985, J. Comput. Phys., 60, 62 [NASA ADS] [CrossRef] [Google Scholar]
  51. Higl, J., Müller, E., & Weiss, A. 2021, A&Amp;A, 646, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Horst, L., Edelmann, P. V. F., Andrássy, R., et al. 2020, A&Amp;A, 641, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. 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]
  54. Hotta, H. 2017, ApJ, 843, 52 [Google Scholar]
  55. Hotta, H., Rempel, M., & Yokoyama, T. 2015, ApJ, 803, 42 [NASA ADS] [CrossRef] [Google Scholar]
  56. 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]
  57. Jardin, S. C. 2012, J. Comput. Phys., 231, 822 [NASA ADS] [CrossRef] [Google Scholar]
  58. Jones, C. A., Kuzanyan, K. M., & Mitchell, R. H. 2009, J. Fluid Mech., 634, 291 [NASA ADS] [CrossRef] [Google Scholar]
  59. Käpylä, P. J. 2011, Astron. Nachr., 332, 43 [CrossRef] [Google Scholar]
  60. Käpylä, P. J. 2019, A&Amp;A, 631, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Käpylä, P. J. 2021, A&Amp;A, 651, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Käpylä, P. J., Korpi, M. J., & Brandenburg, A. 2008, A&Amp;A, 491, 353 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2011, Astron. Nachr., 332, 883 [CrossRef] [Google Scholar]
  64. Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2012, ApJ, 755, L22 [CrossRef] [Google Scholar]
  65. Käpylä, P. J., Mantere, M. J., Cole, E., Warnecke, J., & Brandenburg, A. 2013, ApJ, 778, 41 [CrossRef] [Google Scholar]
  66. 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]
  67. 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]
  68. Kazantsev, A. P. 1968, Sov. J. Exp. Theor. Phys., 26, 1031 [NASA ADS] [Google Scholar]
  69. Khomenko, E., & Collados, M. 2006, ApJ, 653, 739 [Google Scholar]
  70. Kippenhahn, R., Weigert, A., & Weiss, A. 2013, Stellar Structure and Evolution (Springer) [Google Scholar]
  71. Kolmogorov, A. 1941, Akad. Nauk SSSR Dokl., 30, 301 [Google Scholar]
  72. Kupka, F., & Muthsam, H. J. 2017, Living Rev. Comput. Astrophys., 3, 1 [Google Scholar]
  73. Lerbinger, K., & Luciani, J. F. 1991, J. Comput. Phys., 97, 444 [NASA ADS] [CrossRef] [Google Scholar]
  74. LeVeque, R. J. 2002, Finite Volume Methods for Hyperbolic Problems, Cambridge Texts in Applied Mathematics (Cambridge University Press) [Google Scholar]
  75. Londrillo, P., & del Zanna, L. 2004, J. Comput. Phys., 195, 17 [NASA ADS] [CrossRef] [Google Scholar]
  76. Lütjens, H., & Luciani, J.-F. 2010, J. Comput. Phys., 229, 8130 [CrossRef] [Google Scholar]
  77. Masada, Y., Yamada, K., & Kageyama, A. 2013, ApJ, 778, 11 [NASA ADS] [CrossRef] [Google Scholar]
  78. Matthaeus, W. H., & Brown, M. R. 1988, Phys. Fluids, 31, 3634 [NASA ADS] [CrossRef] [Google Scholar]
  79. Meneguzzi, M., Frisch, U., & Pouquet, A. 1981, Phys. Rev. Lett., 47, 1060 [NASA ADS] [CrossRef] [Google Scholar]
  80. Mestel, L. 1999, Stellar Magnetism (Oxford University Press) [Google Scholar]
  81. Miczek, F. 2013, Dissertation, Technische Universität München, Germany [Google Scholar]
  82. Miczek, F., Röpke, F. K., & Edelmann, P. V. F. 2015, A&Amp;A, 576, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  83. Mignone, A., & Del Zanna, L. 2021, J. Comput. Phys., 424, 109748 [NASA ADS] [CrossRef] [Google Scholar]
  84. Minoshima, T., & Miyoshi, T. 2021, J. Comput. Phys., 446, 110639 [CrossRef] [Google Scholar]
  85. Minoshima, T., Kitamura, K., & Miyoshi, T. 2020, ApJS, 248, 12 [NASA ADS] [CrossRef] [Google Scholar]
  86. Miyoshi, T., & Kusano, K. 2005, J. Comput. Phys., 208, 315 [NASA ADS] [CrossRef] [Google Scholar]
  87. Müller, B. 2020, Living Rev. Comput. Astrophys., 6, 3 [CrossRef] [Google Scholar]
  88. Pietarila Graham, J., Cameron, R., & Schüssler, M. 2010, ApJ, 714, 1606 [Google Scholar]
  89. Powell, K. G. 1997, in Upwind and High-Resolution Schemes, eds. M.Y. Hussaini, B. van Leer, J. Van Rosendale (Springer), 570 [CrossRef] [Google Scholar]
  90. 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]
  91. Rempel, M. 2005, ApJ, 622, 1320 [NASA ADS] [CrossRef] [Google Scholar]
  92. Rempel, M. 2014, ApJ, 789, 132 [Google Scholar]
  93. Rempel, M. 2018, ApJ, 859, 161 [Google Scholar]
  94. Riva, F., & Steiner, O. 2022, A&Amp;A, 660, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Rogers, T. M., Lin, D. N. C., McElwaine, J. N., & Lau, H. H. B. 2013, ApJ, 772, 21 [NASA ADS] [CrossRef] [Google Scholar]
  96. Schekochihin, A. A., Cowley, S. C., Taylor, S. F., Maron, J. L., & McWilliams, J. C. 2004, ApJ, 612, 276 [NASA ADS] [CrossRef] [Google Scholar]
  97. Schekochihin, A. A., Iskakov, A. B., Cowley, S. C., et al. 2007, N. J. Phys., 9, 300 [NASA ADS] [CrossRef] [Google Scholar]
  98. Schnack, D., Barnes, D., Mikic, Z., Harned, D. S., & Caramana, E. 1987, J. Comput. Phys., 70, 330 [NASA ADS] [CrossRef] [Google Scholar]
  99. Seta, A., & Federrath, C. 2020, MNRAS, 499, 2076 [NASA ADS] [CrossRef] [Google Scholar]
  100. Shu, C.-W., & Osher, S. 1988, J. Comput. Phys., 77, 439 [NASA ADS] [CrossRef] [Google Scholar]
  101. Smolarkiewicz, P. K., & Charbonneau, P. 2013, J. Comput. Phys., 236, 608 [NASA ADS] [CrossRef] [Google Scholar]
  102. Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137 [NASA ADS] [CrossRef] [Google Scholar]
  103. Strang, G. 1968, SIAM J. Numer. Anal., 5, 506 [CrossRef] [MathSciNet] [Google Scholar]
  104. Sutherland, B. R. 2010, Internal Gravity Waves (Cambridge University Press) [Google Scholar]
  105. Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501 [NASA ADS] [CrossRef] [Google Scholar]
  106. Toro, E. F. 2009, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction (Springer) [CrossRef] [Google Scholar]
  107. Tóth, G. 2000, J. Comput. Phys., 161, 605 [Google Scholar]
  108. Viallet, M., Baraffe, I., & Walder, R. 2011, A&Amp;A, 531, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  109. Viviani, M., Käpylä, M. J., Warnecke, J., Käpylä, P. J., & Rheinhardt, M. 2019, ApJ, 886, 21 [NASA ADS] [CrossRef] [Google Scholar]
  110. Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&Amp;A, 429, 335 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  111. 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]
  112. Yadav, R. K., Christensen, U. R., Wolk, S. J., & Poppenhaeger, K. 2016, ApJ, 833, L28 [Google Scholar]

1

Low-β environments can be found in the outer layers of active stars, like the solar corona.

2

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).

3

However, other source terms can be added to the system depending on the problem at hand. These include energy generation by nuclear reactions, radiative transport of energy in the diffusion limit, neutrino cooling, and parabolic viscous terms.

4

Throughout the paper we use the Lorentz-Heaviside notation: B = b/ .

5

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 energy-conservation properties in simulations of gas dynamics with gravity (Müller 2020; Edelmann et al. 2021).

6

Here, x represents a generic direction.

7

For simplicity, we only consider the homogeneous MHD system in this analysis.

8

Here we describe the 3D algorithm; however, 1D and 2D Cartesian grids can also be used in SLH.

9

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.

10

The physical fluxes are usually computed in Riemann solvers as some variation of the central flux (FL + FR)/2.

11

Deviations in the primitive variables can also be reconstructed if the corresponding equilibrium values are provided at the cell centers and at the cell boundaries.

12

Here the calculation is made over the cell boundary (i + 1/2, j, k).

13

For more details on the implementation of implicit time stepping in SLH, see Miczek (2013) and Miczek et al. (2015).

14

We took these values as representative of the typical conditions found in stellar convection zones close to equipartition regimes (Augustson et al. 2016).

15

For the time-explicit simulations, the CFL time step is reduced by 20%.

16

No snapshots were saved throughout the simulations to minimize the cost of I/O operations.

17

Better entropy-conservation properties can be achieved by directly evolving the specific entropy instead of ρEϕ. However, this approach does not conserve the total energy.

18

More details on how to compute the pressure profile can be found in Edelmann et al. (2021).

19

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.

20

Details on how to compute the density perturbation can be found in Andrassy et al. (2022).

All Tables

Table 1

Time averages of EM/b2/3 [1045 × erg] over the time interval 20 < t/τconv < 70 for the different resolutions, N, and boost factors, b, considered in this study.

All Figures

thumbnail Fig. 1

Wave structure of the MHD system.

In the text
thumbnail Fig. 2

Wave structure of HLLD-type solvers: only the fast magnetosonic (SL, SR), Alfvén , and entropy (SM) waves are considered in the computation of the states across the Riemann fan. The slow magnetosonic waves are discarded.

In the text
thumbnail Fig. 3

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

In the text
thumbnail Fig. 4

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

In the text
thumbnail Fig. 5

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

In the text
thumbnail Fig. 6

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

In the text
thumbnail Fig. 7

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

In the text
thumbnail Fig. 8

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

In the text
thumbnail Fig. 9

Initial setups of Vx and Vy (here rescaled by ℳx) used for simulating the growth of the Kelvin–Helmholtz instability.

In the text
thumbnail Fig. 10

Time evolution of the y-direction 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 dot-dashed 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
thumbnail Fig. 11

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

In the text
thumbnail Fig. 12

Distribution of the sonic Mach number in the Kelvin-Helmholtz instability test at t/tmax =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 Mx.

In the text
thumbnail Fig. 13

Convergence with resolution N of the L1 error associated with EK,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 (dot-dashed lines) solvers. The dashed black line is the second-order scaling.

In the text
thumbnail 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
thumbnail 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 well-balancing, whereas no well-balancing 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
thumbnail 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, dot-dashed: 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 b2/3 to take into account the different energy contents of the flows.

In the text
thumbnail 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 = Vy/a multiplied by 103, while the plots on the right show the vertical magnetic field rescaled by the root mean square value across the plane.

In the text
thumbnail 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 b1/3 scaling.

In the text
thumbnail 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 b2/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 (k3/2) scalings.

In the text
thumbnail Fig. A.1

Distribution of the sonic Mach number in the simulations of the Kelvin-Helmholtz instability at t/tmax = 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 Mx = 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
thumbnail Fig. A.2

Distribution of Βx in the simulations of the Kelvin–Helmholtz instability at t/tmax = 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 grid-scale 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
thumbnail Fig. A.3

Time evolution of Bx 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
thumbnail 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
thumbnail Fig. B.2

Final distribution of ∣B/Bx,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
thumbnail 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 low-Mach-number 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
thumbnail Fig. C.1

Vertical profiles of Vy (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 b2/3 to remove the dependence of the energy injection rate, while Βy is rescaled by the corresponding equipartition value, .

In the text
thumbnail 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 small-scale structures in the velocity field are damped by the Lorentz force, and vertical motions mainly happen in the form of large-scale flows.

In the text
thumbnail Fig. C.3

Time evolution of the maximum (solid line) and mean (dot-dashed) 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 round-off 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
thumbnail 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 (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.