Issue 
A&A
Volume 664, August 2022



Article Number  A24  
Number of page(s)  15  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202142754  
Published online  03 August 2022 
Toward fully compressible numerical simulations of stellar magnetoconvection with the RAMSES code
^{1}
IRSOL Istituto Ricerche Solari “Aldo e Cele Daccò” Locarno, Università della Svizzera Italiana (USI), Via Patocchi 57 – Prato Pernice, 6605 LocarnoMonti, Switzerland
email: jose.canivete@irsol.usi.ch
^{2}
Center for Theoretical Astrophysics and Cosmology, Institute for Computational Science (ICS), University of Zurich, Winterthurerstrasse 190, 8057 Zürich, Switzerland
^{3}
Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA
Received:
25
November
2021
Accepted:
7
June
2022
Context. Numerical simulations of magnetoconvection have greatly expanded our understanding of stellar interiors and stellar magnetism. Recently, fully compressible hydrodynamical simulations of fullstar models have demonstrated the feasibility of studying the excitation and propagation of pressure and internal gravity waves in stellar interiors, which would allow for a direct comparison with asteroseismological measurements. However, the impact of magnetic fields on such waves has not been taken into account yet in threedimensional simulations.
Aims. We conduct a proof of concept for the realization of threedimensional, fully compressible, magnetohydrodynamical numerical simulations of stellar interiors with the RAMSES code.
Methods. We adapted the RAMSES code to deal with highly subsonic turbulence, typical of stellar convection, by implementing a wellbalanced scheme in the numerical solver. We then ran and analyzed threedimensional hydrodynamical and magnetohydrodynamical simulations with different resolutions of a planeparallel convective envelope on a Cartesian grid.
Results. Both hydrodynamical and magnetohydrodynamical simulations develop a quasisteady, turbulent convection layer from random density perturbations introduced over the initial profiles. The convective flows are characterized by smallamplitude fluctuations around the hydrodynamical equilibrium of the stellar interior, which is preserved over the whole simulation time. Using our compressible wellbalanced scheme, we were able to model flows with Mach numbers as low as ℳ ∼ 10^{−3}, but even lower Mach number flows are possible in principle. In the magnetohydrodynamical runs, we observe an exponential growth of magnetic energy consistent with the action of a smallscale dynamo. The weak seed magnetic fields are amplified to mean strengths of 37% relative to the kinetic equipartition value in the highest resolution simulation. Since we chose a compressible approach, we see imprints of pressure and internal gravity waves propagating in the stable regions above and beneath the convection zone. In the magnetohydrodynamical case, we measured a deficit in acoustic and internal gravity wave power with respect to the purely hydrodynamical counterpart of 16% and 13%, respectively.
Conclusions. The wellbalanced scheme implemented in RAMSES allowed us to accurately simulate the smallamplitude, turbulent fluctuations of stellar (magneto)convection. The qualitative properties of the convective flows, magnetic fields, and excited waves are in agreement with previous studies in the literature. The power spectra, profiles, and probability density functions of the main quantities converge with resolution. Therefore, we consider the proof of concept to be successful. The deficit of acoustic power in the magnetohydrodynamical simulation shows that magnetic fields must be included in the study of pressure waves in stellar interiors. We conclude by discussing future developments.
Key words: stars: interiors / convection / waves / magnetohydrodynamics (MHD)
© J. R. Canivete Cuissa and R. Teyssier 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1. Introduction
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 smallscale dynamo operating in the nearsurface turbulent convection zone of the Sun is possibly at the origin of the smallscale magnetic fields permeating the quiet photosphere (Lites et al. 2014; Rempel 2018), while the cyclic regeneration of the solar largescale magnetic field probably stems from the interplay between deepconvection zone turbulence, differential rotation, and magnetic flux transport (Charbonneau 2013, 2020). Moreover, convective cores in A and Btype stars are likely able to develop a dynamo action (Brun et al. 2005; Augustson et al. 2016) and all lowmass stars appear to host dynamoproduced, smallscale, 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 HerzsprungRussell diagram (Hekker & ChristensenDalsgaard 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öhmVitense 1958). However, largescale (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 timestepping constraints imposed by the lowMach number regime, numerical schemes solving the NavierStokes (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 threedimensional, fully compressible, magnetohydrodynamical 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 magnetohydrodynamics (MHD) in presence of gravity by employing a MUSCLHancock 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 highresolution 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 smallamplitude 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 timeexplicit 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 wellbalanced method.
2.1. Wellbalanced scheme
A wellbalanced 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 shallowwater 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 onedimensional EulerPoisson set of equations,
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,
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,
where the equilibrium states of density and pressure, and , exactly satisfy the hydrostatic equilibrium, Eq. (2), while the velocity field u is made up of perturbations only since . With this separation at hand, we can rewrite the EulerPoisson equations as,
where the hydrostatic equilibrium in stationary regimes, that is where ρ′,p′,u = 0, is explicitly preserved.
Wellbalanced 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 wellbalanced version of the evolution equations. Indeed, we extended the wellbalanced methodology also to the equations of ideal MHD, where we assume the equilibrium state to be characterized by . 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,
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, planeparallel, stellar convective region. We chose to follow the initial setup presented by Herwig et al. (2006), where an intershell of a typical lowmass 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 mainsequence initial mass of 2 M_{⊙} and metallicity Z = 0.01, which undergoes its secondtolast 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^{−2} and a monoatomic ideal gas with adiabatic index γ = 5/3 and equation of state
Fig. 1. Initial and equilibrium profiles of density , pressure , and specific entropy . Black dotted lines delimit the convection zone, while the gray areas represent the heating (bottom) and cooling (top) layers. 
Initial conditions of the He shell flash convective region model.
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 wellbalanced scheme. Hence, to set off convection, we introduced random density perturbations δρ of the order 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 , where and ρ_{0} 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^{22} erg s^{−1} cm^{−2}. 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^{7} 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 magnetohydrodynamical (MHD) numerical simulations of the convective envelope described above. In both cases, we used the wellbalanced 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 magnetohydrodynamical cases, we set a constant, homogeneous, and horizontal (along the x axis) magnetic field of strength B_{0} = 10^{3} 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 smallscale 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^{8} 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 magnetohydrodynamical simulations with resolutions N = 64^{3}, 128^{3}, 256^{3}, and 512^{3}. 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,
where L_{CZ} ∼ 6 Mm is the height of the convection zone and v_{conv} ∼ 1.2 × 10^{6} cm s^{−1} 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 , at times t = 100, 400, and 1000 s, respectively, for the HD_512 simulation. Smallscale thermal convective instabilities developing from the initially perturbed cells can be seen at time t = 100 s. These fingerlike 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 largescale structure of the flow is therefore granularlike, 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 largescale 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 v_{z} (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 wellbalanced 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, , where v_{rms} is the rootmeansquared (rms) turbulent threedimensional 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 quasisteady state is reached, characterized by a quasiconstant mean turbulent kinetic energy density. The runs HD_128, HD_256, and HD_512 reach a mean turbulent kinetic energy density of around 10^{16} erg cm^{−3}, while the low resolution one, HD_64, stabilizes around a slightly lower value. We notice long lived fluctuations in the quasisteady state mean kinetic energy density that are due to the episodic and largescale nature of the flow in the convection zone (Meakin & Arnett 2007).
Fig. 4. Time evolution of the mean turbulent kinetic energy density E_{K} 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 lowresolution simulations (HD_64 and HD_128) show a fast and sharp evolution to a quasisteady state, while a smooth and slower growth characterizes the highresolution ones (HD_256 and HD_512). We can qualitatively explain this difference with the transition between initial instabilities and largescale 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 gridsize level. On the other hand, the kinetic energy dominant structures are the largescale 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, largescale 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, , 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 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 colorshade is 50 s. 
When the quasisteady 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 quasisteady state phase between the different runs. The better the spatial resolution, the larger the range where the Kolmogorov cascade is well reproduced. The quasisteady state energy peak instead is found around k ∼ 1 Mm^{−1} for all resolutions. This value corresponds to the typical size of the largescale uprising convective plumes observed in Fig. 3, that is l = 2π/k ∼ 6 Mm.
Fig. 6. Comparison of the kinetic energy power spectra between the four HD simulations in the quasisteady 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. Smallscale dynamo
In this section we present the magnetohydrodynamical 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^{2}/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 quasisteady 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(γ_{K}t), where γ_{K} is the kinematic growth rate, and then it reaches a quasiconstant value.
Fig. 7. Time evolution of the mean turbulent kinetic energy density E_{K} (blue) and of the mean magnetic energy density E_{M} (red) in the convection zone for the four magnetohydrodynamical (MHD) simulations. Dashed lines represent exponential fits obtained with the kinematic growth rates γ_{K} listed in Table 2. The top panel is a zoomin 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 smallscale 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 lowresolution 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 smallscale dynamo simulations of the quiet Sun magnetism, is very puzzling and in contrast with Kazantsev’s dynamo theory prediction of γ_{K} ∼ Δx^{−2/3}.
Kinematic growth rates γ_{K} and ratio between saturated rms magnetic field strength B_{rms} and equipartition value B_{eq} for the MHD simulations.
The kinematic phase comes to an end when the magnetic field strength approaches the equipartition value, . At this stage, the magnetic fields are strong enough to start backreacting into the plasma dynamics by means of the Lorentz force. The mean magnetic energy eventually saturates and each simulation reaches a magnetoconvective quasisteady state. Typical magnetic field strengths in the convection zone at this stage are in the order of B ∼ 10^{8} 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 powerlaw 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 (blue) and magnetic energy (red) for the four magnetohydrodynamical (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 quasisteady 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^{−1} (l ∼ 2 Mm). The transition to superequipartition of the magnetic fields at large k coincides with a suppression of kinetic power at the same scales, which follows from the Lorentzforce feedback. Consequently, the kinetic power spectra profiles deviate from the Kolmogorov power law for k ≳ 3 Mm^{−1} (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^{−1}, 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 lowresolution 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^{−1}) and the crossover scale (k ∼ 10 Mm^{−1}). The power spectra evolution and steady state configuration for the highresolution simulations are in qualitative agreement with the results obtained by Rempel (2014) and Hotta et al. (2015) for smallscale 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 smallscale dynamo.
Fig. 9. Comparison between the power spectra of the turbulent kinetic energy (blue) and magnetic energy (red), and , for the four magnetohydrodynamical (MHD) simulations after saturation of the magnetic field. The average quasisteady 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 v_{z} (middle), and vertical magnetic field B_{z} (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 v_{z} (top) and vertical magnetic field B_{z} (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 largescale 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 smallscale filaments with mixed polarity. The typical thickness of these filaments is l ∼ 0.5 Mm, which corresponds to the wavenumber (k ∼ 3 Mm^{−1}) 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 frozenin 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 quasisteady magnetoconvectional 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 onedimensional velocity dispersion (see, e.g., Shu 1992) as,
where α_{MLT} is the mixinglength parameter^{1} 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 is the mean density profile. Assuming quasiequipartition, we can also predict the vertical profile of the onedimensional magnetic field strength as,
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 , 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 onedimensional 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 threedimensional Mach number, ℳ = v_{rms}/c_{s}, for each one of the MHD simulations. We compare the simulated profiles to a MLT estimate given by , where 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^{−2} 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 mixinglength 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 . 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 smallscale 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 quasisteady 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 smallscale 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 quasisteady 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 magnetoconvection 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^{−3} 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 , ranges from k_{h, min} = 5.71 × 10^{−1} 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 cosinebell 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),
Fig. 17. Wave diagnostic diagrams, also knows as k_{h} − ω diagrams, for the HD_256 (left) and MHD_256 (right) simulations at z = 8.5 Mm. The sampling rate and grid cell size are Δt_{s} = 0.5 s and Δx = 42.0 km, respectively. The strength of the signal is colorcoded. 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 dashdotted 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,
where H_{ρ} and H_{p} are the density and pressure scale heights, respectively. The solutions of setting in Eq. (10) separate two domains where waves can propagate () and one region where waves are evanescent (). 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_{s}k_{h} for the upper boundary. In our case, 1/2(2H_{ρ}) ∼ 0.5 Mm^{−1}, 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 dashdotted 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 (fmodes), , 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. Gravitymodified pressure waves (pmodes) are responsible for the signal in the top left of both panels, that is above the ω = c_{s}k_{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 (gmodes) can be observed at small horizontal wavenumbers (k_{h} ≲ 3 Mm^{−1}) 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 gmodes are excited by the convective plumes perturbing the stratified plasma in the stable region, while pmodes are mainly produced and propagate within the convective region. The result is a horizontal flow with an oscillatory component (gmodes), that resemble waves at the surface of the see, and the propagation of an overpressure in the stable layers (pmodes). 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}, (left), and over the frequency ω, (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 gmode ridge visible in Fig. 17. The local crest at high frequency (ω ∼ 5.5 Hz) corresponds to the pmodes signal, while the single peak at low horizontal wavenumber (k_{h} ∼ 2 Mm^{−1}) 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 k_{h} (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,
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 pmodes 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 pmodes is not necessarily expected. Although the power associated with pressure waves scales with the eight power of the characteristic turbulent velocity (Lighthill 1952), magnetoacoustic 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 gmodes (ω ∼ 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 threedimensional, fully compressible, numerical simulations of magnetoconvection in a Cartesian grid with the RAMSES code. In order to deal with lowMach number and lowamplitude fluctuations typical of convective motions, we adapted the numerical solvers of RAMSES by implementing a wellbalanced scheme.
We ran hydrodynamical and MHD simulations of a He shell flash convective envelope model with different resolutions, ranging from N = 64^{3} to N = 512^{3}. 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^{−3}. We note that, in principle, our wellbalanced 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 timestep as well because of the Courant condition. Regarding the MHD simulations, a turbulent smallscale 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^{3}, 256^{3}, and 512^{3}, while the difference is more pronounced regarding the magnetic fields. In both cases, the lowresolution run (N = 64^{3}) presents important departures from the expected results and it is therefore insufficient to capture the magnetoconvective dynamics of this setup. Indeed, the lowresolution run resolves the minimum pressure scale heights () with only ∼4−5 gridcells, which is not sufficient for an (magneto)hydrodynamical code. We recall that RAMSES employs a secondorder numerical scheme in both hydrodynamical and MHD solvers. Using a higherorder scheme could improve the results even at lowresolutions. 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 magnetoconvection. Our results differ from the twodimensional 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 threedimensional compressible magnetoconvection 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 magnetohydrodynamical 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 smallscale dynamo action in the MHD simulation. Indeed, magnetic fields amplified by the same smallscale dynamo can also be sources of magnetoacoustic 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 magnetoconvection. 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 wavemodes in the core and therefore suppress dipole oscillations (Fuller et al. 2015; Stello et al. 2016). Magnetohydrodynamical 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^{6} 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^{5} 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^{6} 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 SevenLeague 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 timestep 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 higherorder 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 threedimensional, fully compressible numerical simulations of magnetoconvection in the lowMach number regime with the RAMSES code. This allowed us to study the selfconsistent 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
 Aerts, C. 2021, Rev. Mod. Phys., 93, 015001 [Google Scholar]
 Andrassy, R., Herwig, F., Woodward, P., & Ritter, C. 2020, MNRAS, 491, 972 [NASA ADS] [CrossRef] [Google Scholar]
 Augustson, K. C., Brun, A. S., & Toomre, J. 2016, ApJ, 829, 92 [Google Scholar]
 Basu, S., & Hekker, S. 2020, Front. Astron. Space Sci., 7, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Beck, P. G., Bedding, T. R., Mosser, B., et al. 2011, Science, 332, 205 [Google Scholar]
 Beeck, B., Schüssler, M., Cameron, R. H., & Reiners, A. 2015, A&A, 581, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 BöhmVitense, E. 1958, ZAp, 46, 108 [NASA ADS] [Google Scholar]
 Brandenburg, A., & Subramanian, K. 2005, Phys. Rep., 417, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A., Jennings, R. L., Nordlund, Å., et al. 1996, J. Fluid Mech., 306, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Browning, M. K., Brun, A. S., & Toomre, J. 2004, ApJ, 601, 512 [Google Scholar]
 Browning, M. K., Miesch, M. S., Brun, A. S., & Toomre, J. 2006, ApJ, 648, L157 [Google Scholar]
 Brun, A. S., & Palacios, A. 2009, ApJ, 702, 1078 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073 [Google Scholar]
 Brun, A. S., Browning, M. K., & Toomre, J. 2005, ApJ, 629, 461 [Google Scholar]
 Charbonneau, P. 2013, SaasFee Advanced Course, 39, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Charbonneau, P. 2020, Liv. Rev. Sol. Phys., 17, 4 [Google Scholar]
 Deubner, F.L., & Gough, D. 1984, ARA&A, 22, 593 [NASA ADS] [CrossRef] [Google Scholar]
 Edelmann, P. V. F., Ratnasingam, R. P., Pedersen, M. G., et al. 2019, ApJ, 876, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Edelmann, P. V. F., Horst, L., Berberich, J. P., et al. 2021, A&A, 652, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Freytag, B., Steffen, M., & Dorch, B. 2002, Astron. Nachr., 323, 213 [Google Scholar]
 Freytag, B., Steffen, M., Ludwig, H. G., et al. 2012, J. Comput. Phys., 231, 919 [Google Scholar]
 Freytag, B., Liljegren, S., & Höfner, S. 2017, A&A, 600, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fuller, J., Cantiello, M., Stello, D., Garcia, R. A., & Bildsten, L. 2015, Science, 350, 423 [Google Scholar]
 Ghizaru, M., Charbonneau, P., & Smolarkiewicz, P. K. 2010, ApJ, 715, L133 [NASA ADS] [CrossRef] [Google Scholar]
 Gomes, P., & Lopes, I. 2020, MNRAS, 496, 620 [NASA ADS] [CrossRef] [Google Scholar]
 Greenberg, J. M., & Leroux, A. Y. 1996, SIAM J. Numer. Anal., 33, 1 [Google Scholar]
 Hekker, S., & ChristensenDalsgaard, J. 2017, A&ARv, 25, 1 [Google Scholar]
 Herwig, F. 2005, ARA&A, 43, 435 [NASA ADS] [CrossRef] [Google Scholar]
 Herwig, F., & Austin, S. M. 2004, ApJ, 613, L73 [NASA ADS] [CrossRef] [Google Scholar]
 Herwig, F., Freytag, B., Hueckstaedt, R. M., & Timmes, F. X. 2006, ApJ, 642, 1057 [NASA ADS] [CrossRef] [Google Scholar]
 Horst, L., Edelmann, P. V. F., Andrássy, R., et al. 2020, A&A, 641, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Horst, L., Hirschi, R., Edelmann, P. V. F., Andrássy, R., & Röpke, F. K. 2021, A&A, 653, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hotta, H., Rempel, M., & Yokoyama, T. 2014, ApJ, 786, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Hotta, H., Rempel, M., & Yokoyama, T. 2015, ApJ, 803, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Hotta, H., Iijima, H., & Kusano, K. 2019, Sci. Adv., 5, 2307 [Google Scholar]
 Houdek, G., & Dupret, M.A. 2015, Liv. Rev. Sol. Phys., 12, 8 [Google Scholar]
 Käppeli, R., & Mishra, S. 2016, A&A, 587, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J. 2021, A&A, 651, A66 [Google Scholar]
 Käpylä, P. J., Käpylä, M. J., Olspert, N., Warnecke, J., & Brandenburg, A. 2017, A&A, 599, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kupka, F., & Muthsam, H. J. 2017, Liv. Rev. Comput. Astrophys., 3, 1 [Google Scholar]
 Langer, N. 2014, in Magnetic Fields throughout Stellar Evolution, eds. P. Petit, M. Jardine, & H. C. Spruit, 302, 1 [NASA ADS] [Google Scholar]
 Lighthill, M. J. 1952, Proc. R. Soc. London Ser. A, 211, 564 [NASA ADS] [CrossRef] [Google Scholar]
 Lites, B. W., Centeno, R., & McIntosh, S. W. 2014, PASJ, 66, S4 [Google Scholar]
 Meakin, C. A., & Arnett, D. 2007, ApJ, 667, 448 [NASA ADS] [CrossRef] [Google Scholar]
 Miczek, F. 2013, Ph.D. Thesis, Technische Universität München, Germany [Google Scholar]
 Miesch, M. S., Brun, A. S., DeRosa, M. L., & Toomre, J. 2008, ApJ, 673, 557 [NASA ADS] [CrossRef] [Google Scholar]
 Noelle, S., Xing, Y., & Shu, C.W. 2007, J. Comput. Phys., 226, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Nordlund, Å., Stein, R. F., & Asplund, M. 2009, Liv. Rev. Sol. Phys., 6, 2 [Google Scholar]
 Pietarila Graham, J., Cameron, R., & Schüssler, M. 2010, ApJ, 714, 1606 [Google Scholar]
 Priest, E. 2014, Magnetohydrodynamics of the Sun (Cambridge: Cambridge University Press) [Google Scholar]
 Rempel, M. 2005, ApJ, 622, 1320 [NASA ADS] [CrossRef] [Google Scholar]
 Rempel, M. 2014, ApJ, 789, 132 [Google Scholar]
 Rempel, M. 2018, ApJ, 859, 161 [Google Scholar]
 Riva, F., & Steiner, O. 2022, A&A, 660, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rogers, T. M., Lin, D. N. C., McElwaine, J. N., & Lau, H. H. B. 2013, ApJ, 772, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Salaris, M., & Cassisi, S. 2017, R. Soc. Open Sci., 4, 170192 [Google Scholar]
 Salhab, R. G., Steiner, O., Berdyugina, S. V., et al. 2018, A&A, 614, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shu, F. H. 1992, Physics of Astrophysics II. Gas Dynamics, 1st edn. (Mill Valley: University Science Books) [Google Scholar]
 Stein, R. F. 2012, Liv. Rev. Sol. Phys., 9, 4 [Google Scholar]
 Stein, R. F., & Nordlund, Å. 1998, ApJ, 499, 914 [Google Scholar]
 Stello, D., Cantiello, M., Fuller, J., et al. 2016, Nature, 529, 364 [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [CrossRef] [EDP Sciences] [Google Scholar]
 Teyssier, R., Fromang, S., & Dormy, E. 2006, J. Comput. Phys., 218, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Tremblay, P. E., Fontaine, G., Freytag, B., et al. 2015, ApJ, 812, 19 [Google Scholar]
 Veiga, M. H., VelascoRomero, D. A., Abgrall, R., & Teyssier, R. 2019, Commun. Comput. Phys., 26, 1 [CrossRef] [Google Scholar]
 Viallet, M., Meakin, C., Arnett, D., & Mocák, M. 2013, ApJ, 769, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Vigeesh, G., Jackiewicz, J., & Steiner, O. 2017, ApJ, 835, 148 [NASA ADS] [CrossRef] [Google Scholar]
 Woodward, P. R., Herwig, F., & Lin, P.H. 2015, ApJ, 798, 49 [Google Scholar]
Appendix A: LowMach number simulations
In this appendix, we demonstrate that our fully compressible approach is able to deal with lowMach 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^{4} L_{⊙}, which is 10^{3} 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 magnetohydrodynamical simulations with both nominal and low luminosity setups with resolutions N = 64^{3}, 128^{3}, and 256^{3}. We also reduced the amplitude of the initial random density perturbations δρ to δρ/ρ^{eq} ∼ 10^{−5}. 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 quasisteady state around E_{K} ∼ 10^{16} erg cm^{−3}, 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^{14} erg cm^{−3}. 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 E_{K} in the convection zone for simulations with nominal (L = 3.21 × 10^{7} L_{⊙}, blue) and low (L = 3.21 × 10^{4} 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 quasisteady 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^{−3} 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 wellbalanced version of the RAMSES code presented in this paper is able to simulate lowMach number (magneto)convection, down to ℳ ∼ 10^{−3}. Moreover, we have shown once more the good agreement between our numerical simulations and MLT predictions.
All Tables
Kinematic growth rates γ_{K} and ratio between saturated rms magnetic field strength B_{rms} and equipartition value B_{eq} for the MHD simulations.
All Figures
Fig. 1. Initial and equilibrium profiles of density , pressure , and specific entropy . 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 v_{z} (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 E_{K} in the convection zone for the four hydrodynamical (HD) simulations. 

In the text 
Fig. 5. Time evolution of the turbulent kinetic energy power spectra 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 colorshade is 50 s. 

In the text 
Fig. 6. Comparison of the kinetic energy power spectra between the four HD simulations in the quasisteady 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 E_{K} (blue) and of the mean magnetic energy density E_{M} (red) in the convection zone for the four magnetohydrodynamical (MHD) simulations. Dashed lines represent exponential fits obtained with the kinematic growth rates γ_{K} listed in Table 2. The top panel is a zoomin 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 (blue) and magnetic energy (red) for the four magnetohydrodynamical (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 quasisteady 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), and , for the four magnetohydrodynamical (MHD) simulations after saturation of the magnetic field. The average quasisteady 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 v_{z} (middle), and vertical magnetic field B_{z} (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 v_{z} (top) and vertical magnetic field B_{z} (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 onedimensional 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 mixinglength 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 quasisteady 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 quasisteady 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 k_{h} − ω diagrams, for the HD_256 (left) and MHD_256 (right) simulations at z = 8.5 Mm. The sampling rate and grid cell size are Δt_{s} = 0.5 s and Δx = 42.0 km, respectively. The strength of the signal is colorcoded. 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 dashdotted line. 

In the text 
Fig. 18. Comparison between the mean power as a function of frequency ω (left) and of horizontal wavenumber k_{h} (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 E_{K} in the convection zone for simulations with nominal (L = 3.21 × 10^{7} L_{⊙}, blue) and low (L = 3.21 × 10^{4} 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.