© The Authors 2023

1. Introduction

Convection plays a key role in the evolution of stars. In the deep, optically thick layers, convective flows can efficiently transport energy and angular momentum outward, so they determine both the thermal structure and the rotational profile of stars (see, e.g., Maeder 2009; Kippenhahn et al. 2013). Furthermore, because the characteristic spatial scales of convection are much larger than the mean free path in the stellar plasma, stellar convection zones are highly turbulent environments, with Reynolds numbers that can be as high as 10¹⁴ (Jermyn et al. 2022). Turbulent flows quickly mix chemical elements over the relatively short convective turnover timescale, thus profoundly affecting the nuclear energy generation in burning layers of stars and their evolution.

Despite the huge imprint of convection on stars, most one-dimensional (1D) stellar evolution codes still rely on simplified parametrizations of the convective energy transport, such as the popular mixing-length theory (MLT; Prandtl 1925; Böhm-Vitense 1958). On the one hand, these parameterized theories allow 1D models to simulate the evolution of stars over thermal and nuclear timescales, which is still unfeasible in multi-D. On the other hand, the parameters that enter these prescriptions cannot be derived from first principles and need to be calibrated. Usually, their value is tuned so that 1D models can reproduce the global properties of our Sun (see, e.g., Richard et al. 1996), but the universality of this approach has been heavily questioned in the literature (Trampedach et al. 2014; Magic et al. 2015; Joyce & Chaboyer 2018; Sonoi et al. 2019). Moreover, local theories of convection such as the MLT assume that convective mixing stops at the position of the formal convective boundary. More realistically, convective plumes approach the convective boundary with nonzero velocities and give rise to hydrodynamic processes that can entrain some excess mass and entropy from the neighboring stable layer into the convection zone. Evidence of extra mixing occurring at stellar convective boundaries has been provided by a number of observations, including eclipsing binaries (Valle et al. 2016; Claret & Torres 2016), old open clusters (Aparicio et al. 1990), or asteroseismology (Bossini et al. 2015; Aerts 2021). The entrainment of fresh fuel into a burning layer can prolong its lifetime, enlarge convective cores in upper main sequence stars, and determine the structure of supernova progenitors in more massive stars (Müller 2020, and references therein). In stellar evolution codes, mixing at convective boundaries is crudely modeled by means of additional parametrizations, usually in the form of diffusive over-mixing or convective penetration (Anders & Pedersen 2023, and references therein). The uncertainties arising from the usage of such simplistic models limit the predictive power of stellar evolution calculations and have far-reaching consequences for supernova explosions, the formation of stellar populations, and galactic chemical evolution. Although several nonlocal theories of convection have been presented in the literature (see, e.g., Xiong 1978; Kuhfuss 1986; Canuto 1997, 2011; Li & Yang 2007; Garaud et al. 2010), they have not been extensively used in 1D stellar calculations so far.

To overcome the limitations of stellar-evolution models, several research groups have started focusing their efforts in the past two decades on multi-D hydrodynamic modeling of turbulent convection and mixing at convective boundaries in different classes of stellar objects, including core-convective main sequence stars (Gilet et al. 2013; Horst et al. 2020; Higl et al. 2021; Baraffe et al. 2023; Herwig et al. 2023; Andrassy et al. 2023), envelope-convective stars (Pratt et al. 2016; Hotta 2017; Käpylä 2019; Blouin et al. 2023), and more massive stars during late burning stages (Meakin & Arnett 2007; Jones et al. 2017; Cristini et al. 2017; Andrassy et al. 2020; Rizzuti et al. 2023). In this approach, nonlinear hydrodynamic processes are captured self-consistently, which allows parameterized theories of convection and convective boundary mixing to be tested and calibrated for different stellar masses and evolutionary stages.

One more layer of complexity to the problem of stellar convection is, however, represented by the possible presence of magnetic fields, which have been observed both in low- and high-mass stars (Brun & Browning 2017; Keszthelyi 2023, and references therein). The coupling between turbulent fluid motions and magnetic fields can give rise to small-scale dynamos (SSDs; Meneguzzi et al. 1981; Brandenburg & Subramanian 2005; Schekochihin et al. 2007), which amplify magnetic fields on scales smaller than the forcing scale of turbulence. As observed in numerous simulations of solar convection (Vögler & Schüssler 2007; Pietarila Graham et al. 2010; Rempel 2014; Thaler & Spruit 2015; Hotta et al. 2016; Hotta 2017), SSDs can drastically change the morphology of the convective flows, reduce their speed, and alter the dynamics of the overshoot region at the bottom of the solar convection zone as compared to the purely hydrodynamic case.

Other than the Sun, effects of magnetohydrodynamic (MHD) processes on the properties of convective flows have also been investigated in cool (Browning 2008; Käpylä 2021; Bhatia et al. 2022) and upper-main-sequence stars (Brun et al. 2005; Featherstone et al. 2009; Augustson et al. 2016), but very few MHD simulations of late burning stages of massive stars have been run to date (Varma & Müller 2021, 2023; Canivete Cuissa & Teyssier 2022). In contrast to the main sequence phase, these late evolutionary stages are characterized by vigorous convective shells that can entrain a substantial amount of mass on relatively short timescales, possibly giving rise to shell mergers (Ritter et al. 2018; Mocák et al. 2018; Yadav et al. 2020; Andrassy et al. 2020). If an efficient SSD action takes place inside these shells, it could reduce the mass entrainment rate at the convective boundaries and possibly delay or even inhibit the occurrence of the merger events. Numerical simulations of the dynamical amplification of magnetic fields in these layers then become essential for determining the stratification of the supernova progenitor and the fate of the star. Acquiring more insight into dynamo mechanisms in late burning shells of massive stars is also particularly important because magnetic fields can act as seeds to magneto-rotationally powered supernovae (Müller & Varma 2020).

In this paper, we investigated the effects of a small-scale turbulent dynamo acting in a late stellar convective shell using our fully compressible, MHD, SEVEN-LEAGUE HYDRO (SLH) code. In particular, we used an idealized setup whose stratification is close to that of an oxygen-burning shell in a massive star (Andrassy et al. 2022). In this study, we did not intend to perform realistic simulations of such an evolutionary stage. Instead, we checked the efficiency of small-scale dynamo action in generating dynamically relevant magnetic fields on the typical timescales set by convective motions in the oxygen shell and quantified their impact on the evolution of the convection and on the mixing at the convective boundary.

The paper is structured as follows: in Sect. 2, we give a brief description of the numerical methods used to run the simulations needed for this study. The details of the initial stratification are provided in Sect. 3. In Sect. 4, we present the numerical results, including the evolution of the small-scale turbulent dynamo and its effects on the boundary mixing. Finally, in Sect. 5, we draw conclusions and summarize the main results.

2. Methods

2.1. Equations solved

We described the physical problem by means of the fully compressible equations of ideal MHD with time-independent gravity,

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·(\rho \mathit{V})=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& \frac{\partial (\rho \mathit{V})}{\partial t}+\mathrm{\nabla }·[\rho \mathit{V}\otimes \mathit{V}+(p+p_B)\mathbb{I}-\mathit{B}\otimes \mathit{B}]=\rho \mathit{g},\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\partial (\rho e_{\mathrm{tot}})}{\partial t}+\mathrm{\nabla }·[(\rho e_{\mathrm{tot}}+p+p_B)\mathit{V}-\mathit{B}(\mathit{B}·\mathit{V})]=0,\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\partial \mathit{B}}{\partial t}+\mathrm{\nabla }·(\mathit{V}\otimes \mathit{B}-\mathit{B}\otimes \mathit{V})=\mathbf{0},\hfill \end{array}$$(4)

$$\begin{array}{cc}& \frac{\partial (\rho \psi )}{\partial t}+\mathrm{\nabla }·(\rho \psi \mathit{V})=0,\hfill \end{array}$$(5)

where ρ denotes the density, V = (V_x, V_y, V_z) the velocity vector, B = (B_x, B_y, B_z) the magnetic field¹, p_B = |B|²/2 the magnetic pressure, g = (g_x, g_y, g_z) the gravitational acceleration, e_tot = e_int + |V|²/2 + |B|²/(2ρ)+ϕ the total energy per unit mass, e_int the internal energy per unit mass, ϕ the gravitational potential, ψ the mass fraction of a passive tracer, and 𝕀 the unit tensor. The system in Eqs. (1)–(5) is closed by an Equation of State (EoS) for the gas pressure p,

$$\begin{array}{c}\hfill p=p(\rho ,e_{\mathrm{int}}).\end{array}$$(6)

In our simulations, we assumed an ideal gas law,

$$\begin{array}{c}\hfill p(\rho ,e_{\mathrm{int}})=(\gamma -1)\rho e_{\mathrm{int}},\end{array}$$(7)

where γ = 5/3 is the adiabatic index.

We stress that the absence of the viscous and resistive dissipation terms in Eqs. (2)–(4) does not mean that the simulated flows are inviscid and nonresistive. In fact, the numerical methods that we used to solve the equations of ideal MHD (see Sect. 2.2) must add a certain amount of numerical dissipation into the system in order to achieve numerical stability.

2.2. Numerical methods

We solved Eqs. (1)–(5) with the SLH code, which is suited for simulating low-Mach-number (magneto-)convection and excitation of internal gravity waves (IGWs) in the deep interiors of stars (Miczek et al. 2015; Edelmann et al. 2017, 2021; Horst et al. 2020, 2021; Andrassy et al. 2022, 2023; Leidi et al. 2022). SLH makes use of a second-order finite-volume discretization and upwinding techniques that require an approximate solution to the Riemann problem at each cell interface. Here, the pair of Riemann states was reconstructed using the Piecewise-Parabolic-Method (PPM) of Colella & Woodward (1984). Upwind, hyperbolic fluxes were computed at cell interfaces with the low-dissipation version of the HLLD solver (LHLLD) of Minoshima & Miyoshi (2021). LHLLD modifies the stabilizing pressure-diffusion term in the original HLLD solver of Miyoshi & Kusano (2005) to ensure that the magnitude of numerical dissipation (relative to the physical central flux) is independent of the Mach number of the flow, ℳ = |V|/c, where c = (γp/ρ)^1/2 is the sound speed. This correction dramatically reduces the excessive amount of numerical diffusion introduced by shock-capturing methods in simulations of subsonic flows (Miczek et al. 2015; Leidi et al. 2022).

To suppress the development of spurious flows due to grid discretization errors in strongly stratified setups, we used a well-balancing technique (the Deviation method, Berberich et al. 2021; Edelmann et al. 2021). In this method, the vector of conserved quantities, U = (ρ,  ρV,  ρe_tot,  B,  ρψ), is split into a time-independent hydrostatic component, $\stackrel{\sim }{\mathit{U}}$, and a fully nonlinear perturbation, δU. Equations (1)–(5) are then solved by enforcing $\partial \stackrel{\sim }{\mathit{U}}/\partial t=\mathbf{0}$, which is achieved in practice by subtracting the hydrostatic fluxes and source terms from the spatial residuals. Such a measure is necessary because conventional finite volume methods discretize hyperbolic fluxes and gravitational source terms at different locations on the computational grid, so hydrostatic solutions cannot be maintained for long times. If ignored, these discretization errors can dramatically affect the evolution of buoyancy-driven flows and produce grossly inaccurate numerical solutions (Edelmann et al. 2021).

To keep the strength of magnetic monopoles under control, we used a staggered constrained transport (CT) method (Evans & Hawley 1988). Different from the finite volume discretization, here the surface-averaged magnetic field components are evolved at cell interfaces by performing the line integral of the electromotive force along the cell edges. Thanks to this operation, the update on the cell-volume average of ∇ ⋅ B vanishes to machine precision. In SLH, the upwind electromotive force at the cell edges is computed according to the CT-Contact scheme of Gardiner & Stone (2005).

Finally, both cell-centered and staggered quantities were evolved in time with a semi-discrete scheme based on the method of lines. The resulting system of ordinary differential equations was solved using the explicit strong stability preserving (SSP) RK2 method of Shu & Osher (1988). Further details regarding the implementation of the fully unsplit MHD solver in SLH can be found in Leidi et al. (2022).

3. Setup

We used the setup first described in Andrassy et al. (2022), who performed a comparison of five different hydrodynamic codes (SLH, PPMSTAR, MUSIC, FLASH, and PROMPI) on a problem involving turbulent convection, convective boundary mixing, and the excitation of IGWs in an overlying stable layer. The thermodynamic conditions of this test setup resemble those found during oxygen shell burning in a star with an initial mass of 25 M_⊙ (Jones et al. 2017). However, Andrassy et al. (2022) adopted a few simplifications to make the study easily reproducible by other research groups. In particular, the geometry of the shell was plane-parallel, the EoS was that of an ideal gas, neutrino cooling was not included, and detailed nuclear burning was replaced by a time-independent heat source term, whose amplitude was set such that convective motions were driven with root-mean-square velocities characteristic of late evolutionary stages in massive stars (ℳ_rms ≈ 0.04).

We mapped the initial hydrostatic stratification (see Fig. 1) on a 3D Cartesian grid with spatial domain (x, y, z)∈[−L_ref, L_ref]×[L_ref, 3 L_ref]×[−L_ref, L_ref], where L_ref = 4 × 10⁸ cm. We used periodic boundaries in the horizontal x- and z-direction. At the top and bottom boundaries of the domain, instead, we adopted impermeable, stress-free boundary conditions for the velocity field,

$$\begin{array}{c}\hfill \frac{\partial V_x}{\partial y}=\frac{\partial V_z}{\partial y}=V_y=0,\end{array}$$(8)

Fig. 1. Profiles of density, pressure, pseudo-entropy (p/ργ), gravity, mass fraction of a passive tracer (ψ), and heat source term ($\dot{q}$) as a function of the vertical coordinate y at t = 0 s. Here, ρref = 1.82 × 106 g cm−3, pref = 4.64 × 1023 dyne cm−2, gref = 6.37 × 108 cm s−2, ${\dot{q}}_{\mathrm{ref}}=1.76\times {10}^{20}$ erg cm−3 s−1, and Lref = 4 × 108 cm.

we forced the magnetic field to be purely horizontal,

$$\begin{array}{c}\hfill \frac{\partial B_x}{\partial y}=\frac{\partial B_z}{\partial y}=B_y=0,\end{array}$$(9)

and for the scalar quantities we assumed

$$\begin{array}{c}\hfill \frac{\partial \rho }{\partial y}=\frac{\partial p}{\partial y}=\frac{\partial \psi }{\partial y}=0.\end{array}$$(10)

The gravitational acceleration, assumed to be time-independent, points downward in the y-direction and goes to zero at the vertical boundaries according to Eq. (1) of Andrassy et al. (2022). In that work, such a choice for the gravitational acceleration was made to allow the hydrostatic density and pressure profiles to become constant at the domain boundaries, making the problem consistent with the conditions in Eq. (10). Although unrealistic, turning off the gravity at the boundaries does not appreciably alter the stratification of the oxygen shell, which is mostly affected by the aforementioned simplifications, as can be seen in Fig. 1 of Andrassy et al. (2022). For consistency with their model, we decided not to modify the profile of the gravitational acceleration here.

The stratification is isentropic up to approximately y = 2 L_ref, and it smoothly turns subadiabatic in the upper half of the domain. Overall, the grid covers 4.35 pressure scale heights in the vertical direction. To be able to track the time evolution of the mass entrained into the convection zone, we filled the convectively stable layer with a passive tracer at t = 0 s, whose abundance progressively drops to zero across the convective boundary. Further details regarding the initial stratification and the heat source can be found in Andrassy et al. (2022).

To start the dynamo action, we planted an initially horizontal magnetic field into the grid, B_x = 10⁵ G. The strength of the seed field was chosen such that the Lorentz force exerted on the fluid at early times was weak enough to not affect the development of convection.

We judged the numerical convergence of the results obtained in this work by running simulations on grids with 128³, 256³, and 512³ cells. To compute meaningful time-averaged quantities and avoid introducing temporal correlations caused by the turbulent nature of the convective flows, all test cases were run until t_max = 25τ_conv, where τ_conv = 63.36 s is the mean convective turnover timescale in the purely hydrodynamic case, defined according to Eq. (16) of Andrassy et al. (2022). By t = t_max, the growing convective layer is still sufficiently far away from the upper boundary of the spatial domain that the imposed boundary conditions do not appreciably alter the dynamics of the mixing region. Thus, we decided not to extend the simulations beyond t_max. Finally, in order to capture possible differences between the MHD and the purely hydrodynamic case, we ran an additional set of simulations without magnetic fields.

4. Results

4.1. Onset of convection and kinematic stage of the dynamo

To break the initial symmetry, we added a small-amplitude perturbation to the hydrostatic density stratification according to Eq. (6) of Andrassy et al. (2022). The energy injected by the heat source at the base of the box leads to the development of buoyant parcels of hot fluid that rise in the adiabatic layer (see Fig. 2). As soon as these flows cross the boundary of the subadiabatic layer, the buoyant acceleration changes sign (so it points downward in the y-direction) and forces the rising plumes to turn around. The large-scale buoyant fluid elements that are driven by the energy source quickly develop shear instabilities that cascade down to smaller scales, and turbulent convection fully develops by t ≈ τ_conv.

Fig. 2. Development of convection in the MHD simulation of the idealized oxygen shell run on a 5123 grid. The panels show fluctuations in pseudo-entropy (A = p/ργ) in the z = 0 plane at different times, as indicated by the insets. The entropy generated by the heat source at the base of the box (see Fig. 1) is mixed throughout the initially adiabatic layer by turbulent convection. This process slowly increases the entropy content of the convection zone in time. The broad stripe of negative entropy fluctuation visible in the upper half of the domain at early times is due the thermal expansion of the convective layer. The turbulent flows also excite IGWs at the upper convective boundary (lower center panel) which then propagate in the subadiabatic layer.

As shown in Fig. 3, the mean magnetic energy density inside the convective shell²,

$$\begin{array}{c}\hfill {\stackrel{\sim }{E}}_B=\frac{1}{2}{⟨|\mathit{B}|}^2{⟩}_{\mathrm{conv}},\end{array}$$(11)

Fig. 3. Time evolution of the mean magnetic energy density inside the convection zone for the indicated grid resolutions.

increases exponentially in time. The growth rate of the instability is higher on finer grids, which indicates that the amplification process mostly occurs on intermediate or small spatial scales. It could be an SSD, where the magnetic field is randomly stretched at scales smaller than the forcing scale of turbulence. However, two other processes could contribute to the amplification of small-scale magnetic fields in this setup: turbulent induction, which is the stretching of large-scale magnetic fields by a turbulent, small-scale velocity component (Schekochihin et al. 2007), and turbulent cascade of magnetic energy toward smaller scales (Pietarila Graham et al. 2010). In fact, the large-scale fields that are needed to excite the latter two processes, not only are characteristic of large-scale dynamos (Brandenburg 2009; Charbonneau 2013), but they can also be supported by the large-scale velocity structures typical of turbulent convection (Käpylä et al. 2018). Moreover, the imposed boundary conditions (see Sect. 3) allow the integrated horizontal magnetic flux, and consequently the mean horizontal magnetic field, to be preserved in time inside the spatial domain. Therefore, the mean horizontal field takes the value of the chosen seed field, B_x = 10⁵ G, and represents a persistent large scale magnetic field component that could, in principle, contribute to the amplification of magnetic energy via turbulent induction.

To get a better understanding of the underlying mechanisms that amplify the magnetic energy in these simulations, we computed transfer functions T_XYZ(k) in the Fourier space between the kinetic (K) and magnetic (B) energy reservoirs inside the convective layer, following the approach of Pietarila Graham et al. (2010). In particular, T_XYZ(k) represents the energy received (or lost in case of T_XYZ(k) < 0) per unit time and per unit wavenumber at scale k of energy type Y from all scales of energy type X via process Z. The transfer of magnetic energy to the kth component of kinetic energy is determined by the net work done on the fluid by the magnetic tension force,

$$\begin{array}{cc}\hfill T_{\mathrm{BKT}}(\mathit{k})=& \frac{1}{2}\widehat{\mathit{V}}(\mathit{k})·{[\widehat{\mathit{B}·\mathrm{\nabla }\mathit{B}}]}^{\ast }(\mathit{k})\hfill \\ \hfill & +\frac{1}{2}(\widehat{\rho \mathit{V}})(\mathit{k})·[\widehat{\frac{1}{\rho }\mathit{B}·\mathrm{\nabla }\mathit{B}}{]}^{\ast }(\mathit{k})+\mathrm{c}.\mathrm{c}.,\hfill \end{array}$$(12)

and the magnetic pressure force,

$$\begin{array}{cc}\hfill T_{\mathrm{BKP}}(\mathit{k})=& -\frac{1}{4}\widehat{\mathit{V}}(\mathit{k})·[\widehat{{\mathrm{\nabla }|\mathit{B}|}^2}{]}^{\ast }(\mathit{k})\hfill \\ \hfill & -\frac{1}{4}(\widehat{\rho \mathit{V}})(\mathit{k})·[\widehat{\frac{1}{\rho }\mathrm{\nabla }{|\mathit{B}|}^2}{]}^{\ast }(\mathit{k})+\mathrm{c}.\mathrm{c}.,\hfill \end{array}$$(13)

where is the complex conjugate, c.c. is the complex conjugate of the whole expression on the right-hand side, and ^ represents the Fourier projection³. Magnetic energy on scale k is produced or removed via stretching of the magnetic field lines,

$$\begin{array}{c}\hfill T_{\mathrm{KBS}}(\mathit{k})=\widehat{\mathit{B}}(\mathit{k})·{[\widehat{\mathit{B}·\mathrm{\nabla }\mathit{V}}]}^{\ast }(\mathit{k})+\mathrm{c}.\mathrm{c}.,\end{array}$$(14)

and through compression and advection of the magnetic field,

$$\begin{array}{cc}\hfill T_{\mathrm{KBCA}}(\mathit{k})=& -\widehat{\mathit{B}}(\mathit{k})·{[\widehat{\mathit{B}\mathrm{\nabla }·\mathit{V}}]}^{\ast }(\mathit{k})\hfill \\ \hfill & -\widehat{\mathit{B}}(\mathit{k})·{[\widehat{\mathit{V}·\mathrm{\nabla }\mathit{B}}]}^{\ast }(\mathit{k})+\mathrm{c}.\mathrm{c}.\hfill \end{array}$$(15)

Because here we solved the fully compressible MHD equations, T_KBCA(k) includes both the transport of energy within the magnetic energy reservoir and the generation of magnetic energy through fluid compression. These two processes, however, cannot be decoupled (Rempel 2014), which can be seen by expanding the advective flux of magnetic energy as

$$\begin{array}{c}\hfill -\mathrm{\nabla }·(\mathit{V}\frac{{|\mathit{B}|}^2}{2})(\mathit{k})=-\widehat{\mathit{B}}(\mathit{k})·\widehat{[(\mathit{V}·\mathrm{\nabla })\mathit{B}+\frac{\mathit{B}}{2}\mathrm{\nabla }·\mathit{V}{]}^{\ast }}(\mathit{k})+\mathrm{c}.\mathrm{c}.\end{array}$$(16)

To simplify the calculations, instead of computing 3D Fourier projections in Eqs. (12)–(15), we averaged transfer functions T_XYZ(k_h, y_j) obtained at each horizontal plane y_j inside the convection zone,

$$\begin{array}{c}\hfill T_{\mathit{XYZ}}(k_{\mathrm{h}})={⟨T_{\mathit{XYZ}}(k_{\mathrm{h}},y_j)⟩}_{y_j\in (L_{\mathrm{ref}},2L_{\mathrm{ref}})}.\end{array}$$(17)

Here, k_h is the horizontal wavenumber $k_{\mathrm{h}}=\sqrt{k_x^2+k_z^2}$, where

$$\begin{array}{cc}\hfill k_x& =\{\begin{array}{cc}m,\hfill & 0\le m\le \lfloor \frac{N_x-1}{2}\rfloor ,\hfill \\ -N_x+m,\hfill & \lfloor \frac{N_x-1}{2}\rfloor <m<N_x,\hfill \end{array}\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill k_z& =\{\begin{array}{cc}n,\hfill & 0\le n\le \lfloor \frac{N_z-1}{2}\rfloor ,\hfill \\ -N_z+n,\hfill & \lfloor \frac{N_z-1}{2}\rfloor <n<N_z,\hfill \end{array}\hfill \end{array}$$(19)

⌊.⌋ is the floor function, and N_x and N_z are the number of cells in the x- and z- direction, respectively.

Figure 4 shows results from the transfer analysis performed on the grid with 512³ cells. Stretching of the magnetic field lines contributes most of the magnetic energy generation at spatial wavenumbers close to k_h = 50. In these simulations, the typical velocities in the convection zone are considerably subsonic (ℳ_rms ≈ 0.04), so fluid compression due to the ram pressure of the convective flows (p_ram ∼ ℳ²) has a negligible contribution to the generation of magnetic energy (see $T_{\mathrm{KBC}}(k_{\mathrm{h}})=-\widehat{\mathit{B}}(k_{\mathrm{h}})·{[\widehat{\mathit{B}\mathrm{\nabla }·\mathit{V}}]}^{\ast }(k_{\mathrm{h}})$ + c.c. in Fig. 4). Therefore, T_KBCA(k_h) measures mostly the advective transport of magnetic energy to scale k_h from all scales of the magnetic field, similar to the case of incompressible MHD,

$$\begin{array}{c}\hfill T_{\mathrm{KBCA}}(k_{\mathrm{h}})\approx -\widehat{\mathit{B}}(k_{\mathrm{h}})·{[\widehat{\mathit{V}·\mathrm{\nabla }\mathit{B}}]}^{\ast }(k_{\mathrm{h}})+\mathrm{c}.\mathrm{c}.\end{array}$$(20)

Fig. 4. Transfer functions extracted from the convection zone on the 5123 grid and averaged over the kinematic stage of the dynamo. The panel on the left shows the transfer rates from the kinetic (K) to the magnetic (B) energy reservoir, whereas the panel on the right shows magnetic-to-kinetic energy transfer rates. The time averaging was performed such that, at each time t, all the transfer curves were rescaled by the maximum value of TKBS(kh). The resulting curves were then averaged over the time interval t ∈ (1.5τconv, 3τconv). Dashed lines represent negative transfer rates, while solid lines are used for positive rates.

We observe that the magnetic cascade mainly removes magnetic energy from large scales, where T_KBCA < 0, and redistributes it at scales with k_h > 75, where T_KBCA > 0. This process dominates the generation of magnetic energy over stretching only at k_h > 130, which corresponds to a spatial scale of 3.9 times the width of the grid cells, Δx. Work done by fluid motions against the magnetic tension force (−T_BKT) most efficiently transforms kinetic energy at k_h ≈ 30 into magnetic energy. Work done by the magnetic pressure force on the fluid is negligible everywhere except on very large scales. These results allow us to find the range of wavenumbers where the magnetic field is most efficiently stretched by fluid motions. As pointed out by Pietarila Graham et al. (2010), the spatial wavenumber at which magnetic field is generated (q), the one at which it is stretched (k), and the one at which the flow works against magnetic tension (p) form a triadic relation,

$$\begin{array}{c}\hfill \mathit{k}=\mathit{q}-\mathit{p}.\end{array}$$(21)

By considering the most extreme cases in which q and p have the same or the opposite orientation, we estimate that the magnetic field lines are most efficiently stretched at 20 ≲ k_h ≲ 80. As we show in Sect. 4.2, this interval lies at the bottom of the inertial range of the turbulent kinetic energy spectrum. Thus, the amplification of magnetic energy is mostly caused by the action of a small-scale turbulent dynamo, with a minor contribution from the turbulent cascade close to the grid scale.

Further evidence of small-scale turbulent dynamo action can be provided by checking the scaling of the growth rate of the magnetic energy, $\gamma =\partial \ln {\stackrel{\sim }{E}}_B/\partial t$, with the grid resolution. Because in this work we used the Implicit Large Eddy Simulation (ILES) method, the magnetic Prandtl number (Pr_m = ν/η) is likely to be close to or larger than unity (Vögler & Schüssler 2007; Rempel 2014). In this regime of Prandtl numbers, an SSD can only be started if the fluid Reynolds number, Re = V_rmsL_ref/ν, is larger than a critical threshold. The growth rate of the magnetic energy in an unstable SSD should then scale as Re^1/2 (Kazantsev 1968; Schekochihin et al. 2004). Although the effective value of the kinematic viscosity (ν) and resistivity (η) coefficients are determined by the underlying numerical methods used to solve the MHD equations, in the ILES approach the fluid Reynolds number should depend on the spatial resolution as Δx^−4/3 (Cristini et al. 2017), which leads to γ ∝ Re^1/2 ∝ Δx^−2/3. Our study indicates that the growth rate γ follows the predicted theoretical scaling (see Fig. 5). These results can only be confirmed by using an explicit kinematic viscosity coefficient so that Re can be measured directly, which, however, is beyond the scope of this work.

Fig. 5. Growth rate of the mean magnetic energy inside the convection zone (averaged over the kinematic phase of the dynamo) as a function of the grid spacing, Δx. For this analysis, we also simulated the kinematic phase of the dynamo on grids with 3843 and 6403 cells (Δx/Lref = 5.2 × 10−3, and Δx/Lref = 3.1 × 10−3, respectively). Error bars represent three standard deviations, computed over the time series. The expected theoretical scaling for SSD amplification (γ ∝ Re1/2 ∝ Δx−2/3) is represented by the black dashed line.

At early times, the magnetic energy is still subdominant with respect to the kinetic energy content of the flow. Therefore, the Lorentz force does not affect the evolution of the convection, and we do not observe any systematic difference in the velocity field between the MHD and the hydrodynamic simulations. This is the kinematic stage of the dynamo, which lasts for several convective turnover timescales, depending on the resolution of the grid.

4.2. Nonlinear phase of the dynamo

The amplification of the magnetic field due to the action of the small-scale turbulent dynamo proceeds until the Lorentz force becomes strong enough to start having a feedback effect on the flow. Such a change in the evolution of the dynamo happens when the magnetic energy approaches equipartition with the kinetic energy content of the eddies on the small scales of turbulence. Strong, small-scale magnetic fields inhibit the development of shear instabilities that feed the turbulent cascade and drive the dynamo amplification. The stretching of the magnetic field lines happens now on larger scales, where turbulence has not been quenched. Work done against the magnetic tension force by the turbulent convective flows sustains the magnetic field against numerical (resistive) dissipation, and the dynamo reaches saturation. By t/τ_conv ≈ 15, all of the MHD simulations presented here have entered this phase. This stage of the dynamo, however, does not represent a statistical steady state solution of the simulated setup. In fact, the continuous injection of entropy into the system by the heat source and the mixing processes that take place at the convective boundary (see Sect. 4.3) both contribute to the entrainment of material from the overlying stable layer. Therefore, the size, mass, and entropy content of the convective layer keep increasing over time.

The magnetic-to-kinetic energy ratio, shown in Fig. 6, reaches saturation with a mean value of ≈0.22 on the finest grid, with sporadic, intermittent episodes in which it reaches values as high as 0.33. During the nonlinear phase of the dynamo, the mean kinetic energy density inside the convective shell in the MHD simulations,

$$\begin{array}{c}\hfill {\stackrel{\sim }{E}}_{\mathrm{K}}=\frac{1}{2}{⟨\rho |\mathit{V}|}^2{⟩}_{\mathrm{conv}},\end{array}$$(22)

Fig. 6. Time evolution of the mean magnetic-to-kinetic energy ratio inside the convective shell.

is on average 25% lower than that in the hydrodynamic case on the 512³ grid (see Fig. 7). We obtained this result by first computing the time average { ⋅ } of ${\stackrel{\sim }{E}}_{\mathrm{K}}$ over t ∈ (15τ_conv, 25τ_conv), and then by calculating

$$\begin{array}{c}\hfill ϵ=\{{\stackrel{\sim }{E}}_{\mathrm{K},\mathrm{MHD}}\}/\{{\stackrel{\sim }{E}}_{\mathrm{K},\mathrm{HYDRO}}\}-1.\end{array}$$(23)

Fig. 7. Time evolution of the mean kinetic energy density inside the convection zone, ${\stackrel{\sim }{E}}_{\mathrm{K}}$, for the purely hydrodynamic (vermilion) and MHD (light blue) simulations. The time evolution of the mean magnetic energy density in the MHD simulations, ${\stackrel{\sim }{E}}_B$, is also shown.

However, the large temporal fluctuations that characterize ${\stackrel{\sim }{E}}_{\mathrm{K}}$ made it necessary to provide an error estimate on ϵ in order to prove the statistical significance of this result. The error on ϵ is a combination of the statistical uncertainties on $\{{\stackrel{\sim }{E}}_{\mathrm{K},\mathrm{HYDRO}}\}$ and $\{{\stackrel{\sim }{E}}_{\mathrm{K},\mathrm{MHD}}\}$, which we computed as follows. First we obtained the standard deviation σ_K of ${\stackrel{\sim }{E}}_{\mathrm{K}}$ over the selected time series for both the hydrodynamic and MHD setups. Second, we estimated the statistical uncertainty on the mean quantity $\{{\stackrel{\sim }{E}}_{\mathrm{K}}\}$ by taking into account possible temporal correlations introduced by the turbulent nature of the convective flows. According to Fig. 4 of Andrassy et al. (2022), the autocorrelation function of the convective velocity drops to zero after a time shift Δt ≈ τ_conv, which suggests that there is approximately one independent realization of the convective flows per convective turnover. Thus, we approximated the uncertainty associated to $\{{\stackrel{\sim }{E}}_{\mathrm{K}}\}$, ${\stackrel{\sim }{\sigma }}_{\mathrm{K}}$, as

$$\begin{array}{c}\hfill {\stackrel{\sim }{\sigma }}_{\mathrm{K}}=\frac{{\sigma }_{\mathrm{K}}}{\sqrt{N_{\mathrm{to}}}},\end{array}$$(24)

where N_to = 10 is the number of convective turnovers (see also Table 1). Finally, we computed the variance of ϵ,

$$\begin{array}{c}\hfill {\sigma }_{ϵ}^2=(\frac{\partial ϵ}{\partial \{{\stackrel{\sim }{E}}_{\mathrm{K},\mathrm{HYDRO}}\}}{)}^2{\stackrel{\sim }{\sigma }}_{\mathrm{K},\mathrm{HYDRO}}^2+(\frac{\partial ϵ}{\partial \{{\stackrel{\sim }{E}}_{\mathrm{K},\mathrm{MHD}}\}}{)}^2{\stackrel{\sim }{\sigma }}_{\mathrm{K},\mathrm{MHD}}^2.\end{array}$$(25)

We find σ_ϵ ≈ 3%, making a mean relative deviation of 25% statistically significant (i.e., different from zero) by more than 8σ_ϵ.

In the ILES approach, increasing the grid resolution reduces the amount of numerical resistivity introduced into the system, making the turbulent dynamo progressively more efficient. Consequently, the mean magnetic energy density in our simulations increases by a factor of two from the 128³ to the 512³ grid, where the typical strength of the magnetic field is ≈5 × 10⁹ G (see Table 2). On the other hand, ${\stackrel{\sim }{E}}_{\mathrm{K}}$ does not seem to show a significant resolution dependence when averaged over the saturated phase of the dynamo, t ∈ (15τ_conv, 25τ_conv), as can be noted from the values provided in Table 1. Averaging over a wider time window could potentially reveal a statistically significant trend in $\{{\stackrel{\sim }{E}}_{\mathrm{K}}\}$ with increasing spatial resolution. However, this is not possible with our current setup, in which the convective boundary approaches the upper domain boundary at late times, potentially altering the dynamics of the mixing region and producing unreliable results (see also Sect. 3). We also note that a fixed time-averaging window, in principle, samples different evolutionary times of the dynamo depending on grid resolution. In fact, the oxygen shell does not have a statistical steady state solution, and the time at which the dynamo enters its nonlinear regime is resolution dependent. However, as visible in Fig. 7, neither ${\stackrel{\sim }{E}}_{\mathrm{K}}$ nor ${\stackrel{\sim }{E}}_B$ show clear, long-term trends in the saturated phase of the dynamo. This result suggests that the secular evolution of the mean stratification can be neglected for this analysis since it does not seem to affect basic properties of the flows and of the dynamo.

The suppression of small-scale shear instabilities in the convective flows caused by the generated strong magnetic fields can be noted in Fig. 8, where we compare snapshots of the Mach number taken from an MHD and a purely hydrodynamic simulation in the nonlinear phase of the dynamo. In contrast to the hydrodynamic case, where turbulence is essentially isotropic on spatial scales smaller than L_ref, the velocity field in the MHD simulation is characterized by the presence of anisotropic, thread-like structures that extend over a large part of the convective shell.

Fig. 8. Snapshots of the Mach number on the 5123 grid at t/τconv = 22. The panels on the left show results from the purely hydrodynamic simulation, while the panels on the right show the MHD simulation. The upper row shows vertical cuts taken in the z = 0 plane, while the lower row shows a horizontal cut through the y = 1.5 Lref plane. The white-dotted line indicates the initial position of the convective boundary.

A comparison of horizontally averaged velocity profiles (see Fig. 9) shows that horizontal velocities in the MHD case are reduced by as much as 30% in the convective layer, as compared to the simulations without magnetic fields. Vertical velocities, instead, are diminished on average only by 10%. As visible in Fig. 10, the partial suppression of the horizontal mixing in the MHD runs increases the magnitude of the root-mean-square entropy fluctuation inside the convection zone with respect to the hydrodynamic simulations. A larger contrast in the thermal content between up- and downflows in turn increases the efficiency of the convective energy transport. For this reason, despite the mild suppression of vertical velocities in the convection zone caused by the action of the dynamo, we do not observe any significant difference in the vertical enthalpy fluxes,

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{H}}={⟨(e_{\mathrm{int}}+p)V_y⟩}_{x,z},\end{array}$$(26)

Fig. 9. Vertical profiles of the horizontal ($V_{\mathrm{h}}=\sqrt{V_x^2+V_z^2}$) and vertical (Vy) velocity components averaged over t ∈ (15τconv, 25τconv). Here, light blue is used to indicate the quantities extracted from the MHD simulations, whereas vermilion is used for the hydrodynamic simulations. An estimate of the statistical uncertainty associated to the averaged profiles can be inferred by looking at the typical dispersion among the hydrodynamic simulations, which represent a set of three independent (and numerically converged) realizations of the oxygen shell.

between the MHD and the hydrodynamic setups (see Fig. 11). Only in the overshoot layer, where ℱ_H < 0, the MHD simulations are characterized by a smaller unsigned enthalpy flux than their hydrodynamic counterpart. This result may be due to the combined effects of reduced flow speeds and entropy fluctuations at the upper convective boundary in the MHD case, as visible in Figs. 9 and 10, respectively.

Fig. 10. Vertical profiles of the root-mean-square relative pseudo-entropy fluctuation in the MHD (light blue) and purely hydrodynamic simulations (vermilion), averaged over t ∈ (15τconv, 25τconv). The upper convective boundary is, on average, at the position of the peaks visible at y = (2.3 − 2.4) Lref.

Fig. 11. Vertical profiles of the vertical enthalpy flux (ℱH) averaged over t ∈ (15τconv, 25τconv). Here, light blue is used to indicate the quantities extracted from the MHD simulations, whereas vermilion is used for the hydrodynamic simulations.

In Fig. 12, we show the kinetic and magnetic energy spectra computed in the y = 1.5 L_ref plane, averaged over the saturated phase of the dynamo. The kinetic energy spectra resulting from the hydrodynamic simulations converge to the Kolmogorov scaling ($k_{\mathrm{h}}^{-5/3}$, Kolmogorov 1941) in the inertial range. The scale at which the power spectrum deviates from the Kolmogorov law due to the action of numerical dissipation becomes smaller as the resolution is progressively increased, as expected in the ILES approach. The kinetic energy spectra in the MHD simulations, instead, deviate from the hydrodynamic curves already at wavenumbers of k_h ≈ 10, where the dynamo is most efficient at converting kinetic energy into magnetic energy (i.e., it achieves maximum $|T_{\mathrm{BKT}}|/{\widehat{E}}_{\mathrm{K}}$). The observed drop of kinetic energy in the MHD case corresponds to an increase of magnetic power in the inertial range, and the sum of the two energy contributions approximately resembles the kinetic energy spectrum in the hydrodynamic simulations. The wavenumber at which the magnetic-to-kinetic energy ratio becomes greater than unity decreases on finer grids⁴. On the grid with 512³ cells, the break-even point is at k_h = 10, which corresponds to approximately half of the pressure scale height at y = 1.5 L_ref. Unlike the case of forced, isotropic MHD turbulence (Haugen et al. 2004; Schekochihin et al. 2004; Iskakov et al. 2007; Brandenburg 2011), our simulations of stratified convection also show a substantial amount (up to 30%) of magnetic power stored at wavenumbers k_h < 10. This field component is generated by coherent structures in the form of large-scale up- and downflows that stretch the magnetic field lines over a large fraction of the size of the convection zone. The presence of a large-scale field component can be observed in Fig. 13, where we show vertical and horizontal cuts in B_y for all our three grids. A small-scale component with mixed polarity also becomes more noticeable on progressively finer grids with reduced numerical dissipation.

Fig. 12. Magnetic (light blue) and kinetic (black-dotted for the hydrodynamic simulations, and vermilion for the MHD simulations) energy spectra computed in the y = 1.5 Lref plane and averaged over the time interval t ∈ (15τconv, 25τconv). Each panel shows the curves extracted from a different grid. The black-dashed lines indicate the Kolmogorov scaling law, $k_{\mathrm{h}}^{-5/3}$.

Fig. 13. Distribution of By on the z = 0 plane (upper row) and the y = 1.5 Lref plane (lower row) at t/τconv = 22 for the indicated grid resolutions.

Both the horizontal and the vertical component of the magnetic field become stronger as the grid is refined (see Fig. 14). The magnetic field smoothly turns horizontal across the upper convective boundary, where the convective flows overturn due to the negative buoyant acceleration. At the bottom of the shell, the magnetic field is forced to be horizontal in order to retain its solenoidal property given the imposed boundary conditions. We note that the reflecting boundary forces the convective flows to abruptly change direction over a few computational cells, which artificially enhances the stretching and the compression of the magnetic field lines. This process, however, only affects the generation of magnetic energy in a narrow region close to the bottom boundary of the convective shell. In simulations of SSD action in the solar convection zone, Hotta (2017) finds that including part of the underlying stable layer does not appreciably alter the generation of the magnetic field as compared to simulations with closed boundary conditions. This author shows that an imposed steep, positive entropy gradient across the solar overshoot region prevents convective plumes from penetrating into the stable layer deeper than a small fraction of the pressure scale height, at least on the characteristic timescales set by convection. Thus, in those simulations, the bottom boundary of the solar convective region acts like a reflecting wall for the turbulent flows and the magnetic fields. Oxygen-burning shells of massive stars are also characterized by steep, stabilizing entropy gradients at their bottom boundary (Meakin & Arnett 2007; Jones et al. 2017; Varma & Müller 2021). In the model of Jones et al. (2017, which this setup is based on), the square of the Brunt–Väisälä frequency at the silicon-oxygen boundary is several times larger than that at the upper boundary of the oxygen shell (see Fig. 4 of Jones et al. 2017). Because this quantity is directly related to the buoyancy jump, entrainment of material from the underlying stable layer into the convection zone can easily be neglected over the timescales simulated here. All these considerations give us confidence that the chosen boundary conditions are well suited for the simulations presented in this study.

Fig. 14. Vertical profiles of the horizontal ($B_{\mathrm{h}}=\sqrt{B_x^2+B_z^2}$) and vertical (By) magnetic field, averaged over the time interval t ∈ (15τconv, 25τconv).

The strength of the magnetic field is not uniform across the convective shell. As previously discussed, the magnetic energy approaches equipartition with the kinetic energy content of the turbulent eddies in the inertial range. Because the mean flow speed does not vary considerably across the convective shell, the spatial dependence of the magnetic field strength is mostly set by the density stratification. Indeed, the vertical profile of the root-mean-square magnetic field rescaled by its local equipartition value ($B_{\mathrm{eq}}=\sqrt{\rho }{V}_{\mathrm{rms}}$) shows much less dependence on y as compared with the results shown in Fig. 15. The height-dependence of the dynamo action can also be seen in Fig. 16, where we show horizontal averages of the Lorentz work,

$$\begin{array}{c}\hfill W_{\mathrm{L}}={⟨\mathit{V}·[(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}]⟩}_{x,z}.\end{array}$$(27)

Fig. 15. Vertical profiles of the root-mean-square magnetic field divided by the equipartition field, averaged over the time interval t ∈ (15τconv, 25τconv).

Fig. 16. Vertical profiles of the buoyancy (light blue) and Lorentz work (vermilion) averaged over t ∈ (15τconv, 25τconv) for the indicated grid resolutions.

As a reference, we also show the buoyancy work,

$$\begin{array}{c}\hfill W_{\mathrm{b}}={⟨V_y{g}_y\delta \rho ⟩}_{x,z},\end{array}$$(28)

where δρ is the density fluctuation. Buoyancy generates kinetic energy in the whole convective shell except in the overshoot layer, where it is responsible for the deceleration of the convective flows. On all grids, the magnitude of the Lorentz work is maximum at the bottom of the convective shell and progressively drops to zero toward the upper convective boundary. W_L is negative throughout the whole convective layer, meaning that, on average, kinetic energy is everywhere converted into magnetic energy. Moreover, profiles of the Lorentz work approach convergence on the finest grids. This result confirms that most of the conversion of kinetic energy into magnetic energy happens on relatively large spatial scales, which are well resolved even with moderate grid resolutions.

We note that the dynamo does not operate in the subadiabatic layer, although a seed field is present there as well. The turbulent structures created by the nonlinear breaking of IGWs (visible in Fig. 8), which is one of the mechanisms that can excite an SSD in stable stratifications (Skoutnev et al. 2021), are not efficient enough to build a significant magnetic field in these simulations. The magnetic field in the stable layer reaches saturation with average strengths of only two to ten times that of the initial seed field.

4.3. Impact of magnetic fields on the growth of the convective shell

The turbulent convective flows generated in this setup give rise to a rich variety of hydrodynamic processes at the upper convective boundary, including shear instabilities, breaking of surface gravity waves, and convective overshooting. These processes contribute to the entrainment of high-entropy material from the overlying stable layer into the convection zone, which causes the convective shell to grow in time (see also Fig. 8). We computed the mass entrained per unit surface area inside the convection zone by using horizontal averages of the density and the passive tracer ψ as in Andrassy et al. (2022),

$$\begin{array}{c}\hfill M_{\mathrm{e}}(t)={\displaystyle {\int }_{L_{\mathrm{ref}}}^{y_{\mathrm{cb}}(t)}}\overline{\rho }(y,t)\overline{\psi }(y,t)\mathrm{d}y.\end{array}$$(29)

At each time, we assumed that the vertical coordinate of the upper convective boundary, y_cb, was the position of the steepest gradient in $\overline{\psi }$. In Fig. 17 we show the time evolution of the entrained mass for all the simulations run in this study. In the hydrodynamic case, numerical convergence is reached already on the lowest-resolved grid (with 128³ cells) within a maximum relative statistical uncertainty of 5%, which is consistent with the results obtained by Andrassy et al. (2022). Instead, the mass entrained in the MHD runs slightly decreases with increasing the grid resolution. A significant deviation between the MHD and the hydrodynamic simulations is visible only after 15τ_conv, when the dynamo has fully entered its nonlinear phase. By the time the simulations have finished, the best resolved MHD setup has entrained 12% less mass than the hydrodynamic runs, and the mass entrainment rate Ṁ_e has been reduced by ≈20%.

Fig. 17. Time evolution of the mass entrained from the stable layer into the convection zone, rescaled by the total mass contained in the stable layer at t = 0 s. Light blue is used for the MHD simulations, whereas vermilion is used for the purely hydrodynamic runs.

Because MHD effects do not dramatically reduce the mass entrainment rate at the upper convective boundary, finding the mechanisms responsible for the observed discrepancy in M_e between the MHD and the hydrodynamic results is challenging. One possible explanation is that convective flows in the MHD simulations have some of their kinetic energy converted into magnetic energy by the action of the SSD (see the discussion in Sect. 4.2), which reduces the amount of energy available to overcome the buoyancy of the entrained high-entropy material as compared to the hydrodynamic case (Spruit 2015). Furthermore, strong horizontal magnetic fields (see Fig. 14) could considerably reduce the growth rate of the shear instabilities that take place at the convective boundary, which are in part responsible for the mixing. Magnetic fields aligned with the shear flows, in fact, have a stabilizing effect against the growth of Kelvin–Helmholtz instabilities, especially if the Alfvénic Mach number ${\mathcal{M}}_{\mathrm{Alf}}=\sqrt{\rho }|\mathit{V}|/|\mathit{B}|$ is close to unity (Chandrasekhar 1961; Frank et al. 1996). As shown in Fig. 15, the mean magnetic field at the upper convective boundary reaches values as high as 60% of the equipartition field (ℳ_Alf ≳ 1.5), so short-wavelength Kelvin–Helmholtz instabilities are likely to be partly suppressed.

As pointed out by a number of authors (Meakin & Arnett 2007; Andrassy et al. 2020, 2023; Horst et al. 2020), some degree of mixing at stellar convective boundaries can be induced by nonadiabatic effects. In this setup, the mean entropy inside the convection zone increases by the action of the heat source, and eventually overcomes the entropy level of a narrow subadiabatic layer right above the upper convective boundary. This layer becomes negatively buoyant, so it sinks and gets mixed into the convection zone. This process enlarges the size of the convective shell over time as long as the source of entropy generation is active. In the work of Andrassy et al. (2022), it was estimated that by the end of the simulations, ≈60% of Ṁ_e in this setup was due to the heating. This process is expected to operate regardless of the properties of the convective flows, so MHD processes would only be able to affect the remaining 40% of the entrainment rate. In the absence of heating-induced mixing, magnetic fields would then suppress as much as 50% of the mass entrainment.

5. Summary and discussion

We have run 3D simulations of turbulent convection, dynamo amplification and convective boundary mixing in an idealized oxygen-burning shell of a 25 M_⊙ star. In particular, we have searched for possible MHD effects on the boundary mixing and the properties of the convective flows by performing a comparison between an MHD and a purely hydrodynamic setup. The numerical results have been carefully analyzed by means of a convergence study, in which the grid resolution was progressively increased from 128³ to 512³ cells.

Random stretching of the magnetic field lines due to the turbulent motions in the convective shell excites small-scale dynamo (SSD) action on all of the considered grids. The dynamo instability amplifies the seed field by ≈4 orders of magnitude in a few convective turnover timescales. The kinematic phase of the dynamo ends when the magnetic field becomes strong enough to affect the evolution of the flows on the small scales of turbulence. During the saturated stage, the work done by fluid motions against magnetic tension forces sustains the magnetic field against numerical resistive dissipation. The saturated mean magnetic-to-kinetic energy ratio reaches values in the 20 − 30% range. The magnetic field strength in the oxygen shell moderately increases with the grid resolution, and it has characteristic values of ∼10¹⁰ G in the 512³ simulation. Such strong fields partly suppress the small-scale isotropic features in the velocity field typical of turbulent convection in hydrodynamic simulations. The resulting flows present anisotropic, thread-like structures that extend over a large fraction of the convective shell. The magnetic fields generated during the oxygen burning stage can be further amplified if parts of the oxygen shell end up collapsing onto a neutron star. By assuming a simple flux-freezing model, we estimate that the magnetic field strength at the surface of the neutron star would be on the order of 10¹⁵ − 10¹⁶ G, which is in agreement with values inferred from observations of magnetars (Kouveliotou et al. 1998; Woods & Thompson 2006; Olausen & Kaspi 2014).

Vertical and horizontal fluid velocities in the bulk of the convective layer in the MHD simulations are reduced, as compared to the hydrodynamic runs, on average by 10% and 20%, respectively. The fact that the dynamo does not have the same impact on the different fluid velocity components could be related to the transport of thermal energy inside the convection zone. In fact, in order for convection to transport the excess heat deposited by the energy source outward with partly suppressed vertical velocities, the thermal content of the buoyant flows must be enhanced. Because the heat source is unchanged in the MHD simulations, this can only be achieved by reducing the horizontal mixing of entropy between the up- and the downflows. Indeed, we observe that root-mean-square entropy fluctuations are systematically enhanced in the MHD simulations, which in turn increases the thermal contrast between the convective plumes. Consequently, we do not observe a significant impact of magnetic fields on the vertical enthalpy fluxes.

Power spectra computed in the bulk of the convective shell reveal that 30% of the total magnetic energy is stored at spatial wavenumbers k_h < 10. Such a large-scale field component is generated by large-scale convective cells, which can efficiently stretch the magnetic field lines on length scales comparable to the size of the convection zone. The kinetic energy spectra in the MHD simulations deviate from the Kolmogorov law ($k_{\mathrm{h}}^{-5/3}$) in the inertial range, where the efficiency of the dynamo is maximum. On the finest grid, the magnetic energy becomes greater than the kinetic energy at a spatial wavenumber of k_h ≈ 10, corresponding to half of a pressure scale height in the convection zone.

The mass entrained into the convection zone in the MHD case is smaller by 12% than that in the hydrodynamic setup. The partial suppression of the mixing at the convective boundary correlates with the average strength of the magnetic field in the convection zone. It is possible that the reduction of the kinetic energy of the convective flows caused by the action of the small-scale turbulent dynamo and the presence of magnetic fields with strengths up to 60% of the equipartition value at the convective boundary contribute to the partial suppression of the mixing. By the end of the simulations, the mass entrainment rate is reduced by 20% with respect to the hydrodynamic simulations.

In our simulations, SSDs seem to have only a mild effect on the growth of the convective oxygen shell. This is consistent with the findings of Varma & Müller (2021), who ran global MHD simulations of an oxygen shell in an 18 M_⊙ star. Overall, these authors do not observe any significant impact of the small-scale turbulent dynamo on the properties of the convective flows as compared to a non-MHD simulation. Those results, however, may have been affected by low effective Reynolds numbers, which are typical for global simulations of turbulent convection at moderate grid resolutions. Our “box-in-a-star” approach, used in combination with special numerical solvers optimized for tackling stratified, subsonic magnetoconvection, allows us to achieve much higher effective resolution than that obtained by Varma & Müller (2021). We find that the small-scale turbulent dynamo reduces the kinetic energy content of the convective shell by 25% on average and significantly changes the topology of the velocity field with respect to the purely hydrodynamic problem. These results are particularly important in the context of the “perturbation-aided” explosion mechanism, whose efficiency is set by both the magnitude of the convective velocities and the typical spatial scales of convection in the burning shells of the supernova progenitor (Müller et al. 2017; Couch et al. 2020).

One point of concern is related to the usage of the Implicit Large Eddy Simulation (ILES) approach in this study, which gives rise to effective magnetic Prandtl numbers Pr_m = ν/η close to or even larger than unity. Such large values of Pr_m are an overestimation of the actual conditions found in oxygen-burning shells. Based on the analytic expressions of the viscous and resistive coefficients provided in Augustson et al. (2019), we estimate that realistic magnetic Prandtl numbers in the oxygen shell considered here range from Pr_m = 0.001 to Pr_m = 0.1. It is well known that certain properties of the small-scale turbulent dynamo are sensitive to the value of Pr_m. In particular, the strength of the saturated magnetic field is often found to be an increasing function of Pr_m (Schekochihin et al. 2004; Brandenburg 2011; Käpylä et al. 2018). Unfortunately, no general consensus has been reached so far as to the behavior of SSDs at low Pr_m, which is likely setup dependent. More MHD simulations of the oxygen-burning phase with realistic Prandtl numbers are therefore required in order to properly establish the impact of magnetic fields on the dynamics of the convective shell. Finally, we did not include a nuclear-burning network in our simplified setup, so further investigation on possible indirect effects of magnetic fields on the nuclear energy generation and nucleosynthesis is certainly needed.

Acknowledgments

The work of G.L. and F.K.R. was supported by the German Research Foundation (DFG) through the grant RO 3676/3-1. P.V.F.E. was supported by the US Department of Energy through the Los Alamos National Laboratory (LANL). LANL is operated by Triad National Security, LLC, for the National Nuclear Security Administration of the US Department of Energy (Contract No. 89233218CNA000001). This work has been assigned a document release number LA-23-27290. G.L., R.A., J.H., and F.K.R. acknowledge support by the Klaus Tschira Foundation. The authors would like to thank Hirschi R., Rizzuti F., and Varma V. (Keele University) for fruitful discussions and hosting a visit to Keele University.

References

All Tables

All Figures

	Fig. 1. Profiles of density, pressure, pseudo-entropy (p/ργ), gravity, mass fraction of a passive tracer (ψ), and heat source term ($\dot{q}$) as a function of the vertical coordinate y at t = 0 s. Here, ρref = 1.82 × 106 g cm−3, pref = 4.64 × 1023 dyne cm−2, gref = 6.37 × 108 cm s−2, ${\dot{q}}_{\mathrm{ref}}=1.76\times {10}^{20}$ erg cm−3 s−1, and Lref = 4 × 108 cm.
In the text

Fig. 2. Development of convection in the MHD simulation of the idealized oxygen shell run on a 5123 grid. The panels show fluctuations in pseudo-entropy (A = p/ργ) in the z = 0 plane at different times, as indicated by the insets. The entropy generated by the heat source at the base of the box (see Fig. 1) is mixed throughout the initially adiabatic layer by turbulent convection. This process slowly increases the entropy content of the convection zone in time. The broad stripe of negative entropy fluctuation visible in the upper half of the domain at early times is due the thermal expansion of the convective layer. The turbulent flows also excite IGWs at the upper convective boundary (lower center panel) which then propagate in the subadiabatic layer.

In the text

	Fig. 3. Time evolution of the mean magnetic energy density inside the convection zone for the indicated grid resolutions.
In the text

Fig. 4. Transfer functions extracted from the convection zone on the 5123 grid and averaged over the kinematic stage of the dynamo. The panel on the left shows the transfer rates from the kinetic (K) to the magnetic (B) energy reservoir, whereas the panel on the right shows magnetic-to-kinetic energy transfer rates. The time averaging was performed such that, at each time t, all the transfer curves were rescaled by the maximum value of TKBS(kh). The resulting curves were then averaged over the time interval t ∈ (1.5τconv, 3τconv). Dashed lines represent negative transfer rates, while solid lines are used for positive rates.

In the text

Fig. 5. Growth rate of the mean magnetic energy inside the convection zone (averaged over the kinematic phase of the dynamo) as a function of the grid spacing, Δx. For this analysis, we also simulated the kinematic phase of the dynamo on grids with 3843 and 6403 cells (Δx/Lref = 5.2 × 10−3, and Δx/Lref = 3.1 × 10−3, respectively). Error bars represent three standard deviations, computed over the time series. The expected theoretical scaling for SSD amplification (γ ∝ Re1/2 ∝ Δx−2/3) is represented by the black dashed line.

In the text

	Fig. 6. Time evolution of the mean magnetic-to-kinetic energy ratio inside the convective shell.
In the text

	Fig. 7. Time evolution of the mean kinetic energy density inside the convection zone, ${\stackrel{\sim }{E}}_{\mathrm{K}}$, for the purely hydrodynamic (vermilion) and MHD (light blue) simulations. The time evolution of the mean magnetic energy density in the MHD simulations, ${\stackrel{\sim }{E}}_B$, is also shown.
In the text

Fig. 8. Snapshots of the Mach number on the 5123 grid at t/τconv = 22. The panels on the left show results from the purely hydrodynamic simulation, while the panels on the right show the MHD simulation. The upper row shows vertical cuts taken in the z = 0 plane, while the lower row shows a horizontal cut through the y = 1.5 Lref plane. The white-dotted line indicates the initial position of the convective boundary.

In the text

Fig. 9. Vertical profiles of the horizontal ($V_{\mathrm{h}}=\sqrt{V_x^2+V_z^2}$) and vertical (Vy) velocity components averaged over t ∈ (15τconv, 25τconv). Here, light blue is used to indicate the quantities extracted from the MHD simulations, whereas vermilion is used for the hydrodynamic simulations. An estimate of the statistical uncertainty associated to the averaged profiles can be inferred by looking at the typical dispersion among the hydrodynamic simulations, which represent a set of three independent (and numerically converged) realizations of the oxygen shell.

In the text

	Fig. 10. Vertical profiles of the root-mean-square relative pseudo-entropy fluctuation in the MHD (light blue) and purely hydrodynamic simulations (vermilion), averaged over t ∈ (15τconv, 25τconv). The upper convective boundary is, on average, at the position of the peaks visible at y = (2.3 − 2.4) Lref.
In the text

	Fig. 11. Vertical profiles of the vertical enthalpy flux (ℱH) averaged over t ∈ (15τconv, 25τconv). Here, light blue is used to indicate the quantities extracted from the MHD simulations, whereas vermilion is used for the hydrodynamic simulations.
In the text

	Fig. 12. Magnetic (light blue) and kinetic (black-dotted for the hydrodynamic simulations, and vermilion for the MHD simulations) energy spectra computed in the y = 1.5 Lref plane and averaged over the time interval t ∈ (15τconv, 25τconv). Each panel shows the curves extracted from a different grid. The black-dashed lines indicate the Kolmogorov scaling law, $k_{\mathrm{h}}^{-5/3}$.
In the text

	Fig. 13. Distribution of By on the z = 0 plane (upper row) and the y = 1.5 Lref plane (lower row) at t/τconv = 22 for the indicated grid resolutions.
In the text

	Fig. 14. Vertical profiles of the horizontal ($B_{\mathrm{h}}=\sqrt{B_x^2+B_z^2}$) and vertical (By) magnetic field, averaged over the time interval t ∈ (15τconv, 25τconv).
In the text

	Fig. 15. Vertical profiles of the root-mean-square magnetic field divided by the equipartition field, averaged over the time interval t ∈ (15τconv, 25τconv).
In the text

	Fig. 16. Vertical profiles of the buoyancy (light blue) and Lorentz work (vermilion) averaged over t ∈ (15τconv, 25τconv) for the indicated grid resolutions.
In the text

	Fig. 17. Time evolution of the mass entrained from the stable layer into the convection zone, rescaled by the total mass contained in the stable layer at t = 0 s. Light blue is used for the MHD simulations, whereas vermilion is used for the purely hydrodynamic runs.
In the text