© J. R. Canivete Cuissa and R. Teyssier 2022

1. Introduction

Thermal turbulent convection is one of the fundamental processes in stellar physics. It is responsible for the outward transport of the energy generated in the core, but it also affects the structure, the dynamics, and the evolution of the star. In our Sun, for example, the convective envelope shapes and influences the observable features and activity on the solar surface and in the overlying atmosphere (Nordlund et al. 2009; Stein 2012). Core convection in massive stars might impact the star’s lifetime by bringing fresh fuel into the core as the convective cells overshoot into the stable radiative zones above it (Salaris & Cassisi 2017), and convective mixing in asymptotic giant branch (AGB) stars provides a rich environment for nucleosynthesis (Herwig 2005).

Stellar magnetism is also linked to turbulent convective motions in stellar interiors. A small-scale dynamo operating in the near-surface turbulent convection zone of the Sun is possibly at the origin of the small-scale magnetic fields permeating the quiet photosphere (Lites et al. 2014; Rempel 2018), while the cyclic regeneration of the solar large-scale magnetic field probably stems from the interplay between deep-convection zone turbulence, differential rotation, and magnetic flux transport (Charbonneau 2013, 2020). Moreover, convective cores in A- and B-type stars are likely able to develop a dynamo action (Brun et al. 2005; Augustson et al. 2016) and all low-mass stars appear to host dynamo-produced, small-scale, surface magnetic fields (Langer 2014).

Finally, turbulent convective envelopes and cores also cause the excitation and propagation of a rich spectrum of oscillation modes in stellar interiors (Houdek & Dupret 2015). With the advent of asteroseismology, it is possible to determine properties of stellar internal structures from the observation of global oscillations across the Herzsprung-Russell diagram (Hekker & Christensen-Dalsgaard 2017; Basu & Hekker 2020; Aerts 2021). It is therefore crucial to understand and characterize the formation and propagation of these modes, in particular in the presence of magnetic fields.

The classical modeling of convection is based on mixing length theory (MLT, Böhm-Vitense 1958). However, large-scale (magneto-)hydrodynamical simulations of solar and stellar convection are nowadays employed to include multidimensional dynamical processes such as overshooting, oscillations, and dynamo action (see, e.g., Freytag et al. 2002, 2017; Browning et al. 2004, 2006; Brun et al. 2004; Herwig et al. 2006; Meakin & Arnett 2007; Miesch et al. 2008; Brun & Palacios 2009; Ghizaru et al. 2010; Rempel 2014; Hotta et al. 2014, 2015; Beeck et al. 2015; Tremblay et al. 2015; Augustson et al. 2016; Käpylä et al. 2017; Salhab et al. 2018; Edelmann et al. 2019; Andrassy et al. 2020; Horst et al. 2020, 2021; Käpylä 2021).

Such simulations can be challenging for common numerical techniques as stellar convective flows are often highly subsonic and characterized by small perturbations around a hydrostatic equilibrium. In order to alleviate the time-stepping constraints imposed by the low-Mach number regime, numerical schemes solving the Navier-Stokes (or MHD) equations in the anelastic approximation are typically favored (see Kupka & Muthsam 2017, for a review).

Although internal gravity waves are preserved in this approach, the physics of pressure waves are precluded and an artificial viscosity is required to achieve numerical stability. Fully compressible simulations are therefore necessary to examine the properties and the dynamics of the complete spectrum of excited waves in stellar interiors, as well as the coupling between these modes (Beck et al. 2011). The computational cost of these simulations is higher than the anelastic ones, but Horst et al. (2020, 2021) recently demonstrated the feasibility of this approach in two and three dimensions.

The next logical step is to extend the fully compressible framework to magnetized stars. Such simulations would allow to study the interaction between convective turbulent motions, magnetic fields, and global oscillations. In particular, it would enable a numerical analysis of the imprints of magnetic fields in global oscillations, which could be used to assess the presence of strong magnetic fields hidden in stellar interiors from asteroseismologic measurements (Fuller et al. 2015; Gomes & Lopes 2020).

Therefore, we test the feasibility of three-dimensional, fully compressible, magneto-hydrodynamical numerical simulations of stellar convection. This work is hence a proof of concept: we intend to present and validate our approach by performing a convergence study and comparing qualitatively the results to previous works on stellar convection, dynamo action, and wave propagation. The physics and the numerical setup are therefore intentionally simplistic, and they will be enhanced and addressed in future works.

This paper is organized as follows: in Sect. 2 we briefly describe the code, the numerical techniques, and define the initial conditions of our simulations. The results are presented and discussed in Sect. 3, while in Sect. 4 we summarize the main aspects of our work and we give an outlook for future developments.

2. Numerics

We ran our numerical simulations with the Adaptive Mesh Refinement (AMR) code RAMSES (Teyssier 2002), which solves the compressible equations of ideal magneto-hydrodynamics (MHD) in presence of gravity by employing a MUSCL-Hancock scheme with constrained transport on a finite volume Cartesian mesh (Teyssier et al. 2006; Fromang et al. 2006). We used a HLLC and a HLLD approximate Riemann solver for purely hydrodynamical and MHD simulations, respectively. The advantage of this type of solvers is that they resolve contact discontinuities explicitly, which is particularly important in highly subsonic flows to reduce numerical diffusivity. RAMSES makes intensive use of the Message Passing Interface (MPI) library and it can therefore be used on massively parallel architectures. It is designed for high-resolution numerical simulations of a wide range of astrophysical problems, such as cosmology, galaxy and structure formation and evolution, and star formation.

However, the features of convective flows in stars pose a number of supplementary challenges that numerical schemes need to overcome. The pressure gradients in stellar interiors are balanced by gravity so that the whole system can attain a stationary configuration, that is hydrostatic equilibrium. Turbulent thermal convection is then often characterized by small-amplitude perturbations close to the equilibrium profile. As a consequence, we sought a numerical scheme that is able to preserve the hydrostatic equilibrium profiles in nonperturbed setups. Moreover, the typical velocities arising from turbulent convection are expected to be highly subsonic, which can be troublesome for time-explicit schemes. Therefore, we had to adapt the numerical scheme of RAMSES to deal with highly subsonic, close to hydrostatic equilibrium turbulent flows, and we achieved that by implementing a well-balanced method.

2.1. Well-balanced scheme

A well-balanced scheme numerically ensures the dynamical preservation of an equilibrium state by implicitly including it in the underlying discrete set of equations (Greenberg & Leroux 1996). Such methods have been mainly developed in the context of shallow-water simulations (see, e.g., Noelle et al. 2007), but they are also becoming popular to tackle astrophysical problems (see, e.g., Freytag et al. 2012; Käppeli & Mishra 2016; Veiga et al. 2019; Edelmann et al. 2021).

For the sake of simplicity, let us consider the one-dimensional Euler-Poisson set of equations,

$$\begin{array}{cc}\hfill \frac{\partial \rho }{\partial t}+\frac{\partial (\rho u)}{\partial x}& =0,\hfill \\ \hfill \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+\frac{1}{\rho }\frac{\partial p}{\partial x}& =-\frac{\partial \varphi }{\partial x},\hfill \\ \hfill \frac{\partial p}{\partial t}+u\frac{\partial p}{\partial x}+\gamma p\frac{\partial u}{\partial x}& =0,\hfill \end{array}$$(1)

where ρ is the density, u is the velocity, p is the thermal pressure, ϕ is the gravitational potential, and γ is the adiabatic index of the ideal gas equation of state. The hydrostatic equilibrium equation then reads,

$$\begin{array}{c}\hfill \frac{\partial p}{\partial x}=-\rho \frac{\partial \varphi }{\partial x}·\end{array}$$(2)

Standard finite volume methods struggle to accurately describe convective flows in stratified fluids because a discrete version of Eq. (2) may not be preserved. Hence, truncation errors introduced in the supposedly equilibrium states are interpreted as propagating waves by the numerical solvers. These spurious waves can conceal the small perturbations of interest, disrupt the equilibrium states, and even undermine the numerical stability of the simulation. The numerical resolution required to prevent the truncation errors from spoiling the results would make the simulations easily impractical.

A possible solution is to impose the hydrostatic equilibrium given by Eq. (2) directly in the set of equations governing the dynamics, that is Eq. (1). We separate the primitive variables of the problem into equilibrium, stationary states and dynamical perturbations,

$$\begin{array}{c}\hfill \rho =\stackrel{\sim }{\rho }+{\rho }^{\prime },p=\stackrel{\sim }{p}+p^{\prime },u=u^{\prime },\end{array}$$(3)

where the equilibrium states of density and pressure, $\stackrel{\sim }{\rho }$ and $\stackrel{\sim }{p}$, exactly satisfy the hydrostatic equilibrium, Eq. (2), while the velocity field u is made up of perturbations only since $\stackrel{\sim }{u}=0$. With this separation at hand, we can rewrite the Euler-Poisson equations as,

$$\begin{array}{cc}\hfill \frac{\partial {\rho }^{\prime }}{\partial t}+\frac{\partial \left(\rho u\right)}{\partial x}& =0,\hfill \\ \hfill \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+\frac{1}{\rho }\frac{\partial p^{\prime }}{\partial x}& =-\frac{{\rho }^{\prime }}{\rho }\frac{\partial \varphi }{\partial x},\hfill \\ \hfill \frac{\partial p^{\prime }}{\partial t}+u\frac{\partial p}{\partial x}+\gamma p\frac{\partial u}{\partial x}& =0,\hfill \end{array}$$(4)

where the hydrostatic equilibrium in stationary regimes, that is where ρ′,p′,u = 0, is explicitly preserved.

Well-balanced schemes have already been used in the context of stellar convection by Hotta et al. (2019), Horst et al. (2020, 2021), and Edelmann et al. (2021), for example. For this work, we modified both the hydrodynamical and MHD numerical solvers of RAMSES to account for a well-balanced version of the evolution equations. Indeed, we extended the well-balanced methodology also to the equations of ideal MHD, where we assume the equilibrium state to be characterized by $\stackrel{\sim }{\mathit{B}}=0$. Hence, magnetic fields are treated as pure perturbations, B = B′, and they do not modify the initial hydrostatic equilibrium, Eq. (2). Moreover, we solve for the conservation equation of specific entropy s = p/ρ^γ instead of total energy, that is,

$$\begin{array}{c}\hfill \frac{\partial s}{\partial t}+u\frac{\partial s}{\partial x}=0.\end{array}$$(5)

As a consequence, the total energy of the system may not be exactly conserved. On the other hand, we can accurately follow the dynamics of entropy perturbations which stem from the convective flows in adiabatic stratifications. The validity of this approach is guaranteed by the highly subsonic nature of typical stellar convective flows, where shocks are essentially absent (Woodward et al. 2015).

2.2. Simulation setup

We aim to test our code on a simple but realistic, plane-parallel, stellar convective region. We chose to follow the initial setup presented by Herwig et al. (2006), where an intershell of a typical low-mass AGB star near the He shell flash luminosity peak is modeled. The simulation box covers a cubic domain of dimensions x, y, z ∈ [0, L], where L = 11.0 Mm. This box is located at a radius of 7.51 Mm of a stellar model, with main-sequence initial mass of 2 M_⊙ and metallicity Z = 0.01, which undergoes its second-to-last thermal pulse (Herwig & Austin 2004).

The box is divided vertically in three domains which approximate the stellar stratification: a bottom stable region that extends up to z = 1.64 Mm, a convection zone in the middle, and an upper stable region from z = 7.77 Mm to the top boundary. We initialized each domain with a polytropic stratification in hydrostatic equilibrium. The properties of the stratification at the bottom of each domain, as well as the respective polytropic indices, are given in Table 1. The corresponding equilibrium profiles for density, pressure, and specific entropy are shown in Fig. 1. We also assumed constant gravity directed downward along the vertical axis (z) with value g = 10^7.7 cm s⁻² and a mono-atomic ideal gas with adiabatic index γ = 5/3 and equation of state

$$\begin{array}{c}\hfill p=(\gamma -1)\rho (e-\frac{1}{2}{u}^2),\end{array}$$(6)

Fig. 1. Initial and equilibrium profiles of density $\stackrel{\sim }{\rho }$, pressure $\stackrel{\sim }{p}$, and specific entropy $\stackrel{\sim }{s}$. Black dotted lines delimit the convection zone, while the gray areas represent the heating (bottom) and cooling (top) layers.

where e is the specific total energy. For more information about the stellar model, the reader can refer to Herwig et al. (2006).

If the simulations are initialized with the equilibrium profiles of Fig. 1, the code is able to preserve the initial state indefinitely up to machine precision thanks to the implemented well-balanced scheme. Hence, to set off convection, we introduced random density perturbations δρ of the order $\delta \rho /\stackrel{\sim }{\rho }\sim {10}^{-2}$ in every grid cell within 1 Mm from the bottom of the convection zone.

Turbulent convective motions are then sustained by nuclear reactions heating at the bottom and radiative cooling at the top, which we modeled by adding and removing energy in two bands of thickness Δz = 0.5 Mm at the bottom and top, respectively, as shown in gray in Fig. 1. We note that Herwig et al. (2006) did not include radiative cooling, therefore our setup resembles that of surface convection. We assumed constant and equal volumetric heating and cooling rates given by $\dot{e}={\dot{ϵ}}_0{\rho }_0$, where ${\dot{ϵ}}_0=2\times {10}^{10}\mathrm{erg}{\mathrm{g}}^{-1}{\mathrm{s}}^{-1}$ and ρ₀ is the density at the base of the convection zone given in Table 1. Consequently, the energy flux at the bottom of the convection zone is F = ėΔz = 1.17 × 10²² erg s⁻¹ cm⁻². This value accounts for the integrated amount of energy released according to the stellar model and corresponds to a stellar luminosity of L = 3.21 × 10⁷ L_⊙. We note that the stellar luminosity is often boosted by several orders of magnitude in anelastic numerical simulations to balance the effect of high artificial viscosities (see, e.g., Rogers et al. 2013). Since we are using a fully compressible approach, we can stick to more realistic values of L.

We ran hydrodynamical (HD) and magneto-hydrodynamical (MHD) numerical simulations of the convective envelope described above. In both cases, we used the well-balanced MHD solver of RAMSES to allow for a direct comparison between the two. The hydrodynamical simulations are initialized with B = (0, 0, 0) G, ensuring that the magnetic field strength will remain zero everywhere throughout the simulation. For the magneto-hydrodynamical cases, we set a constant, homogeneous, and horizontal (along the x axis) magnetic field of strength B₀ = 10³ G in the convection zone alone. In the stable domains instead, B = (0, 0, 0) G. This configuration preserves exactly the constraint ∇ ⋅ B = 0. Moreover, the initial magnetic field is weak enough for a turbulent small-scale dynamo in the convection zone to amplify it significantly, as the expected equipartition value of the magnetic field strength for this setup is B_eq ∼ 10⁸ G. We adopted fixed boundary conditions in the vertical direction by forcing the hydrodynamical variables to follow the polytropic equilibrium profiles, while the magnetic fields are forced to be zero in the ghost cells. The lateral boundary conditions are periodic for all variables. We ran both hydrodynamical and magneto-hydrodynamical simulations with resolutions N = 64³,  128³,  256³, and 512³. We identify the different simulations by their type and resolution (e.g., HD_256). In this study, we did not use the AMR capabilities of RAMSES. Therefore the grid is equidistant in each direction with cell sizes ranging from 21 km to 172 km.

3. Results and discussion

In this section, we present and discuss the results of our simulations. We focus on the onset of convection, on the turbulent amplification of the magnetic energy, and on the general properties of the convective flows and magnetic fields. Finally, we investigate the propagation of pressure and internal gravity waves in the stable regions.

3.1. Onset of convection

We start by analyzing the onset of convection in the hydrodynamical simulations, HD_64, HD_128, HD_256, and HD_512, which we ran for a total of 2000 s physical time. We define the convective turnover time scale as,

$$\begin{array}{c}\hfill {\tau }_{\mathrm{conv}}=\frac{2L_{\mathrm{CZ}}}{v_{\mathrm{conv}}}\sim 1000\mathrm{s},\end{array}$$(7)

where L_CZ ∼ 6 Mm is the height of the convection zone and v_conv ∼ 1.2 × 10⁶ cm s⁻¹ is the characteristic convective velocity (see Sect. 3.3). Hence, our hydrodynamical simulations cover ∼2 convective turnover times.

Figure 2 shows three vertical sections of the temperature fluctuations, defined as $T^{\prime }=T-\stackrel{\sim }{T}$, at times t = 100,  400, and 1000 s, respectively, for the HD_512 simulation. Small-scale thermal convective instabilities developing from the initially perturbed cells can be seen at time t = 100 s. These finger-like structures rise through the convection zone at slightly higher temperatures than the equilibrium values, on the order of ∼0.1%. The growth of these instabilities also excite small amplitude pressure waves that sweep the convection zone and reach the top stable layer. These waves induce tiny perturbations in the equilibrium profiles, that are not visible in the left panel of Fig. 2 but seed a second cascade of thermal convective instabilities. This cascade is driven by the external radiative cooling layer, as we can see in the middle panel of the same figure, and it is characterized by descending turbulent filaments with negative temperature fluctuations around the equilibrium profile (see, e.g., Viallet et al. 2013).

Fig. 2. Temperature fluctuations T′ in a vertical section of the HD_512 simulation at t = 100 s (left), t = 400 s (middle), and t = 1000 s (right). The vertical section is taken at y = 3.43 Mm. Black dotted lines indicate the boundaries of the convection zone.

After t ∼ 1000 s, we attain a fully developed convection zone, which is maintained thereafter by the balance between external cooling and heating. Indeed, when the rising plasma heated at the bottom reaches the top of the convection zone, it is cooled down and initiates a turbulent downflow (see right panel of Fig. 2). The large-scale structure of the flow is therefore granular-like, with large and hot slowly uprising plumes surrounded by narrow and cold filamentary downflows.

In Fig. 3 we show a horizontal section of the vertical velocity v_z and of the temperature fluctuations T′ near the top of the convection zone at t = 1000 s. A granulation pattern, typical of surface convection (see, e.g., Stein & Nordlund 1998), is clearly observable with the central granule having an approximate diameter of ∼5 Mm. The large-scale convective motions buffet the stable plasma in the bottom and top stable layers, exciting there different oscillation modes that are visible as horizontally extending temperature fluctuations in the central and right panels of Fig. 2. We address the nature of these oscillations in Sect. 3.4.

Fig. 3. Horizontal sections of the vertical velocity vz (top) and temperature fluctuation T′ (bottom) in the upper layers of the convection zone (z = 7.4 Mm). The snapshot is taken from the HD_512 simulation at t = 1000 s.

It is important to notice that the dynamical fluctuations characterizing the convective region are small if compared to the equilibrium profiles values. Figures 2 and 3 show temperature fluctuations with typical values around 3 orders of magnitude smaller than the equilibrium values in the convection zone. Similar amplitudes are also valid for the density and pressure perturbations, while the typical vertical velocities at the top of the convection zone are 2 orders of magnitude smaller than the local sound speed c_s. Therefore, the well-balanced scheme implemented in the code and presented in Sect. 2.1 is crucial to correctly capture the dynamics of this problem.

Figure 4 shows, for the four hydrodynamical simulations, the time evolution of the mean turbulent kinetic energy density, $E_{\mathrm{K}}=\frac{1}{2}\rho v_{\mathrm{rms}}^2$, where v_rms is the root-mean-squared (rms) turbulent three-dimensional velocity field strength. Since the two stable layers present no significant turbulent motions, we restrict the computation of E_K to the convection zone alone. The onset of convection is represented by an initial exponential growth, which lasts between ∼250 s and ∼1000 s, depending on the resolution. Then, a quasi-steady state is reached, characterized by a quasi-constant mean turbulent kinetic energy density. The runs HD_128, HD_256, and HD_512 reach a mean turbulent kinetic energy density of around 10¹⁶ erg cm⁻³, while the low resolution one, HD_64, stabilizes around a slightly lower value. We notice long lived fluctuations in the quasi-steady state mean kinetic energy density that are due to the episodic and large-scale nature of the flow in the convection zone (Meakin & Arnett 2007).

Fig. 4. Time evolution of the mean turbulent kinetic energy density EK in the convection zone for the four hydrodynamical (HD) simulations.

The growth rate of the mean kinetic energy density seems to decrease with increasing resolution: the low-resolution simulations (HD_64 and HD_128) show a fast and sharp evolution to a quasi-steady state, while a smooth and slower growth characterizes the high-resolution ones (HD_256 and HD_512). We can qualitatively explain this difference with the transition between initial instabilities and large-scale convective flows. In fact, the thermal convective instabilities at the origin of convection form at the smallest scales in the simulation, since the random perturbations in density are introduced at the grid-size level. On the other hand, the kinetic energy dominant structures are the large-scale ones, that is plumes and downflows. This results in an apparent slower growth since the smaller the initial scale of the turbulence, the longer it will take to form steady state, large-scale convection.

The transition process between small and large scales can also be observed in Fig. 5, where we show the time evolution of the turbulent kinetic energy density power spectra, ${\widehat{E}}_{\mathrm{K}}$, for the four hydrodynamical simulations. We notice that in the first stages of the simulations (light blue), the peaks of the various power spectra are found at large wavenumber k, that is at small scales in physical space. As the simulations evolve, the total kinetic energy densities grow and the peaks of the power spectra shifts toward small wavenumbers.

Fig. 5. Time evolution of the turbulent kinetic energy power spectra ${\widehat{E}}_{\mathrm{K}}$ for the four HD simulations. The time evolution is represented by the color grading: light blues correspond to early times, while the darkest shade corresponds to t = 2000 s. The interval between each color-shade is 50 s.

When the quasi-steady state is reached (dark blue profiles), the energy cascade follows a Kolmogorov power law with index −5/3 from the steady state energy peak up to the dissipation scale in all panels of Fig. 5. The same can be observed in Fig. 6, where we compare the average kinetic energy power spectra during the quasi-steady state phase between the different runs. The better the spatial resolution, the larger the range where the Kolmogorov cascade is well reproduced. The quasi-steady state energy peak instead is found around k ∼ 1 Mm⁻¹ for all resolutions. This value corresponds to the typical size of the large-scale uprising convective plumes observed in Fig. 3, that is l = 2π/k ∼ 6 Mm.

Fig. 6. Comparison of the kinetic energy power spectra ${\widehat{E}}_{\mathrm{K}}$ between the four HD simulations in the quasi-steady state regime. The profiles are obtained by averaging the power spectrum profiles over the last 10 snapshots of each simulation, that is between t  =  1500 s and t  =  2000 s.

3.2. Small-scale dynamo

In this section we present the magneto-hydrodynamical simulations MHD_64, MHD_128, MHD_256, and MHD_512.

Figure 7 shows the time evolution of the mean turbulent kinetic energy density E_K in blue, and of the mean magnetic energy density E_M = B²/8π in red, for the four simulations. Both energy densities are computed in the convection zone alone. The length of the simulations decreases with resolution because of the increase of computational cost. Nevertheless, we ran each simulation long enough for the magnetic field to be amplified and attain saturation. Indeed, while convection develops and reaches a quasi-steady state in a much shorter time scale, we observe in all four simulations a slower kinematic and a saturation phase: first the magnetic energy grows exponentially as E_M ∼ exp(γ_Kt), where γ_K is the kinematic growth rate, and then it reaches a quasi-constant value.

Fig. 7. Time evolution of the mean turbulent kinetic energy density EK (blue) and of the mean magnetic energy density EM (red) in the convection zone for the four magneto-hydrodynamical (MHD) simulations. Dashed lines represent exponential fits obtained with the kinematic growth rates γK listed in Table 2. The top panel is a zoom-in on the first 20 000 s of the turbulent kinetic energy evolution.

The different slopes show how the exponential growth depends on resolution, which is compatible with the action of a small-scale dynamo. In the first row of Table 2 we show the kinematic growth rates γ_K of the different simulations. We find that the kinematic growth rate scales with resolution as γ_K ∼ Δx^−1.2 if we include all four simulations. However, the low-resolution simulation (MHD_64) is not capturing all the dynamical scales of convection, as we have seen in Sect. 3.1. If we consider only runs MHD_128, MHD_256, and MHD_512, the growth rate is γ_K ∼ Δx^−1.3. This result, which is consistent with what Pietarila Graham et al. (2010), Rempel (2014), and Riva & Steiner (2022) found for small-scale dynamo simulations of the quiet Sun magnetism, is very puzzling and in contrast with Kazantsev’s dynamo theory prediction of γ_K ∼ Δx^−2/3.

The kinematic phase comes to an end when the magnetic field strength approaches the equipartition value, $B_{\mathrm{eq}}=\sqrt{4\pi \rho }{v}_{\mathrm{rms}}$. At this stage, the magnetic fields are strong enough to start back-reacting into the plasma dynamics by means of the Lorentz force. The mean magnetic energy eventually saturates and each simulation reaches a magneto-convective quasi-steady state. Typical magnetic field strengths in the convection zone at this stage are in the order of B ∼ 10⁸ G. The ratios between the rms magnetic field strength B_rms and the equipartition magnetic field strength B_eq during the saturation phase are shown in Table 2 and grow with increasing spatial resolution.

We show the time evolution of the turbulent kinetic and magnetic energy power spectra in Fig. 8. We notice that the magnetic field strength is more efficiently amplified at small scales. During the kinematic phase (light red) in fact, the larger the spatial resolution, the more the peak of the magnetic energy power spectrum is shifted toward large k. In an ideal Kazantsev’s dynamo, the magnetic energy power spectra should follow a power-law with index 3/2 during the kinematic phase (Brandenburg & Subramanian 2005). In our case, the magnetic energy power spectra follow the Kazantsev’s prediction only during the very early times of the kinematic phase. However, as they approach the saturation phase, they do not seem to follow such a profile.

Fig. 8. Time evolution of the power spectra of the turbulent kinetic energy ${\widehat{E}}_{\mathrm{K}}$ (blue) and magnetic energy ${\widehat{E}}_{\mathrm{M}}$ (red) for the four magneto-hydrodynamical (MHD) simulations. The time evolution is represented by the color grading: light colors correspond to early times, while the darkest shade corresponds to the last snapshot of the simulation. The time interval between each shade is 5000 s for MHD_64, 1000 s for MHD_128, 500 s for MHD_256, and 300 s for MHD_512. The average quasi-steady state kinetic energy power spectrum for the hydrodynamical simulation HD_512 is shown in green.

As the simulations evolve, the magnetic fields become dominant over the kinetic energy density at large k for all the resolutions except the MHD_64. The crossover scale is found around k ∼ 10 Mm⁻¹ (l ∼ 2 Mm). The transition to super-equipartition of the magnetic fields at large k coincides with a suppression of kinetic power at the same scales, which follows from the Lorentz-force feedback. Consequently, the kinetic power spectra profiles deviate from the Kolmogorov power law for k ≳ 3 Mm⁻¹ (l ∼ 0.5 Mm). We can appreciate the difference to the hydrodynamical simulation HD_512 power spectrum which is plotted in green in Fig. 5. Approaching the saturation phase (dark red), the magnetic energy spectrum peaks shift toward smaller wavenumbers and finally stabilizes around k ∼ 3 Mm⁻¹, that corresponds to magnetic structures in physical space of size l ∼ 0.5 Mm.

In Fig. 9 we show the kinetic and magnetic power spectra for the last snapshot of the four simulations. We notice that MHD_128, MHD_256, and MHD_512 yield very similar results in the respective energy injection and inertial ranges, while the low-resolution simulation, MHD_64, presents lower profiles for both kinetic and magnetic spectra at all scales. The kinetic energy power spectra of the MHD simulations clearly deviate from the hydrodynamical analog, HD_512, shown in green. The suppression of kinetic energy is evident at around the injection scale (k ∼ 1 Mm⁻¹) and the crossover scale (k ∼ 10 Mm⁻¹). The power spectra evolution and steady state configuration for the high-resolution simulations are in qualitative agreement with the results obtained by Rempel (2014) and Hotta et al. (2015) for small-scale dynamo studies in the solar convection zone. Therefore, the exponential growth of the magnetic field energy in our simulations is compatible with the action of a turbulent small-scale dynamo.

Fig. 9. Comparison between the power spectra of the turbulent kinetic energy (blue) and magnetic energy (red), ${\widehat{E}}_{\mathrm{K}}$ and ${\widehat{E}}_{\mathrm{M}}$, for the four magneto-hydrodynamical (MHD) simulations after saturation of the magnetic field. The average quasi-steady state kinetic energy power spectrum for the hydrodynamical simulation HD_512 is shown in green.

In Fig. 10 we show a vertical section of the temperature fluctuations T′, vertical velocity v_z, and vertical magnetic field B_z for the MHD_512 simulation in the saturated phase. To enhance the visibility of the perturbations, we scaled each quantity by their mean rms value at each height z. However, in the stable regions where the magnetic field is very weak, spurious signals arise where B_z, rms ∼ 0 G. Therefore we masked out the regions where B_z, rms < 1 G in the third panel. For the same snapshot, we also plot a horizontal section of vertical velocity v_z and vertical magnetic field B_z close to the top the convection zone (z = 7.4 Mm) in Fig. 11.

Fig. 10. Vertical sections of the temperature fluctuations T′ (left), vertical velocity vz (middle), and vertical magnetic field Bz (right). To enhance the visibility of the perturbations, each quantity is scaled by its mean rms value at each height z. The sections are taken from the simulation MHD_512 at t = 5000 s and y = 2.34 Mm. Black dotted lines indicate the boundaries of the convection zone.

Fig. 11. Horizontal sections of the vertical velocity vz (top) and vertical magnetic field Bz (bottom) near the top of the convection zone (z = 7.40 Mm). The sections are taken from the MHD_512 simulation at t = 5000 s.

The large-scale structures of the convective flow are the same ones we found in the hydrodynamical simulations: large and slowly ascending convective plumes surrounded by cold and narrow downflows, as we can see in the left and middle panels of Fig. 10 and in the top panel of Fig. 11. Moreover, as for the hydrodynamical simulations, we observe oscillating patterns in temperature and vertical velocity in the stable layers that hint to the presence of waves propagating there.

In the right panel of Fig. 10, we see that the vertical magnetic field is characterized by small-scale filaments with mixed polarity. The typical thickness of these filaments is l ∼ 0.5 Mm, which corresponds to the wavenumber (k ∼ 3 Mm⁻¹) of the peaks in the magnetic energy power spectra seen in Figs. 5 and 6. Strong magnetic fields are found preferentially in intergranular downflows, as it can be seen in the bottom panel of Fig. 11 for the top of the convection zone. We also notice the presence of magnetic fields in thin layers above and below the convection zone in the right panel of Fig. 10. The magnetic field strength quickly decays to zero in these layers, but its presence demonstrates that the plasma is overshooting into the stable regions. Indeed, the frozen-in magnetic field generated in the turbulent convection zone could never be found in the stable layers without overshooting.

3.3. Velocity and magnetic fields properties

Once the simulations attain a quasi-steady magneto-convectional state, we can infer statistical properties of the convective flows and magnetic fields. In this section, all the results for each different simulation represent averages over all snapshots where the saturation phase has been reached. In particular, we consider the last 7 snapshots for the MHD_64 simulation, the last 10 for MHD_128, 7 for MHD_256, and 6 for MHD_512.

From a MLT perspective and assuming isotropical turbulence, we can predict the vertical profile of the one-dimensional velocity dispersion ${\sigma }_{\mathrm{MLT}}^{1\mathrm{D}}$ (see, e.g., Shu 1992) as,

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{MLT}}^{1\mathrm{D}}={\left(\frac{2(\gamma -1){\alpha }_{\mathrm{MLT}}{Q}_{\mathrm{T}}}{\gamma \overline{\rho }}\right)}^{1/3},\end{array}$$(8)

where α_MLT is the mixing-length parameter¹ which we fix to α_MLT = 1.0, Q_T is the convective flux, which we assume to be equal to the stellar energy flux (see Sect. 2.2), and $\overline{\rho }$ is the mean density profile. Assuming quasi-equipartition, we can also predict the vertical profile of the one-dimensional magnetic field strength as,

$$\begin{array}{c}\hfill B_{\mathrm{MLT}}^{1\mathrm{D}}={\alpha }_{\mathrm{B}}\sqrt{4\pi \overline{\rho }}{\sigma }_{\mathrm{MLT}}^{1\mathrm{D}},\end{array}$$(9)

where α_B = B_rms/B_eq is the ratio between the rms magnetic field strength and the equipartition value found in the MHD simulations and shown in Table 2.

Figure 12 shows the vertical profiles of the rms and mean components of the vertical and horizontal velocities, defined as v_z and $v_{\mathrm{h}}=\frac{1}{\sqrt{2}}\sqrt{v_x^2+v_{\mathit{y}}^2}$, respectively. The rms vertical velocity grows as we approach the top of the convection zone, as predicted by MLT since the density profile decreases with z. The maximal amplitude of the rms profile is reached just before the artificial cooling layer. Near the two boundaries and in the bottom stable region, the rms vertical velocity amplitude is very weak, while in the top stable layer we see the imprints of vertical oscillations. On the other hand, the rms horizontal velocity profile peaks near the boundaries of the convection zone. At the top, the uprising convective plumes hit the top stable layer and are cooled down by the artificial cooling. Therefore, the plasma is pushed horizontally toward the intergranular downflows, where it begins its descent. In the deep convection zone, the descending cool plasma encounters the bottom stable stratification and the artificial heating region, thus it gets heated up and it is channeled horizontally into the convective cells. Moreover, we observe large amplitudes of horizontal velocities (in both rms and mean components) in the top stable layer. These are also imprints of waves propagating horizontally in the stable stratification.

Fig. 12. Vertical profiles of the vertical (left) and horizontal (right) components of the velocity field for the four MHD simulations. The rms profiles are shown in continuous lines, while dashed lines represent the mean profiles. Dotted vertical lines denote the boundaries of the convection zone. The MLT prediction for the one-dimensional dispersion velocity given by Eq. (8) is shown in green.

Qualitatively, the velocity profiles are comparable to what Viallet et al. (2013) obtained for the convective envelope of a 5 M_⊙ red giant and Miesch et al. (2008), Hotta et al. (2014, 2015) for the solar convection zone. Moreover, in Fig. 13 we show the vertical profiles of the average three-dimensional Mach number, ℳ = v_rms/c_s, for each one of the MHD simulations. We compare the simulated profiles to a MLT estimate given by ${\mathcal{M}}_{\mathrm{MLT}}={\sigma }_{\mathrm{MLT}}^{3\mathrm{D}}/c_{\mathrm{s}}=\sqrt{3}{\sigma }_{\mathrm{MLT}}^{1\mathrm{D}}/c_{\mathrm{s}}$, where ${\sigma }_{\mathrm{MLT}}^{1\mathrm{D}}$ is given in Eq. (8). We can once more appreciate the extreme subsonic nature of the convective flows, with the maximum mean Mach number being in the order of ℳ ∼ 10⁻² near the top of the convection zone. In addition, we show that our code yields satisfying results also with lower Mach numbers in Appendix A.

Fig. 13. Comparison between Mach number profiles derived from the MHD simulations and from mixing-length theory (MLT) with αMLT = 1.0. Black dotted lines denote the boundaries of the convection zone, while the gray areas represent the artificial heating (left) and cooling (right) regions.

We find that the velocity profiles are in good agreement with MLT theory and that they converge with resolution in the convection zone. In the upper stable layer however, the Mach number profiles from the different simulations do not overlap. The growth of kinetic energy in this layer, seen also in Fig. 12, is due to trapped internal gravity waves. Ideally, these waves would be free to exit the simulation box, but since the top boundary conditions are fixed, the top stable layer acts as a resonant cavity. We will explore outflow boundary conditions with proper characteristic tracing in future work. Internal gravity waves can transport energy at the smallest scales in the simulation (see Sect. 3.4), therefore the better the spatial resolution, the higher the kinetic energy that can be provided to such modes.

Figure 14 is the magnetic analog of Fig. 12. Just as for the velocity, we define the vertical and horizontal magnetic field components as B_z and $B_{\mathrm{h}}=\frac{1}{\sqrt{2}}\sqrt{B_x^2+B_{\mathit{y}}^2}$. The rms magnetic magnetic vertical profiles are reversed with respect to the velocity ones, as the peaks are near the bottom of the convection zone. Similar profiles are found by Rempel (2014), Hotta et al. (2014, 2015) for magnetic fields generated by a small-scale dynamo in the solar convection zone. The magnetic field profiles qualitatively follow the MLT prediction and also converge with resolution, but the difference between the different simulations is more noticeable. In the top and bottom stable regions the magnetic field is essentially zero apart from two thin layers above and beneath the convection zone where, as we have seen in Fig. 10, magnetic field is transported by overshooting convective plasma.

Fig. 14. Vertical profiles of the vertical (left) and horizontal (right) components of the magnetic field for the four MHD simulations. The rms profiles are shown in continuous lines, while dashed lines represent the mean profiles. Dotted vertical lines denote the boundaries of the convection zone. The MLT prediction for the magnetic field strength according to Eq. (9) and with αB = 0.36 is shown in green.

Figure 15 shows the probability density functions (PDFs) of the vertical and horizontal components of the velocity field at two different heights. Near the top of the convection zone (z = 7.0 Mm), the vertical velocity is characterized by an asymmetric PDF, with the probability peak being at low positive values and a long tail toward strong downflows. Such a PDF stems from the granular structure of the convective flow, where the large and slowly uprising plumes dominate the volume but the strongest velocities are found in the intergranular downflows. Closer to the bottom (z = 3.0 Mm), the qualitative shape of the vertical velocity PDF is similar, but the peak is closer to zero and the probability of having large negative velocities is lower. This follows from the decrease of the average rms velocities with depth due to the stellar stratification (see Fig. 12 and Eq. (8)). The horizontal velocity PDFs appear to follow a chi distribution with two degrees of freedom at both heights, because v_x and v_y are Gaussian distributed. Moreover, their peaks are found at the same values of the vertical ones. Once more, the distribution in the deep convection zone is narrower because of stellar stratification.

Fig. 15. Probability density functions (PDF) of the vertical (left column) and horizontal (right column) components of the velocity field in the four MHD simulations after reaching the quasi-steady state. The PDFs are computed over a horizontal section at z = 7.0 Mm (top row) and at z = 3.0 Mm (bottom row), corresponding to top and bottom convection zone, respectively.

The PDFs of the vertical magnetic field are shown in the left panels of Fig. 16. The profiles are symmetric at both heights in the convection zone, which is compatible with the origin of these magnetic fields being due to a turbulent small-scale dynamo with no preferred direction. The profiles are narrower with respect to the velocity ones, showing a higher intermittency of the magnetic fields (Brandenburg et al. 1996). The tails are more extended at z = 3.0 Mm since in the deeper layers of the convection zone the mean magnetic field amplitude is larger (see Fig. 14). Similar conclusions can be inferred from the right panels of Fig. 16, where we show the PDFs of the horizontal magnetic field, B_h.

Fig. 16. Probability density functions (PDF) of the vertical (left column) and horizontal (right column) components of the magnetic field in the four MHD simulations after reaching the quasi-steady state. The PDFs are computed over a section at z = 7.0 Mm (top row) and at z = 3.0 Mm (bottom row), corresponding to top and bottom convection zone, respectively.

Both in Figs. 15 and 16 the PDFs appear to converge with resolution. Moreover, qualitatively similar results have been obtained by Browning et al. (2004), Miesch et al. (2008), Rempel (2014), and Hotta et al. (2014, 2015) in the context of solar convection. This validates our numerical treatment of magneto-convection flows with a very different code.

3.4. Waves

Our fully compressible approach allows us to study the excitation and propagation of acoustic (or pressure) waves, in addition to internal gravity waves. Therefore, in order to characterize the oscillatory phenomena visible in the stable regions above and below the convection zone in both hydrodynamical and MHD simulations (see Figs. 2 and 10), we performed a Fourier analysis on a time series of N_s = 9000 snapshots of a horizontal section of the vertical velocity v_z taken at z = 8.5 Mm, that is in the top stable layer and ∼0.7 Mm above the convection zone. We used the simulations HD_256 and MHD_256 and we started the analysis at t = 5000 s in both cases. The sampling rate is Δt_s = 0.5 s, hence we cover a frequency domain ranging from ω_min = 2π/T = 1.40 × 10⁻³ Hz, where T = N_sΔt_s = 4500 s, to the Nyquist frequency ω_Nyquist = π/Δt_s = 6.28 Hz. The horizontal wavenumber instead, which we define as $k_{\text{h}}^2=k_x^2+k_y^2$, ranges from k_h,  min = 5.71 × 10⁻¹ Mm to k_{h,  Nyquist} = 73.1 Mm given the resolution and the size of the box. Before performing a Fourier transform over the time sequence, we apodized the data by a cosine-bell function to avoid artifacts coming from the nonperiodicity of the signal. We did not need to apply the same procedure for the spatial Fourier transform since the lateral boundaries of the box are periodic.

The results of the Fourier analysis performed over the hydrodynamical simulation HD_256 and MHD simulation MHD_256 are presented as diagnostic diagrams, or k_h − ω diagrams, in Fig. 17. In such diagrams, a single wave form would appear as a dark spot at the relative frequency ω and horizontal wavenumber k_h, while more complex oscillations and mode families can give rise to rays or ridges. In a stratified, compressible medium, the behavior of a generic wave is characterized by the local dispersion relation (see, e.g., Priest 2014),

$$\begin{array}{c}\hfill k_z^2=\frac{({\omega }^2-{\omega }_{\mathrm{ac}}^2)}{c_{\mathrm{s}}^2}-\frac{({\omega }^2-N_{\mathrm{BV}}^2)k_{\mathrm{h}}^2}{{\omega }^2},\end{array}$$(10)

Fig. 17. Wave diagnostic diagrams, also knows as kh − ω diagrams, for the HD_256 (left) and MHD_256 (right) simulations at z = 8.5 Mm. The sampling rate and grid cell size are Δts = 0.5 s and Δx = 42.0 km, respectively. The strength of the signal is color-coded. The dispersion relations for sound waves and surface modes are shown in dashed and dotted lines, respectively. The Brunt–Väisälä frequency is depicted by a dash-dotted line.

where k_z is the vertical wavenumber and c_s is the sound speed. The acoustic cutoff frequency ω_ac (see, e.g., Deubner & Gough 1984) and the Brunt–Väisälä frequency N_BV are defined as,

$$\begin{array}{cc}& {\omega }_{\mathrm{ac}}^2=\frac{c_{\mathrm{s}}^2}{4H_{\rho }^2}(1-2\frac{\mathrm{d}{H}_{\rho }}{\mathrm{d}z}),\hfill \end{array}$$(11)

$$\begin{array}{cc}& N_{\mathrm{BV}}^2=g(\frac{1}{H_{\rho }}-\frac{1}{\gamma H_{\mathrm{p}}}),\hfill \end{array}$$(12)

where H_ρ and H_p are the density and pressure scale heights, respectively. The solutions of setting $k_z^2=0$ in Eq. (10) separate two domains where waves can propagate ($k_z^2>0$) and one region where waves are evanescent ($k_z^2<0$). For k_h > 1/(2H_ρ), the boundaries in the frequency domain between propagating and evanescent regions can be approximated by the Brunt–Väisälä frequency at the lower boundary and by the dispersion relation for sound waves in a homogeneous, compressible medium ω = c_sk_h for the upper boundary. In our case, 1/2(2H_ρ) ∼ 0.5 Mm⁻¹, so this approximation is valid for most of the parameter space displayed in Fig. 17. In the same figure, we show these two relations as dash-dotted and dashed lines, respectively, where c_s and N_BV are computed using our equilibrium profiles. Moreover, we also show the dispersion relation for surface gravity waves (f-modes), $\omega =\sqrt{g{k}_{\mathrm{h}}}$, with a dotted line. For more details on the local dispersion relation of waves in a gravitationally stratified background, the reader can refer to Vigeesh et al. (2017).

The two panels of Fig. 17 appear very similar qualitatively. In both cases, we observe several different modal components. Gravity-modified pressure waves (p-modes) are responsible for the signal in the top left of both panels, that is above the ω = c_sk_h relation. At the top of the diagram we observe a reflection of the signal due to an aliasing phenomenon, also reported by Herwig et al. (2006). Individual ridges due to internal gravity waves (g-modes) can be observed at small horizontal wavenumbers (k_h ≲ 3 Mm⁻¹) and just below the Brunt–Väisälä frequency, while at smaller scales they cannot be resolved individually and give rise to a horizontal band along ω = N_BV. The g-modes are excited by the convective plumes perturbing the stratified plasma in the stable region, while p-modes are mainly produced and propagate within the convective region. The result is a horizontal flow with an oscillatory component (g-modes), that resemble waves at the surface of the see, and the propagation of an over-pressure in the stable layers (p-modes). The turbulent convective noise is responsible for the strong smeared signal in the bottom left of both panels which extends through the evanescent zone. This shows once more that the convective flow is able to overshoot into the stable layers and there perturb vertically the plasma. Similar results have been obtained by Herwig et al. (2006) for the same hydrodynamical setup as used here and by Meakin & Arnett (2007) for an oxygen shell burning 23 M_⊙ stellar model.

To better expose the difference between the two plots of Fig. 17, we show in Fig. 18 the power averaged over the horizontal wavenumber k_h, ${⟨{\widehat{v}}_z⟩}_k$ (left), and over the frequency ω, ${⟨{\widehat{v}}_z⟩}_{\omega }$ (right), for both simulations. The curves are in both panels qualitatively similar, but the hydrodynamical simulation display slightly higher values. The two peaks at low frequency (ω ∼ 0.1 Hz and ω ∼ 1 Hz) correspond to the convective noise signal and to the g-mode ridge visible in Fig. 17. The local crest at high frequency (ω ∼ 5.5 Hz) corresponds to the p-modes signal, while the single peak at low horizontal wavenumber (k_h ∼ 2 Mm⁻¹) is mainly due to the convective noise signal.

Fig. 18. Comparison between the mean power as a function of frequency ω (left) and of horizontal wavenumber kh (right) for the HD_256 and MHD_256 simulations at z = 8.5 Mm. The bottom panels show the relative difference between the two signals in percentage. The green dashed line shows a smoothed interpolation of the data points.

The bottom panels of Fig. 18 show the relative difference between the two curves computed as,

$$\begin{array}{c}\hfill \mathrm{Difference}=\frac{{⟨{\widehat{v}}_z⟩}_i^{\mathrm{MHD}}-{⟨{\widehat{v}}_z⟩}_i^{\mathrm{HD}}}{{⟨{\widehat{v}}_z⟩}_i^{\mathrm{HD}}},\end{array}$$(13)

where i = k_h,  ω, respectively. We see that waves in the MHD_256 simulation have, in average, less power than in HD_256 for all frequencies and horizontal wavenumbers. In total, the relative difference between the two is −11.5%. By using the Brunt–Väisäla frequency and the sound waves dispersion relation to separate the different contributions, we identify a deficit of 13.1% in power related to internal gravity waves and 16.3% to acoustic waves. Moreover, the relative difference varies with ω and k_h. We can identify two dips of ∼ − 20% in the left plot: one at very low and another one at large frequency. The first one (ω ∼ 0.1 Hz) is related to the convective noise signal, which appear to be stronger in the hydrodynamical simulation, while the second one (ω ∼ 5.5 Hz) refers to the p-modes signal. If the first dip can be easily explained with the difference of convective kinetic energy between the two simulations due to the feedback of magnetic fields (see Figs. 8 and 9), the deficit associated with p-modes is not necessarily expected. Although the power associated with pressure waves scales with the eight power of the characteristic turbulent velocity (Lighthill 1952), magneto-acoustic waves can be excited by means of Lorentz forces. Therefore, the loss in power due to a lower turbulent velocity could be compensated by the action of the magnetic fields. The two smallest differences are found in correspondence with the peak of g-modes (ω ∼ 1 Hz) and around ω ∼ 2 Hz, where the MHD run display a small bump associated with a slightly stronger signal in the evanescent region. This signal is mysterious in origin and we speculate could be caused by the magnetic fields in the top convection zone.

A detailed analysis of the causes behind these differences are beyond the scope of this paper. However, there are clear hints indicating the impact of magnetic fields on the generation and propagation of wave modes in the upper stable layer.

4. Summary and conclusions

In this paper we presented a proof of concept for three-dimensional, fully compressible, numerical simulations of magneto-convection in a Cartesian grid with the RAMSES code. In order to deal with low-Mach number and low-amplitude fluctuations typical of convective motions, we adapted the numerical solvers of RAMSES by implementing a well-balanced scheme.

We ran hydrodynamical and MHD simulations of a He shell flash convective envelope model with different resolutions, ranging from N = 64³ to N = 512³. A turbulent convection zone naturally develops from random perturbations in the initial density profile, and the emerging rms velocity field is compatible with MLT predictions. In the Appendix A we show that our code is able to deal with flows characterized by Mach numbers as low as ℳ ∼ 10⁻³. We note that, in principle, our well-balanced scheme allows to model even lower Mach number flows. The only limitation regards the computational cost. Indeed, as we lower the characteristic convective velocity of the simulation, we lower the time-step as well because of the Courant condition. Regarding the MHD simulations, a turbulent small-scale dynamo effectively amplifies a weak seed magnetic field up to ∼18−37% of the equipartition value depending on the resolution. Once saturated, the magnetic fields dominate the dynamics of the convection zone at small scales. In the two stable layers, we saw indications of convective overshoot and we identified propagating pressure and internal gravity waves.

The resulting properties of the convection region and of the magnetic fields converge by increasing the simulation resolution, as we can see from the power spectra analysis (Figs. 6 and 9), from the vertical profiles (Figs. 12 and 14), and from the PDF distributions (Figs. 15 and 16) of velocity and magnetic fields. In particular, we observe an almost statistically equivalent result for the velocity field for the resolutions N = 128³,  256³, and 512³, while the difference is more pronounced regarding the magnetic fields. In both cases, the low-resolution run (N = 64³) presents important departures from the expected results and it is therefore insufficient to capture the magneto-convective dynamics of this setup. Indeed, the low-resolution run resolves the minimum pressure scale heights ($H_{\mathrm{p}}^{\min }\sim 0.7\mathrm{Mm}$) with only ∼4−5 grid-cells, which is not sufficient for an (magneto-)hydrodynamical code. We recall that RAMSES employs a second-order numerical scheme in both hydrodynamical and MHD solvers. Using a higher-order scheme could improve the results even at low-resolutions. Only for the velocity field in the top stable layer we do not observe a convergence with resolution (see Figs. 12 and 13). There, the fixed top boundary condition creates a resonant cavity which traps the waves generated at the interface with the convection zone. A more realistic treatment of the boundary conditions will address this problems in future works.

Qualitatively, we find similar results to several other studies on convection and magneto-convection. Our results differ from the two-dimensional simulations obtained by Herwig et al. (2006) with the same initial model. The reason for this discrepancy is the artificial cooling layer that we introduced at the top of the convection zone to mimic radiative losses. Indeed, the cooling layer causes the intergranular downflows which drive the convective motions in our simulations (see Fig. 2), whereas in Herwig et al. (2006) convection is induced by hot upflow fingers generated in the heated bottom layers. This explains why our results are more similar to surface solar convection simulations (see, e.g., Rempel 2014; Hotta et al. 2014, 2015) than to core convection simulations (see, e.g., Browning et al. 2004; Brun et al. 2005). In any case, our goal was not to reproduce realistic results regarding a particular stellar convection scenario, but to demonstrate the feasibility of simulating three-dimensional compressible magneto-convection with the RAMSES code, and in that respect we consider the proof of concept successful.

The analysis carried out with respect to the propagation of oscillations in the stable region above the convection zone revealed a deficit of power in the magneto-hydrodynamical simulation. In particular, because of our fully compressible approach, we could measure a deficit of 16.3% in power associated with pressure waves. This difference with respect to the purely hydrodynamical case cannot be readily explained by the reduction of convective kinetic energy due to a small-scale dynamo action in the MHD simulation. Indeed, magnetic fields amplified by the same small-scale dynamo can also be sources of magneto-acoustic waves. Therefore, we provide clear indications that magnetic fields play a non trivial and important role in the excitation and propagation of acoustic waves in stellar interiors.

In the future, we aim to expanding our methodology to global models of stellar magneto-convection. In particular, we are interested in the propagation of waves within stars. Our compressible approach, similarly to Horst et al. (2020), allows for the correct treatment of both pressure and internal gravity modes excited by the convective motions. However, the coupling between magnetic fields and these oscillations is of particular interest for asteroseismological observations. For example, the magnetic greenhouse effect could trap wave-modes in the core and therefore suppress dipole oscillations (Fuller et al. 2015; Stello et al. 2016). Magneto-hydrodynamical numerical simulations of fully compressible convection, such as the ones presented in this paper, represent a viable solution to study these phenomena.

The price to pay though is the high computational cost required. Indeed, the MHD_512 simulation required 2.5 × 10⁶ core h with 4608 cores to run for 6300 s physical time, while for the MHD_256 simulation we used 1152 cores and a total of 1.7 × 10⁵ core h to simulate 9000 s physical time. That corresponds approximately to 6 and 9 convective turnovers times, respectively. Horst et al. (2020) estimate a need of 44 × 10⁶ core h to simulate a full 3 M_⊙ stellar model for 700 h physical time (∼17 convective turnover times) on a spherical grid of size 960 (r)×360 (ϑ)×720 (φ) with the compressible Seven-League Hydro (SLH) code (Miczek 2013). Our MHD simulations on a Cartesian grid would increase the computational time cost required, because of the higher number of degrees of freedom and because of the higher resolution needed given the geometry of the grid. However, even taking into account the extra costs, we estimate our numerical simulations to be rather competitive. Moreover, the AMR capabilities of the RAMSES code could help to reduce this gap. Another possibility is to alleviate the stringent time-step constraints imposed by the highly subsonic convective flow using the reduced speed of sound technique (Rempel 2005; Hotta et al. 2015).

An area of improvement concerns energy conservation. As stated in Sect. 2.1, we solve for the conservation of specific entropy instead of total energy, since entropy is the fundamental quantity governing the dynamics of convectively unstable systems. However, this choice implies that energy conservation is not ensured by the numerical scheme. A direct consequence is that the numerical energy flux in the convection zone is larger than the theoretical one by tens of percent. We believe that, in the future, this problem can be addressed by employing higher-order schemes to solve for the total energy conservation and at the same time accurately capturing the dynamics of low amplitude entropy perturbations. Moreover, we plan to improve our simulation setup by avoiding discontinuities in the heating and cooling functions as well as in the equilibrium profiles. These discontinuities are sources of truncation errors that can lead to spurious energy fluxes.

In conclusion, we proved the feasibility of three-dimensional, fully compressible numerical simulations of magneto-convection in the low-Mach number regime with the RAMSES code. This allowed us to study the self-consistent amplification of magnetic fields in stellar convection zones together with their interaction with the generation and propagation of pressure and internal gravity waves in the adjacent stellar cores, radiative zones, or stellar atmospheres.

Acknowledgments

J.R.C.C. acknowledges support by SNF under grant ID 200020_182094. The simulations have been run on the Piz Daint and Eiger machines at the Swiss National Supercomputer Centre (CSCS) under projects uzh5, uzh3, and s1006. We would like to sincerely thank the anonymous referee, O. Steiner, G. Vigeesh, and E. Quataert for very enlightening discussions and comments.

References

Appendix A: Low-Mach number simulations

In this appendix, we demonstrate that our fully compressible approach is able to deal with low-Mach number flows. For this purpose, we lower the volumetric heating and cooling rates of the model presented in Sect. 2.2. Doing so, we reduce the stellar energy flux and luminosity, and consequently also the characteristic convective velocity and typical Mach number.

We adjusted the volumetric heating and cooling rates so that the stellar luminosity reaches L = 3.21 × 10⁴ L_⊙, which is 10³ times lower than the nominal value used in this paper. According to MLT theory and Eq. (8), the characteristic convective velocity is reduced by a factor 10. We ran magneto-hydrodynamical simulations with both nominal and low luminosity setups with resolutions N = 64³,  128³, and 256³. We also reduced the amplitude of the initial random density perturbations δρ to δρ/ρ^eq ∼ 10⁻⁵. We ran the different simulations for 2.5 convective time scales, which corresponds to 2 500 s physical time for the nominal luminosity simulations and 25 000 s for the low luminosity ones.

Figure A.1 shows the time evolution of the mean turbulent kinetic energy for both nominal and low luminosity simulations. We see that the nominal luminosity runs quickly reach a quasi-steady state around E_K ∼ 10¹⁶ erg cm⁻³, just as in Fig. 4. Therefore, the amplitude of the initial perturbations does not affect the properties of convection. The low luminosity simulations instead all stabilize around E_K ∼ 10¹⁴ erg cm⁻³. This value validates the MLT prediction of a reduction of convective velocities by a factor ∼10.

Fig. A.1. Time evolution of the mean turbulent kinetic energy density EK in the convection zone for simulations with nominal (L = 3.21 × 107 L⊙, blue) and low (L = 3.21 × 104 L⊙, purple) stellar luminosities. The convective time scales are assumed to be τconv = 1 000 s for the nominal simulations (see Eq. 7) and τconv = 10 000 s for the low luminosity ones.

In Fig. A.2 we plot the mean Mach number profiles for the six simulations. We used outputs from the last 1.5 τ_conv of each simulation to compute the mean profiles, so that the quasi-steady state is already attained. The nominal luminosity ones are very similar to the ones shown in Fig. 13, while the low ones only reach ℳ ∼ 10⁻³ in the convection zone. In both cases the simulations are in good agreement with the MLT theory and the results converge with resolution in the convection zone.

Fig. A.2. Mach number profiles derived from the nominal and low luminosity simulations. MLT predictions with αMLT = 1.0 are shown for both scenarios. Black dotted lines denote the boundaries of the convection zone, while the gray areas represent the artificial heating (left) and cooling (right) regions.

In conclusion, we proved that the well-balanced version of the RAMSES code presented in this paper is able to simulate low-Mach number (magneto-)convection, down to ℳ ∼ 10⁻³. Moreover, we have shown once more the good agreement between our numerical simulations and MLT predictions.

All Tables

All Figures

	Fig. 1. Initial and equilibrium profiles of density $\stackrel{\sim }{\rho }$, pressure $\stackrel{\sim }{p}$, and specific entropy $\stackrel{\sim }{s}$. Black dotted lines delimit the convection zone, while the gray areas represent the heating (bottom) and cooling (top) layers.
In the text

	Fig. 2. Temperature fluctuations T′ in a vertical section of the HD_512 simulation at t = 100 s (left), t = 400 s (middle), and t = 1000 s (right). The vertical section is taken at y = 3.43 Mm. Black dotted lines indicate the boundaries of the convection zone.
In the text

	Fig. 3. Horizontal sections of the vertical velocity vz (top) and temperature fluctuation T′ (bottom) in the upper layers of the convection zone (z = 7.4 Mm). The snapshot is taken from the HD_512 simulation at t = 1000 s.
In the text

	Fig. 4. Time evolution of the mean turbulent kinetic energy density EK in the convection zone for the four hydrodynamical (HD) simulations.
In the text

	Fig. 5. Time evolution of the turbulent kinetic energy power spectra ${\widehat{E}}_{\mathrm{K}}$ for the four HD simulations. The time evolution is represented by the color grading: light blues correspond to early times, while the darkest shade corresponds to t = 2000 s. The interval between each color-shade is 50 s.
In the text

	Fig. 6. Comparison of the kinetic energy power spectra ${\widehat{E}}_{\mathrm{K}}$ between the four HD simulations in the quasi-steady state regime. The profiles are obtained by averaging the power spectrum profiles over the last 10 snapshots of each simulation, that is between t  =  1500 s and t  =  2000 s.
In the text

	Fig. 7. Time evolution of the mean turbulent kinetic energy density EK (blue) and of the mean magnetic energy density EM (red) in the convection zone for the four magneto-hydrodynamical (MHD) simulations. Dashed lines represent exponential fits obtained with the kinematic growth rates γK listed in Table 2. The top panel is a zoom-in on the first 20 000 s of the turbulent kinetic energy evolution.
In the text

Fig. 8. Time evolution of the power spectra of the turbulent kinetic energy ${\widehat{E}}_{\mathrm{K}}$ (blue) and magnetic energy ${\widehat{E}}_{\mathrm{M}}$ (red) for the four magneto-hydrodynamical (MHD) simulations. The time evolution is represented by the color grading: light colors correspond to early times, while the darkest shade corresponds to the last snapshot of the simulation. The time interval between each shade is 5000 s for MHD_64, 1000 s for MHD_128, 500 s for MHD_256, and 300 s for MHD_512. The average quasi-steady state kinetic energy power spectrum for the hydrodynamical simulation HD_512 is shown in green.

In the text

	Fig. 9. Comparison between the power spectra of the turbulent kinetic energy (blue) and magnetic energy (red), ${\widehat{E}}_{\mathrm{K}}$ and ${\widehat{E}}_{\mathrm{M}}$, for the four magneto-hydrodynamical (MHD) simulations after saturation of the magnetic field. The average quasi-steady state kinetic energy power spectrum for the hydrodynamical simulation HD_512 is shown in green.
In the text

Fig. 10. Vertical sections of the temperature fluctuations T′ (left), vertical velocity vz (middle), and vertical magnetic field Bz (right). To enhance the visibility of the perturbations, each quantity is scaled by its mean rms value at each height z. The sections are taken from the simulation MHD_512 at t = 5000 s and y = 2.34 Mm. Black dotted lines indicate the boundaries of the convection zone.

In the text

	Fig. 11. Horizontal sections of the vertical velocity vz (top) and vertical magnetic field Bz (bottom) near the top of the convection zone (z = 7.40 Mm). The sections are taken from the MHD_512 simulation at t = 5000 s.
In the text

	Fig. 12. Vertical profiles of the vertical (left) and horizontal (right) components of the velocity field for the four MHD simulations. The rms profiles are shown in continuous lines, while dashed lines represent the mean profiles. Dotted vertical lines denote the boundaries of the convection zone. The MLT prediction for the one-dimensional dispersion velocity given by Eq. (8) is shown in green.
In the text

	Fig. 13. Comparison between Mach number profiles derived from the MHD simulations and from mixing-length theory (MLT) with αMLT = 1.0. Black dotted lines denote the boundaries of the convection zone, while the gray areas represent the artificial heating (left) and cooling (right) regions.
In the text

Fig. 14. Vertical profiles of the vertical (left) and horizontal (right) components of the magnetic field for the four MHD simulations. The rms profiles are shown in continuous lines, while dashed lines represent the mean profiles. Dotted vertical lines denote the boundaries of the convection zone. The MLT prediction for the magnetic field strength according to Eq. (9) and with αB = 0.36 is shown in green.

In the text

	Fig. 15. Probability density functions (PDF) of the vertical (left column) and horizontal (right column) components of the velocity field in the four MHD simulations after reaching the quasi-steady state. The PDFs are computed over a horizontal section at z = 7.0 Mm (top row) and at z = 3.0 Mm (bottom row), corresponding to top and bottom convection zone, respectively.
In the text

	Fig. 16. Probability density functions (PDF) of the vertical (left column) and horizontal (right column) components of the magnetic field in the four MHD simulations after reaching the quasi-steady state. The PDFs are computed over a section at z = 7.0 Mm (top row) and at z = 3.0 Mm (bottom row), corresponding to top and bottom convection zone, respectively.
In the text

Fig. 17. Wave diagnostic diagrams, also knows as kh − ω diagrams, for the HD_256 (left) and MHD_256 (right) simulations at z = 8.5 Mm. The sampling rate and grid cell size are Δts = 0.5 s and Δx = 42.0 km, respectively. The strength of the signal is color-coded. The dispersion relations for sound waves and surface modes are shown in dashed and dotted lines, respectively. The Brunt–Väisälä frequency is depicted by a dash-dotted line.

In the text

	Fig. 18. Comparison between the mean power as a function of frequency ω (left) and of horizontal wavenumber kh (right) for the HD_256 and MHD_256 simulations at z = 8.5 Mm. The bottom panels show the relative difference between the two signals in percentage. The green dashed line shows a smoothed interpolation of the data points.
In the text

	Fig. A.1. Time evolution of the mean turbulent kinetic energy density EK in the convection zone for simulations with nominal (L = 3.21 × 107 L⊙, blue) and low (L = 3.21 × 104 L⊙, purple) stellar luminosities. The convective time scales are assumed to be τconv = 1 000 s for the nominal simulations (see Eq. 7) and τconv = 10 000 s for the low luminosity ones.
In the text

	Fig. A.2. Mach number profiles derived from the nominal and low luminosity simulations. MLT predictions with αMLT = 1.0 are shown for both scenarios. Black dotted lines denote the boundaries of the convection zone, while the gray areas represent the artificial heating (left) and cooling (right) regions.
In the text