3D Simulations of Magnetoconvection in a Rapidly Rotating Supernova Progenitor

We present a first 3D magnetohydrodynamic (MHD) simulation of oxygen, neon and carbon shell burning in a rapidly rotating 16 M_sun core-collapse supernova progenitor. We also run a purely hydrodynamic simulation for comparison. After 180s (15 and 7 convective turnovers respectively), the magnetic fields in the oxygen and neon shells achieve saturation at 10^{11}G and 5 x 10^{10}G. The strong Maxwell stresses become comparable to the radial Reynolds stresses and eventually suppress convection. The suppression of mixing by convection and shear instabilities results in the depletion of fuel at the base of the burning regions, so that the burning shell eventually move outward to cooler regions, thus reducing the energy generation rate. The strong magnetic fields efficiently transport angular momentum outwards, quickly spinning down the rapidly rotating convective oxygen and neon shells and forcing them into rigid rotation. The hydrodynamic model shows complicated redistribution of angular momentum and develops regions of retrograde rotation at the base of the convective shells. We discuss implications of our results for stellar evolution and for the subsequent core-collapse supernova. The rapid redistribution of angular momentum in the MHD model casts some doubt on the possibility of retaining significant core angular momentum for explosions driven by millisecond magnetars. However, findings from multi-D models remain tentative until stellar evolution calculations can provide more consistent rotation profiles and estimates of magnetic field strengths to initialise multi-D simulations without substantial numerical transients. We also stress the need for longer simulations, resolution studies, and an investigation of non-ideal effects.


INTRODUCTION
In recent years, multi-dimensional effects during the advanced convective burning stages of massive stars have received significant interest for multiple reasons, and have been studied extensively by means of hydrodynamic simulations.It has been recognised that seed instabilities from convection play an important dynamical role in core-collapse supernova explosions of massive stars (Couch & Ott 2013;Müller 2015;Couch et al. 2015;Müller et al. 2017;Bollig et al. 2021;Vartanyan et al. 2022).There is also the question of whether convective boundary mixing by turbulent entrainment and shell mergers may lead to structural changes in the pre-collapse structure of supernova progenitors and compared to current spherically symmetric stellar evolution models (e.g., Young et al. 2005;Meakin & Arnett 2007;Müller et al. 2016;Jones et al. 2017;Cristini et al. 2017Cristini et al. , 2019;;Kaiser et al. 2020;Rizzuti et al. 2023) and affect the nucleosynthesis outcomes from massive stars (Ritter et al. 2018).Finally, multi-dimensional simulations of late-stage convective burning are starting to shed light on angular momentum transport and magnetic evolution inside massive stars (Yoshida et al. 2021;Varma & Müller 2021;McNeill & Müller 2022), which is particularly relevant for our understanding of neutron star birth spin rates and magnetic ★ E-mail: v.r.vejayan@keele.ac.uk fields and hyperenergetic supernova explosions that are probably driven by rotation and magnetic fields.
Three-dimensional simulations of convection during advanced burning stages in massive stars have so far largely disregarded two important aspects of real stars -rotation and magnetic fields.The effects of rotation had only been touched upon by the seminal work of Kuhlen et al. (2003), and studies in axisymmetry (2D) of Arnett & Meakin (2010); Chatzopoulos et al. (2016), while 3D simulations have only started to explore rotation in recent years (McNeill & Müller 2022;Yoshida et al. 2021;Fields 2022).Similarly, magnetic fields during the advanced burning stages have only been considered recently by Varma & Müller (2021), and their work was limited to the non-rotating case.
Outside the context of advanced burning stages in massive stars, magnetohydrodynamic (MHD) simulations have been used extensively as a means to study convection over the years, primarily in the context of the Sun and solar-like stars (for reviews see, e.g., Brun & Browning 2017;Charbonneau 2020).Given the high quality of spatially and temporally resolved solar data (e.g., García 2009;Raouafi et al. 2023), these simulations often aim to explain more detailed observational time-varying features on the solar surface (Rincon & Rieutord 2018;Toriumi & Wang 2019;Petrovay 2020) and envelope convection, which includes understanding the formation of the solar rotation profile (e.g.Featherstone & Miesch 2015;Brun et al. 2017;Hotta et al. 2022;Camisassa & Featherstone 2022).Studies of stars more massive than the sun are currently limited to just a handful of core dynamo simulations of A and B-type stars (e.g., Brun et al. 2005;Featherstone et al. 2009;Augustson et al. 2016).Recently, 3D simulations (Raynaud et al. 2020;Masada et al. 2022) have also shown that the convective dynamo can be an important ingredient in amplifying magnetic fields in newly forming neutron stars.
Simulations of magnetoconvection during the late burning stages, both in rotating and non-rotating stars, are a necessity for several reasons.Even in slowly rotating massive stars, magnetic fields have been shown to impact the dynamics of the subsequent neutrinodriven explosions (Obergaulinger et al. 2014;Müller & Varma 2020;Matsumoto et al. 2020Matsumoto et al. , 2022;;Varma et al. 2023).For the magnetorotational explosion scenario (e.g., Burrows et al. 2007;Winteler et al. 2012;Mösta et al. 2014;Mösta et al. 2018;Obergaulinger & Aloy 2020a,b;Kuroda et al. 2020;Aloy & Obergaulinger 2020;Powell et al. 2023; see also early work on explosions driven by millisecond magnetars, e.g., Usov 1992;Duncan & Thompson 1992), a better understanding of the interplay between convection, rotation, and magnetic fields in supernova progenitors is even more critical.In this mechanism, a rapidly rotating core and very strong initial magnetic fields are required to launch the very energetic explosion.Such magnetorotational explosions are thought to explain rare, unusually energetic "hypernovae" with energies of up to ∼10 52 erg (Woosley & Heger 2006).
The magnetorotational explosion mechanism is linked to the problem of rotation and magnetism in massive stars.Initial conditions for magnetorotational explosion simulations currently come from "1.5D" stellar evolution models that assume shellular rotation and include effective recipes for magnetic field generation and angular momentum transport by hydrodynamic and magnetohydrodynamic processes (Heger et al. 2000(Heger et al. , 2005;;Woosley & Heger 2006).
There are still many open questions about the treatment of rotation and magnetic fields in stellar evolution models.Aside from purely hydrodynamic instabilities (Heger et al. 2000(Heger et al. , 2005;;Maeder et al. 2013), the interaction of rotation and magnetic fields is a critical issue.Since convective regions are usually assumed to rotate rigidly (as this is the only allowed rotational state in thermal equilibrium; Landau & Lifshitz 1969), attention has usually focused on angular momentum transport and dynamos in non-convective regions.The dynamo mechanism often implemented in these 1.5D stellar models to generate magnetic fields relies on sufficiently strong differential rotation in convectively stable regions of the star to stretch poloidal magnetic fields into toroidal fields.The dynamo loop is closed by the development of a pinch-type (Pitts-Tayler) instability (Tayler 1973;Pitts & Tayler 1985).This mechanism developed by Spruit (1999Spruit ( , 2002) ) is often referred to as the Tayler-Spruit dynamo.Recently, Fuller et al. (2019) have tried to improve the Tayler-Spruit dynamo mechanism, arguing that the Tayler instability saturates via turbulent dissipation of unstable magnetic field perturbations.This mechanism has a smaller energy dissipation rate and thus allows for stronger magnetic fields and more efficient angular momentum transport than the traditional Tayler-Spruit dynamo.
Other attempts to understand magnetic stellar evolution models have been to derive scaling relationships in convective regions (Jermyn et al. 2020a,b) and to explore the role of the magnetorotational instability (MRI) Wheeler et al. (2015); Takahashi & Langer (2021), driven, in part, by global 3D simulations such as that in Braithwaite & Nordlund (2006); Jouve et al. (2015); Meduri et al. (2019).These simulations have suggested that the interaction between the different instabilities and flows can be quite intricate, and may induce not only the pinch instability but can also be strongly affected by MRI and magnetic buoyancy.
Since 1.5D stellar evolution models implementing the Tayler-Spruit dynamo predict magnetic fields that are rather weak and predominantly toroidal, the general notion has long been that field amplification processes after the collapse are critical in magnetorotational explosions (Obergaulinger et al. 2009;Sawai et al. 2013;Sawai & Yamada 2016;Masada et al. 2015;Guilet & Müller 2015;Guilet et al. 2015;Mösta et al. 2015;Rembiasz et al. 2016a,b;Raynaud et al. 2020;Reboul-Salze et al. 2021), although this has recently been challenged (Obergaulinger & Aloy 2017, 2021;Aloy & Obergaulinger 2020).In particular, for sufficiently strong seed fields in the progenitor, the initial field strengths and geometry could have a significant impact on the development of magnetorotational explosions after collapse (Bugli et al. 2020;Aloy & Obergaulinger 2020), making an understanding of the pre-collapse magnetic fields in 3D indispensable.
In this study, we present a first simulation of rotating magnetoconvection during the final phases of shell burning using the ideal MHD approximation.This simulation constitutes a first step beyond spherically symmetric prescriptions in stellar evolution models to predict the magnetic field strength and geometry, as well as its role in angular momentum transport encountered in the inner shells of massive stars at the pre-supernova stage.We also compare to a corresponding non-magnetic model of the same progenitor to gauge the feedback of magnetic fields on the convective flow and rotation profiles.
Our paper is structured as follows.In Section 2, we describe the numerical methods, progenitor model, and initial conditions used in our study.The results of the simulations are presented in Section 3. We first focus on the strength and geometry of the emerging magnetic field and then analyse the impact of magnetic fields on the convective flows and rotation, with a focus on the turbulent mixing and angular momentum transport within and between the burning shells.We summarise our results and discuss their implications in Section 4.

NUMERICAL METHODS AND SIMULATION SETUP
We simulate oxygen, neon, and carbon shell burning with and without magnetic fields in a rapidly rotating 16 M ⊙ solar-metallicity helium star from Woosley & Heger (2006) with a strong differential rotation profile calculated using the stellar evolution code Kepler.The same progenitor model has previously been used in the PROMETHEUS rotating shell convection simulation of McNeill & Müller (2022).The structure of the stellar evolution model at the time of mapping to 3D is illustrated in Figure 1.
For our 3D simulations we employ the Newtonian magnetohydrodynamic (MHD) version of the CoCoNuT code as described in Müller & Varma (2020); Varma & Müller (2021).The MHD equations are solved in spherical polar coordinates using the HLLC (Harten-Lax-van Leer-Contact) Riemann solver (Gurski 2004;Miyoshi & Kusano 2005).The divergence-free condition ∇ • B = 0 is maintained using a modification of the original hyperbolic divergence cleaning scheme of Dedner et al. (2002) that allows for a variable cleaning speed while still maintaining total energy conservation as described in Varma & Müller (2021) (building on similar ideas by Tricco et al. 2016) where g is the gravitational acceleration,  t is the total (gas and magnetic) pressure, I is the Kronecker tensor,  h is the hyperbolic cleaning speed,  is the damping time scale for divergence cleaning, and  nuc and   are energy and mass fraction source terms from nuclear reactions.This system conserves the volume integral of a modified total energy density ê, which also contains the cleaning field ψ, where  is the mass-specific internal energy.
The simulations are conducted on a grid with 400 × 128 × 256 zones in radius , colatitude , and longitude  with an exponential grid in  and uniform spacing in  and .To reduce computational costs, we excise the non-convective inner core up to 3, 000 km and replace the excised core with a point mass.The grid extends to a radius of 40, 000 km and includes a small part of the silicon shell, the entire convective oxygen, neon and carbon shells.Our simulations cover the full sphere (4 solid angle).
In the MHD simulation, we impose a homogeneous magnetic field with   = 10 7 G parallel to the grid axis as initial conditions.
We implement reflecting and periodic boundary conditions in  and , respectively.For the hydrodynamic variables, we use hydrostatic extrapolation (Müller 2020) at the inner and outer boundary, and impose an effectively slip-free inner boundary.Different from the hydrodynamic simulations of Yoshida et al. (2021) and McNeill & Müller (2022), we do not contract the inner boundary to follow the contraction and collapse of the core.The model is instead meant to be representative of the physical principles governing late-stage, steady-state magnetoconvection with rapid rotation, rather than a pre-collapse model.
The inner and outer boundary conditions for the magnetic fields are less trivial.In simulations of magnetoconvection in the Sun, various choices such as vertical boundary conditions (  =   = 0), radial boundary conditions (  =   = 0), vanishing tangential electric fields or currents, perfect-conductor boundary conditions, or extrapolation to a potential solution have been employed (e.g., Thelen & Cattaneo 2000;Rempel 2014;Käpylä et al. 2020).Since our domain boundaries are separated from the convective regions by shell interfaces with significant buoyancy jumps, we opt for the simplest choice of boundary conditions and merely fix the magnetic fields in the ghost zones to their initial values for a homogeneous vertical magnetic field.We argue that due to the buffer regions at our radial boundaries, and the lack of rotational shear (due to the slip-free boundary conditions), our choice of magnetic boundary conditions should not have a significant impact on the dynamically relevant regions of the star.
Similar to the non-rotating magnetoconvection simulations done in Varma & Müller (2021), our models will not (and are not intended to) provide an exact representation of the pre-collapse state of the particular 16 M ⊙ star that we are simulating.We would expect, e.g., that for the particular 16 M ⊙ model, the burning rate and hence the convective velocities would increase until the onset of collapse due to the contraction of the convective oxygen shell.As a consequence of accelerating convection and flux compression, the magnetic fields will likely also be somewhat higher at the onset of collapse.The model is rather meant to reveal the physical principles governing late-stage magnetoconvection in rapidly rotating massive stars, and to be representative of the typical conditions in the burning shells with the understanding that there are significant variations in convective Mach number and shell geometry at the onset of collapse (Collins et al. 2018), which will also be reflected in the magnetic field strengths in the interiors of magnetorotational supernova progenitors.

Evolution of the magnetic fields
We simulate two rapidly rotating 16 M ⊙ models, one with and one without magnetic fields.The magnetic model is initiated with a homogeneous magnetic field of 10 7 G.We then allow the geometry of the magnetic field to evolve naturally under the influence of rapid rotation and convection.
The evolution of the root mean square (RMS), volume-averaged magnetic fields in the three convective burning shells we simulate -the oxygen, neon and carbon burning shell -are shown in Figure 2. We see an initial period of exponential growth of magnetic field strength in each shell before a plateau forms after ≈200s in each shell.The field strengths in the oxygen and neon shells both appear to follow a very similar trajectory, achieving a peak at ≈190s, before a gradual decline sets in.In each of the shells, convection takes a different amount of time to fully develop, which explains the slight delay from the start of the simulation to the beginning of the exponential field growth.In particular, the carbon shell has a much longer convective turnover timescale  c than the other two shells with an initial value  c ≈ 300 s compared to ≈15 s and 25 s for the oxygen and neon shell, respectively, which considerably delays the growth of Evolution of the volume averaged (solid) magnetic field strength within the oxygen (black), neon (blue) and carbon (red) shells.The dashed line approximates the expected growth of the magnetic fields strength in the oxygen shell via the Ω dynamo mechanism described in Section 3.1.The dotted line shows the expected magnetic field strength saturation due to differential rotation in the oxygen shell based on Equation 9. magnetic fields.The growth of the magnetic fields in the carbon shell already becomes apparent after ≈60 s, even without convection being fully developed.This is due to field amplification from strong differential rotation which develops at the base of the shell, and turbulent fluctuations that developed alongside the convective plumes.
Due to the development of convection in the shells, coupled with rapid differential rotation (maximum rotation rate of Ω ≈ 0.104 rad s −1 ), we expect the field amplification to be dominated by the Ω dynamo mechanism, which is often proposed as the mechanism that sustains the solar magnetic field (Charbonneau 2020).The mechanism stretches the poloidal magnetic fields into toroidal fields via differential rotation (Ω-mechanism), and the toroidal field is then stretched into a poloidal field due to convective motions (-mechanism), completing the cycle and amplifying the seed field.To test if this is the case, we plot the expected growth of the Ω dynamo for the oxygen shell in Figure 2.
To this end, we approximate the growth of the magnetic field via the Ω mechanism in the oxygen shell by assuming the magnetic field evolves via the simplified evolution equation: which has a solution for the magnetic field growth of the form  rms =  0  Γ Ω Δ , where  0 is the initial field strength.We take the growth rate of the Ω dynamo to be Γ Ω = (/)(Ω c ) 1/2 as presented in Vishniac (2005) based on dimensional arguments, where  is the convective velocity,  is the radial extent of the convective zone, Ω the rotation rate and   is the convective turnover timescale.Since  c ∼ /, this effectively amounts to a growth time scale  Ω of the order of the geometric mean of the rotation period  = 2Ω −1 and the convective turnover time, In the evolution plotted in Figure 2, we first calculate the RMS averaged values  and Ω in the oxygen shell, as well as angular averages the radial extent of the oxygen shell, .The convective turnover time is calculated using these averaged quantities,  c = /.These averaged quantities are used to then determine the growth rate, Γ Ω to be used at each time step, to evolve the magnetic field.
The expected growth rate of the Ω dynamo follows the growth of the magnetic field in the oxygen shell very closely for the first ≈140 s, after which time the magnetic field growth in the simulation slows down, and eventually decays after hitting a peak magnetic field strength of ≈ 10 11 G and 7 × 10 10 G for the oxygen and neon burning shells, respectively.We see that the expected growth rate from an Ω dynamo also decreases at later times due to two factors.First, convective velocities drop due to suppression of convection by strong magnetic stresses.Second the rotation rate drops due to large angular momentum fluxes.We will discuss these effects further in Sections 3.2 and 3.3.
These effects stop the magnetic field from being amplified via the Ω dynamo.Since the convection dies down, it is reasonable to expect that late-time field amplification and saturation is determined by differential rotation alone.Interestingly, the saturation of the field appears to be well described by an amplification mechanism that is driven by the MRI.For the MRI, Akiyama et al. (2003) argue that the saturation field is roughly given by where  is the density, Ω the rotation rate,  the radius and  ln Ω  ln  quantifies the amount of differential rotation present.Akiyama et al. (2003) derive Equation ( 9) by assuming that saturation of the magnetic field is achieved in the star when the characteristic mode scale  mode ≈  A (d ln Ω/d ln ) −1 is equal to the local radius , since the wavelength of the mode cannot be larger than the physical size of the unstable region.Here,  A is the Alfvén velocity ( A = / √︁ 4).On dimensional grounds, one can expect Equation (9) to hold not just specifically for the MRI, but more broadly for amplification mechanisms driven solely by differential rotation in the ideal MHD regime with negligible resistivity.Figure 2 shows that the magnetic field in the oxygen shell saturates at a very similar level to Equation ( 9).
The fastest-growing MRI mode has the wavenumber where    =  0  / √  is the Alfvén velocity in the vertical direction (Balbus & Hawley 1991;Rembiasz et al. 2016a).In our simulation, this corresponds to a wavelength of ≈ 10 6 cm at early times when the field growth is dominated by the Ω dynamo.After the field saturates, and we expect differential rotation to play a key role in maintaining the magnetic fields, we find  MRI ≈ 10 8 cm, which is well resolved in our simulated grid(dr ≤ 6×10 6 cm within the oxygen and neon shells).Although buoyancy forces can impact the MRI modes (Obergaulinger et al. 2009), the regions we are considering are well-mixed with a flat entropy gradient, so it is justified to ignore buoyancy.
The strong magnetic fields also result in very rapid redistribution of angular momentum, which slows the rotation rate of the oxygen and neon shells dramatically.Since the magnetic field saturation depends on the rate of rotation, this in turn leads to a drop in the average magnetic field strength by over 50% in these shells, as we see in Figure 2. The consequences of this will be discussed in more detail in the next sections.
Aside from the equilibrium field strength, it is worth investigating the field geometry that has naturally developed in the saturation state.To this end, we show radial profiles of the RMS averaged field strength and of the dipole field strength in Figure 3, as well as zonal averaged plots of B r , B Θ and B Φ in Figure 4.The dipole field is calculated by extracting just the ℓ = 1 component of the spherical Close to the end of the simulation, the RMS field strength appears almost flat throughout the simulated domain, varying only between ≈5 × 10 9 G and 10 10 G.The dipole component of the magnetic field reaches about one-third of the total field strength in the inner oxygen and neon burning shells, which are no longer convective at this point.Further out, however, the dipole is weaker by comparison to the RMS-averaged field.The slowly rotating and still convective carbon shell may be concentrated in smaller scale field structures that are similar to the non-rotating convective shell presented in Varma & Müller (2021).However, as the carbon shell has only completed 1-2 convective turnovers, it is difficult to say if this structure will be maintained at later times.
We attempt to better visualise the evolution of the magnetic field geometry in Figure 4 by plotting the zonal averages of the three field components, B r , B Θ and B Φ at two times, up to 15, 000 km.The left halves of each plot show the averages at ≈ 250s, soon after magnetic field saturation is first achieved and convection is suppressed in the Neon and Oxygen shells.The right halves are at the same time chosen for Figure 3, at ≈ 480s, near the end of the simulation.We find that during the convective dynamo phase, the field strength is concentrated in small-scale structures, but once convective and rotational velocities are reduced, the small-scale magnetic field structures appear to diffuse away, leaving a mostly dipolar magnetic field.Some smaller-scale structures are still present, however, it isn't clear if these will still be present at even later times.
The carbon shell behaves somewhat differently from its neighbours as the shell is much larger, and already more slowly rotating at the beginning of the simulation (Figure 1).Due to the slower rotation and density of the shell, the magnetic fields in the carbon shell saturate at a lower field strength.But unlike the two inner shells, the magnetic stresses here always remain below the radial kinetic stresses (Figure 5), so convection continues unimpeded by the magnetic fields, and coupled with differential rotation, sustains a relatively constant magnetic field strength.Unfortunately, due to the very long convective turnover times, we are only able to resolve about one convective turnover in this shell.Pushing this simulation further becomes untenable as the convection in the carbon shell has begun to interact strongly with our outer domain boundary.

Impact of magnetic fields on convection and energy generation
As we already briefly mentioned above, the amplification of the magnetic fields in our simulation leads to a very rapid suppression of convection, as well as fast transport of angular momentum out of the affected shells.Here, we attempt to understand the consequences of these dynamical changes by comparing the MHD simulation to a purely hydrodynamical simulation of this progenitor.
As the magnetic field grows, it eventually becomes strong enough to affect the bulk flow in the convection zones.To illustrate this, we compare the spherically-averaged diagonal components of the kinetic (Reynolds) and magnetic (Maxwell) stress tensors    and    .   and    are computed as where angled brackets denote volume-weighted averages.Note that we do not subtract the mean rotational flow for    here.In Figure 5 we present the stresses in the MHD model at ≈180 s, where the Maxwell stresses begin to be comparable to the radial Reynolds stress, and near the final time-step of the simulation at ≈480 s.
In Figure 5(a), the radial and meridional magnetic stresses are comparable to the radial kinetic stresses in the innermost regions of the star.This corresponds to pseudo-equipartition of these stress components throughout the oxygen and neon burning shells out to a mass coordinate of ≈2.8 M ⊙ .The magnetic stresses then generate backreaction against the convective flows in these shells.As shown in Figure 5(b) at a later time in the simulation, the backreaction greatly suppresses the convective velocities in these shells, lowering the radial kinetic stresses by several orders of magnitude.We plot the RMS angle-averaged radial kinetic energy near the end of the simulation (≈480 s) in Figure 6 for both the MHD model (top row) and for the purely hydrodynamic case (bottom row).At this time, the suppression of convective motions in the oxygen and neon shells by strong magnetic fields cause the radial kinetic energy in the inner 2.8 M ⊙ to be about three orders of magnitude lower than what is seen in the hydrodynamic model (≈10 16 − 10 17 g cm −1 s −2 compared to ≈10 19 − 10 20 g cm −1 s −2 ).As mentioned above, the carbon shell has only had time for about one convective turnover, i.e., convection is yet to reach a fully developed state.At the end of our simulation, the carbon shell is still convective, and the kinetic stresses in this shell remain higher than the magnetic stresses.
At early times, we see that the angular Reynolds stresses are the dominant components due to the rapid rotation.This along with the sharp gradients in these stresses that develop means that the shear instabilities can efficiently mix material more efficiently than in the underlying 1D stellar evolution models, where there is little mixing beyond the convective zones on dynamical time scales.Shear mixing outside the convective regions initially plays a significant role in the MHD model as well.The   stress component, which is indicative of radial motions that contribute to turbulent mixing, stays high outside the convective regions in the MHD model initially (Figure 5(a)).However, in the MHD model, the transport of angular momentum flattens the rotation profile (as discussed in detail in Section 3.3) , which we can see from the change in    and    , significantly reducing the shear mixing at late times.In our purely hydrodynamic model, however, the Reynolds stresses remain roughly the same throughout our simulation, leading to continuous enhanced mixing compared to the MHD counterpart.We find that this enhanced shear mixing compared to the initial expectation from the 1D stellar evolution model means that the burning occurs at different regions, outside the initial location of the convection zones.
The consequence of this can be seen by analysing how the mass fractions of several key elements evolve.In Figure 6, we compare the mass fractions of silicon, oxygen and neon in the MHD simulation and the purely hydrodynamic simulation.We also plot the radial kinetic energy at the end of the simulations, to further stress the differences in turbulent mixing when magnetic fields are introduced.The plots are limited to the inner 3 M ⊙ of the enclosed mass to focus on the oxygen and neon burning shells, which are most strongly affected by magnetic fields.
The radial profiles of all three elemental mass fractions in the hydrodynamic and MHD model evolve very similarly at early times (up to ≈120 s in Figure 6), when the magnetic fields are not strong enough to significantly affect mixing.At later times, however, the differences in mixing become quite apparent.The plots of the silicon mass fraction show that large fractions of silicon are mixed outwards to an enclosed mass of ≈2.65 M ⊙ in the hydrodynamic model, while there is little mixing beyond ≈2.5 M ⊙ in the MHD simulation.
The inhibition of mixing also means that material that gets burnt at the base of a shell is less efficiently replenished by fresh material or not replenished at all.This effect is seen in the oxygen shell as it burns material in a relatively narrow region around 1.85 M ⊙ .In the MHD simulation, a sharp drop in the oxygen and neon mass fractions develops in this region, which gets steeper at later times due to the rapid burning (mostly of oxygen), which is no longer replenished by convective and shear mixing.It is clear that this is caused by the very strong magnetic fields and a significant reduction in turbulent mixing by comparing the mass fractions to the purely hydrodynamic simulation.Without magnetic fields, there is no sharp drop in mass fraction at the bottom of the burning shell.The profiles remain smoother within the shell and are even smoothed beyond the boundaries by shear mixing.The oxygen and neon mass fractions in the oxygen burning region even increase over time, as the rapid We also find reduced mixing in the neon burning shell between ≈2.20 M ⊙ and 2.40 M ⊙ .Here, instead of a sharp drop, we see the steep gradient of neon, where neon is consumed, move outwards after the initial transient where material is mixed inwards.This has the effect of shifting the location of the base of the burning shell.We show this phenomenon more clearly in Figure 7, which presents angleaveraged radial profiles of the energy generation rate at 110 s, 220 s and 450 s (dotted, dashed and solid lines, respectively).We show their evolution in both mass and radial coordinates in Figures 7(a) and (b) respectively, from the inner boundary of our simulation 1.75 M ⊙ to an enclosed mass of 3.0 M ⊙ , and radius of 3000 km to 12000 km for both the hydrodynamic model and the MHD model.The energy generation profiles for both models show three clear peaks, which correspond to the three burning shells (oxygen, neon and carbon).We see that both simulations quickly deviate from each other, with the energy generation rate in the oxygen and neon shells in the hydrodynamic model receding backwards, and getting stronger at later times, while the opposite is seen in the MHD simulation.We note that the relative change in the position of these energy generation peaks are largely similar in both mass and radial coordinates.
Figure 7 clearly shows that the second energy generation peak in the MHD model, which corresponds to neon shell burning, moves from ≈2.35 M ⊙ (≈5400 km) at early times to ≈2.50 M ⊙ (≈5800 km) later on, compatible with the change in the neon mass fraction profiles in Figure 6.Due to the lack of turbulent mixing, the peak of neon burning moves radially outward as neon gets burnt at its initial position.However, since the neon burning rate is extremely temperature-sensitive (∝  50 , Woosley et al. 2002), the energy generation rate drops significantly when burning is moved to a slightly cooler region.
We see a similar effect for the oxygen shell.At early times, before the magnetic field heavily suppresses turbulent mixing, both the MHD and hydrodynamic model mix material rapidly.The convective boundary mixing (CBM), particularly at the lower boundary of the oxygen shell, entrains material, increasing the size of the oxygen shell, causing it to move inwards (both in mass and in radius).We see the peak energy generation move from ≈1.95 M ⊙ (≈3900 km) to ≈1.85 M ⊙ (≈3500 km).
After turbulent convection is suppressed in the MHD simulation, however, the nuclear energy generation peak in the oxygen shell behaves similarly to the neon burning shell, i.e., it moves outward in mass and radius to where oxygen is burned at a lower temperature.The hydrodynamic model, however, continues to entrain material from beneath the oxygen shell, moving the peak of its energy generation deeper into the core.oxygen burning is also a very temperaturesensitive reaction (∝  33 , Woosley et al. 2002), so these shifts in the burning shells lead to noticeable changes in the total energy generation rate over time.The consequences of the suppression of mixing in the MHD model are seen more clearly in Figure 8.Here, we plot the total (volumeintegrated) energy generation in the oxygen, neon and carbon burning shells with time for both the hydrodynamic and MHD models.As we see in Figure 7, the lack of mixing in the MHD model leads to the location of the oxygen and neon shells moving radially outwards with time, to cooler regions of the star.This causes a gradual decrease in energy generation after the magnetic fields in these shells first achieve saturation at ≈200 s.The increased mixing in the hydrodynamic model, on the other hand, causes a subsequent increase in peak energy generation rate as the shells move deeper into the star.
Interestingly, after ≈220 s of similar energy generation in the carbon shell in both models, the energy generation rate starts to increase in the MHD model and decreases in the hydrodynamic case.This is likely due to the slight change in the position of the carbon burning shell for the hydrodynamic model while it remains mostly stationary in radius in the MHD model.From the profiles in Figure 7, this is likely caused by the radial expansion of the neon shell in the hydro-dynamic model as its energy generation rate increases, pushing the carbon shell further outwards.In the MHD model, by contrast, the energy generation rate in the neon shell has dropped over time and hence there the carbon shell is not driven outward.

Evolution of rotation and angular momentum transport
In addition to its effect on turbulent mixing, the development of strong magnetic fields leads to rapid redistribution of angular momentum.The evolution equation for the angular momentum can be obtained by taking the cross product of the position vector r with the fluid momentum equation (Shu 1992;Mestel 1999).When including magnetic stresses in the momentum equation and integrating over spherical shells, one obtains (cf.Charbonneau & MacGregor 1993;Spruit 2002), where ∇  is the radial component of the divergence operator and angled brackets denote averages over solid angle.We then perform a Reynolds/Favre decomposition (Favre 1965) around a base state with constant angular velocity Ω , Ω = ⟨Ω   2 sin 2 ⟩ ⟨ 2 sin 2 ⟩ = ṽ   sin  ĩ , on spheres, as in McNeill & Müller (2022).Here ĩ = ⟨ 2 sin 2 ⟩/ ρ.Note that we use hats and primes for volume-weighted Reynolds averages and their fluctuating components: where  = sin  d d is the solid angle element.We denote massweighted averages and fluctuations with tildes or angled brackets and double primes: Applying the usual rules for Favre averages, ⟨⟩ = ρ X, ⟨ X Ỹ ⟩ = ρ X Ỹ and ⟨ X ′′ ⟩ = 0, we get We see from the decomposed angular momentum Equation ( 15) that angular momentum is transported by four distinct flux terms.In the order they are listed above we have an advective term, a meridional circulation term and turbulent transport, as well as an additional magnetic stress term.From this equation, aside from the usual hydrodynamic terms, the radial flux of angular momentum also depends on the strength of the magnetic fields.
The evolution of the angle-averaged rotation rates in the MHD model and purely hydrodynamic model are depicted in Figure 9.The top row shows the rotation rates, Ω, over time for both models, and the bottom row presents angle-averaged specific angular momenta.We note that in comparison to its hydrodynamic counterpart, as expected, the MHD model exhibit more efficient angular momentum transport.At the innermost region of our domain, the rotation rate is lowered by over an order of magnitude due to outward transport of angular momentum, from over 0.10 rad s −1 initially to ≈0.01 rad s −1 ).The rotation profile is flattened considerably.The angular momentum is taken up by carbon shell outside 2.8 M ⊙ , whose rotation rate increases.We note that Figure 9 only shows a limited inner portion of the total simulated carbon shell.
The hydrodynamic model displays much less smooth rotation profiles than the MHD model.It is noteworthy that at the various convective shell boundaries, the rotation profile shows significant dips.Figure 9(b) shows dips in rotation at ≈1.85 M ⊙ , which is the base of the oxygen shell, between the oxygen and neon shells at a location that varies over time between ≈2.20 M ⊙ and ≈2.40 M ⊙ , and finally between the neon and carbon shells at ≈2.80 M ⊙ .It is particularly interesting that some of these dips even reach negative values, i.e., there are shells with net retrograde rotation, which remain quite stable.On closer inspection of the rotation profile of the MHD model, we see that these dips in the rotation profile also begin to form early on in this case (Figure 9 (a) at 120s).However, the rotation profile is quickly smoothed once angular momentum transport due to magnetic fields become efficient.
Unlike in the MHD model, where angular momentum is simply transported outward, the redistribution of angular momentum in the hydrodynamic model appears more complicated (Figure 9 (d)).Effectively, positive angular momentum is transported into convective shells from the shell boundaries.This leads to an interesting non-monotonic rotation profile, where the fastest rotation rate is not reached at the inner boundary, but instead at 2.4-2.5 M ⊙ at late times, as well as the aforementioned dips between the convective shells.Redistribution of angular momentum also increases the rate of rotation in the inner carbon shell (around 2.8 M ⊙ ) and induces strong (radial) differential rotation there.
To better understand the emerging rotation pattern in the hydrodynamic model and the counterintuitive phenomenon of retrograde rotation, we consider zonal and temporal averages of the meridional velocity and the rotation rate in Figure 10.In Figure 10  The left halves of Figure 10 (a) and (b) show the zonal flow averages of the quantities, while the right halves show snapshots of the two quantities on a meridional slice.From Figure 10(b), it is clear that the retrograde rotation occurs mostly at the base of the carbon shell, and to a lesser extent, at the base of the oxygen shell.The snapshot on the right half of Figure 10(b) shows that regions of retrograde rotation form at the base of the neon shell as well, but these are more transient and do not show up in the zonal averages.We also observe that for the hydrodynamic model, although the retrograde rotation appears to form as a shell at the base of the carbon shell, in general, the rotation pattern that forms is not shellular.At the base of both the oxygen and carbon shell, we rather see indications of anti-solar differential rotation with faster rotation near the poles.
Such a rotation pattern has been observed before in simulations of surface convection zones.Featherstone & Miesch (2015) attribute the development of different rotation profiles, in part, to a link between the differential rotation and the meridional circulation in the convection zone.Anti-solar rotation profiles, are attributed to inward angular momentum transport, which establishes a single-celled meridional circulation profile throughout the convection zone.This circulation transports angular momentum polewards, spinning up the poles relative to the equator.Although our model exhibits antisolarlike rotation profiles at the base of the carbon (8400 km) and oxygen (≈3400 km) shells in Figure 10(b), the corresponding meridional circulation in Figure 10(a) does not exhibit the single-celled structure expected from Featherstone & Miesch (2015).We instead find that the meridional flows are not clearly structured, but appear more similar to a case with multi-celled circulation.
We see an analogous meridional circulation flow develop in the inner regions of the carbon shell of our model (above 8400km), with two large cells of material (Figure 10(a)).One noticeable difference shown in our model is that the circulation velocities do not drop off towards the poles, instead, they remain very strong.This may be part of the reason why the model exhibits a stable retrograde shell of material at 8400 km that extends across almost all latitudes as shown in Figure 10(b), rather than being confined to the equatorial region as in the surface convection models presented in Featherstone & Miesch (2015).
One major difference that our simulations have that likely alters the meridional flows, is the proximity of other convective shells.While in surface convection simulations, the shell is bounded only by a radiative core, in our model, we have three adjacent convective shells.Although these shells are initially separated by thin radiative zones in the 1D stellar structure input model, the rapid rotation and turbulent convection increases the amount of mixing that occurs at the convective shell boundaries in the 3D model, causing the convective shells to start interacting.This adds additional complications to the transport of angular momentum, and in turn, to the meridional circulation.
From extensive studies in solar physics, it is often expected that the stable differential rotation in our Sun is maintained by rotating turbulent convection (e.g., Ruediger 1989;Hotta et al. 2022), due to the interplay of buoyancy and inertial forces.We suspect that, similarly, the complicated rotation profile developed in our purely hydrodynamic model occurs due to an interaction between the Coriolis and buoyancy forces.Similar retrograde rotation patterns have been found in studies of rotating solar-like stars, and depend on the Rossby number of the flow (e.g., Brown et al. 2008;Guerrero et al. 2013;Brun et   strong enough to make such an effect plausible, we analyse how the Rossby number compared to the MHD model, which develops a flat rotation profile.We define the Rossby number of any given burning shell as where  conv is the convective velocity, Ω shell , is the rotation rate, and Δ is the radial extent of each burning shell.The time evolution of the Rossby number in the oxygen, neon and carbon shell is shown in Figure 11.We find that the Rossby number for the oxygen and neon shell for both models starts at ≈0.4 (and remains largely unchanged for the hydrodynamic model), indicating that the flows are strongly shaped by the Coriolis force.The Rossby numbers in the MHD model clearly reflects that the magnetic fields starts to strongly alter the bulk flows.At ≈180s the magnetic fields become strong enough to start rapidly transporting angular momentum outwards, slowing the rotation of the star and hence increasing the Rossby number.The Rossby number continues to increase until ≈200 s in the neon shell and ≈240 s in the oxygen shell, when the suppression of convective flows by magnetic stresses becomes dominant, lowering the Rossby number.
The carbon shells of both models clearly have not reached a convective steady state and hence not even a transient steady state in the Rossby number can be discerned.The MHD case initially transports angular momentum out of the star more efficiently.We see in Figure 9(c) that angular momentum from the rapidly rotating inner regions (inner 2.8 M ⊙ ), is transported out to the carbon shell (2.8 − 3.0 M ⊙ ).Focusing on the lines at 215 s in Figure 9(a), we see that at these times, rotation rate in the inner carbon shell increases in the MHD case compared to the purely hydrodynamic model.This leads to the deviation in Rossby number initially seen from ≈150 s, causing the Rossby number in the MHD model to be lowered compared to the hydrodynamic model.After ≈260 s, the Rossby number in the hydrodynamic model starts decreasing gradually.We associate this with the decrease in energy generation compared to the MHD model seen after ≈220 s in Figure 8, which would have the effect of decreasing the convective velocity in the carbon shell.The opposite trend is seen  2017), retrograde rotation is sometimes seen when buoyancy forces dominate the flow (i.e.Ro ≥ 1), where solar-like stars develop retrograde rotation at the equators and faster rotation at the poles (anti-solar rotation), and faster rotation at the equator for Ro ≤ 1. Figure 11 shows that although we develop strong retrograde rotation in the hydrodynamic model, the Coriolis force dominates its flow in each shell (i.e.Ro ≤ 1).
We note, however, a key difference between these models.The anti-solar rotation profiles in low-mass stars develop throughout the convective zone, whereas we find the retrograde motion to be largely concentrated at the base of the burning shells.The different phenomenology could be explained by the rather disparate conditions in surface convection zones in solar-like stars and convection during advanced burning stages.Solar-like stars have an inner radiative zone and a single convective shell above it, and radiation diffusion plays a role both for the internal structure of the convection zone and especially for the structure of the convective boundaries.Our progenitor has three interacting convective shells, radiative effects are unimportant (and not included in the model), and the structure of the boundary is determined by turbulent entrainment.
We suspect that the cause of the retrograde rotation in our simulation is similar to what is described in Aurnou et al. (2007) from hydrodynamic simulations of ice giants, and for solar-like simulations in Camisassa & Featherstone (2022), where convective rolls exhibit a preferred "tilt" in the positive  direction.The tilted flow structures create a correlation between flows moving inward (outward) and those moving in the negative (positive)  direction.Due to the strong turbulent mixing between shells, however, it is difficult to see in our models whether this tilted flow structure truly arises.For low Rossby convection in the solar case, this usually results in net angular momentum transport away from the rotation axis, which tends to speed up the equator.Since the usual Rossby number characterises the flow in the convection zone globally, it alone would not give us insight into dynamics at the convective boundaries or between convective shells.To understand the interplay between the Coriolis and buoyancy forces at the convective boundaries, we instead plot the angle-averaged magnitude of the buoyancy and Coriolis forces (  B and  C ) per unit mass of the hydrodynamic model in Figure 12.
Here  is the gravitational acceleration, ρ the RMS averaged density,  is the RMS averaged fluctuation from the average density,  is the angular velocity vector of the rotation and v is the velocity.For simplicity, we assume rotation is confined to   as in our initial conditions, i.e.,  points in the -direction, giving us only a component for  C pointing away from the rotation axis.We plot the RMS average of the absolute value of  C in Figure 12.This is to allow for greater clarity in comparing the two forces, since the retrograde rotation would lead to regions of negative   .We plot these forces at 57 s, 115 s and 190 s (dotted, dashed and solid lines, respectively).These times were chosen to represent approximate times before the development of retrograde rotation in the hydrodynamic model, when retrograde motion initially begins at the base of the carbon shell, and when it begins at the base of the oxygen shell.
For our simulation, we find the convective boundary overshooting at lower convective boundaries leads to buoyancy forces dominating the Coriolis force in between convective shells (≈1.85 M ⊙ and 2.8 M ⊙ , see Figure 12), which is likely why retrograde motion in our models are confined to the shell boundaries.This is further supported by the fact that we see the "stable" retrograde shell start to form at the base of the oxygen shell (≈1.85 M ⊙ ) when the buoyancy force surpasses the Coriolis force.The evolution of the Coriolis force in Figure 12 shows that this effect develops due to the transport of angular momentum away from convective boundaries into the convective shells.The magnetic model initially shows a similar force ratio, however, the ratio of buoyancy to Coriolis forces is soon reduced by the rapid rise of magnetic field strength and the subsequent suppression of convection, and hence buoyancy force.

CONCLUSIONS
We investigated the evolution of magnetic fields during advanced convective burning stages in massive stars and their backreaction on the flow in a simulation of the oxygen, neon and carbon shells in a rapidly rotating 16 M ⊙ progenitor shortly before core collapse.For comparison, we conducted a purely hydrodynamic simulation of the same progenitor as well.The simulations were run for about 8 minutes of physical time, corresponding to about 32 convective turnovers in the oxygen shell.
Rapid differential rotation and convection initially amplify the magnetic fields exponentially via the Ω-dynamo.However, strong magnetic stresses eventually dominate the radial kinetic stresses.The backreaction of the fields on the flow stops the exponential growth and suppresses convection in the oxygen and neon burning shells.These shells are effectively turned into convectively stable shells by the strong magnetic stresses and continue to burn fuel at the base of the shells without mixing of fuel and ashes.The magnetic field reaches saturation in the oxygen and neon shells after 180 s (corresponding to ≈12 and 7 convective turnovers respectively).It peaks at 10 11 G in the oxygen shell but decays to 3 × 10 10 G by the end of the simulation.In the carbon shell, the field appears to saturate at 10 10 G, but this shell has only completed about one turnover during our entire simulation, so a steady state likely has not been reached.
The strong magnetic fields that develop also transport angular momentum much more efficiently than in the purely hydrodynamic model.Already within the short duration of this simulation, the structure transitions from strong differential rotation into a nearly uniform rotation profile, with significant spin-down of the inner shells.
In the purely hydrodynamic model shell convection is sustained and strong differential rotation is maintained.However, it develops a much more complicated rotation profile than in the underlying 1.5D stellar evolution models.
The emerging rotation profile shows sharp drops at convective boundaries and, during some phases, even shells with retrograde rotation.We hypothesise that this is due to an instability that occurs when rapid changes of the rotational and convective velocities occur at the convective boundaries, coupled with strong meridional flows towards the poles.The regions with retrograde rotation are conspicuously associated with spikes in the local Rossby number, i.e., the ratio of the RMS-averaged buoyancy and Coriolis force.To better understand this phenomenon, we require further studies with different progenitors, varying Rossby number flows, and different grid geometries.
The transition of the oxygen and neon shells to slowly and rigidly rotation, non-convective region significantly reduces turbulent mixing.While the hydrodynamic model rapidly mixes new material into regions that are burning, the MHD model exhibits sharp drops in oxygen and neon mass fractions in the narrow burning regions.One consequence of this difference is that the hydrodynamic model entrains material deeper into the star, moving the peak of the energy generation of the oxygen shell radially inwards, while the location of peak energy generation of the same shell in the MHD simulation moves outwards.Due to the strong temperature sensitivity of oxygen burning (∝  33 ), this small change in shell position leads to a noticeable change in energy generation between the two models, resulting in increasing nuclear energy generation in the hydrodynamic model and inhibition of nuclear energy generation in the MHD model at late times.
Our results have important implications for core-collapse supernova modelling.For this particular rotating progenitor model, we predict pre-collapse fields of ≈2 × 10 10 G in the oxygen shell, similar to what we find for the non-rotating case in Varma & Müller (2021).Our rotating model exhibits a more gradual drop in field strength with radius.With relatively strong seed fields, we expect less of a delay until magnetic fields can contribute to become relevant for shock revival, i.e., by providing an additional "boost" to neutrino heating, as seen in Müller & Varma (2020); Varma et al. (2023).Due to the suppressed convective flows, the perturbation-aided mechanism (Couch & Ott 2013;Müller 2015) may be less effective, however, asymmetries seeded by the strong magnetic fields may be enough to deliver a similar effect (Varma et al. 2023).Perhaps most importantly, the very rapid redistribution of angular momentum transport from the inner shells casts doubt on the viability of a fast magnetorotational explosion powered by a "millisecond magnetar".For the right conditions to develop, a mechanism would be required to spin up the proto-neutron star during or after the core collapse for a magnetorotational explosion to be launched.However, there is still work to be done until the findings from simulations of magnetoconvection in rotating stars can be incorporated into models of magnetically-or magnetorotationally-driven explosions.For example, future simulations will need to include the core and self-consistently follow its contraction and incipient collapse to provide initial conditions for supernova simulations.
However, multi-D simulations of rotating massive stars face a much more fundamental challenge.The MHD model, and to some extent the hydrodynamic model, rapidly diverges from the initial structure of the stellar evolution model.Current stellar evolution models are clearly far from the actual quasi-steady state conditions that would emerge under the influence of rotation, convection and magnetic fields.Ideally, 3D simulations should cover significantly longer time scales to follow the relaxation of the structure into equilibrium and then study the subsequent evolution on secular time scales, but this is clearly beyond current computational resources.It is therefore very important to make 1D stellar evolution models and MHD models more consistent with each other to minimise deleterious effects from big initial transients that limit the fidelity of 3D simulations.This will require improved formalisms for stellar evolution with rotation and magnetic fields (Takahashi & Langer 2021).Developing the appropriate methodology for solving the problem of stellar evolution with rotation and magnetic fields by a combination of 1D and 3D modelling is bound to remain an extraordinary and exciting challenge.
Figure2.Evolution of the volume averaged (solid) magnetic field strength within the oxygen (black), neon (blue) and carbon (red) shells.The dashed line approximates the expected growth of the magnetic fields strength in the oxygen shell via the Ω dynamo mechanism described in Section 3.1.The dotted line shows the expected magnetic field strength saturation due to differential rotation in the oxygen shell based on Equation9.

Figure 3 .
Figure 3.The angle-averaged RMS value (black) and the dipole component (red) of the radial magnetic field component as a function of mass coordinate at a time of 480 s.

Figure 4 .Figure 5 .
Figure 4. Zonal averages of the magnetic field components, B  , B  and B  , shown in panels (a), (b) and (c) respectively.The left halves of each image are presented at a time of ≈ 250 s, and the right halves at ≈ 480s .The averages are plotted up to a radial coordinate of 15, 000 km, which includes the entire oxygen and neon shells, and the innermost part of the carbon shell.

Figure 6 .
Figure 6.Profiles of silicon, oxygen and neon mass fractions as a function of enclosed mass at a range of time snapshots.The top row is for the MHD case, and the bottom row for the purely hydrodynamic model.The RMS averaged radial kinetic energy at ≈476s is plotted in black for both MHD and hydrodynamic models.The grey bands represent the locations of the convective oxygen, neon and carbon shells when the simulation begins

Figure 7 .Figure 8 .
Figure 7. Profiles of the nuclear energy generation rate at 110 s (dotted), 220 s (dashed) and 450 s (solid) for both the MHD (black) and hydrodynamic (red) simulations.The energy generation rates are presented both as a function of mass in panel (a) and as a function of radius in panel (b).
(a), the meridional velocity plotted is an average over  (zonal average) and time of |v  + v  |, while Figure 10 (b) plots the same averages of the angular velocity, Ω =   /( sin ), where  is the radius.Both figures plot the cube root of their original values to retain the direction of the flow, but reduce the dynamic range.The positive (red) meridional velocity represents clockwise motion, with the negative (blue) flows showing counter-clockwise flows (viewed from the North pole).

Figure 9 .
Figure 9. Profiles of angular velocity (top row) and angular momentum (bottom row) as a function of enclosed mass for different times during the simulation.Plots on the left column are for the MHD case, and the right column for the purely hydrodynamic model.The plots are limited to mass coordinates inside 3 M ⊙ where the most clear dynamical differences occur.The shaded bands represent the locations of the convective oxygen, neon and carbon shells when the simulation begins.

Figure 10 .Figure 11 .
Figure 10.The left halves of panels (a) and (b) depict zonal and time average quantities from 250 s to 400 s, while the right halves are meridional snapshots at 300 s.The cube root of meridional velocity is shown in panel (a) and the rotation rate in panel (b).The positive meridional velocities represent clockwise motion.Dashed-dotted black lines mark the approximate initial radii of the base of the oxygen, neon and carbon shells.
al. 2017; Camisassa & Featherstone 2022).To confirm that our model is in the relevant regime where Coriolis forces are Profiles of the Coriolis force (green) and buoyancy force (blue) at 57 s, 115 s, and 191 s for the purely hydrodynamic model.