Magnetized oscillatory double-diffusive convection

We study the properties of oscillatory double-diffusive convection (ODDC) in the presence of a uniform vertical background magnetic field. ODDC takes place in stellar regions that are unstable according to the Schwarzschild criterion and stable according to the Ledoux criterion (sometimes called semiconvective regions), which are often predicted to reside just outside the core of intermediate-mass main sequence stars. Previous hydrodynamic studies of ODDC have shown that the basic instability saturates into a state of weak wave-like convection, but that a secondary instability can sometimes transform it into a state of layered convection, where layers then rapidly merge and grow until the entire region is fully convective. We find that magnetized ODDC has very similar properties overall, with some important quantitative differences. A linear stability analysis reveals that the fastest-growing modes are unaffected by the field, but that other modes are. Numerically, the magnetic field is seen to influence the saturation of the basic instability, overall reducing the turbulent fluxes of temperature and composition. This in turn affects layer formation, usually delaying it, and occasionally suppressing it entirely for sufficiently strong fields. Further work will be needed, however, to determine the field strength above which layer formation is actually suppressed in stars. Potential observational implications are briefly discussed.


INTRODUCTION
Oscillatory double-diffusive convection (ODDC hereafter) takes place in regions of stars or planets that have a destabilizing temperature stratification and a stabilizing composition stratification. These would be identified as convectively unstable according to the Schwarzschild criterion, but convectively stable according to the Ledoux criterion 1 . In stellar evolution calculations that use the Ledouxcriterion for convection, regions with this type of stratification most commonly appear just above the convective core of intermediate-mass and high-mass stars. ODDC can in that case enhance the transport of chemical species between the envelope and the core, which provides additional fuel to the nuclear reactions and prolongs the Main Sequence phase (see, e.g. Moore & Garaud 2016, and references therein).
In its basic form, ODDC is a small-scale, wave-like type of thermal convection. It was identified to be a double-diffusive instability by Kato (1966), and was hence named ODDC by Spiegel (1969). A key question is how much mixing ODDC can cause. The answer turns out to be quite complicated, because it depends on whether the instability remains in its basic small-scale form, or whether it transforms into a state of large-scale layered convection, where mixing is much more efficient. Layered convection had been postulated to exist in giant planets and stars by Stevenson (1985) and Spruit (1992) based on its presence in geophysical analogs of ODDC (e.g. in the polar oceans and volcanic lakes on Earth, see the review by Radko 2013). But it was first clearly identified in direct numerical simulations (DNS hereafter) of ODDC at astrophysically relevant input parameters only recently by Rosenblum et al. (2011).
In that work, Rosenblum et al. (2011) demonstrated that layering (i.e. transition to layered convection) can spontaneously occur in ODDC, and that it is caused by a mean-field instability first discovered by Radko (2003) in a related oceanographic context. Mean-field instabilities are defined mathematically as instabilities of the Reynolds-averaged equations, rather than of the original Navier-Stokes equations. Physically speaking, they can be viewed as secondary instabilities that only develop after a particular primary instability has saturated into a state of homogeneous turbulence. In addition, mean-field instabilities have characteristic lengthscales that are much larger than the typical turbulent eddy scale. Understanding their development requires knowledge (or modeling) of relevant turbulent fluxes such as the heat and composition fluxes (and others as needed) as functions of the large-scale properties of the fluid such as the background temperature gradient, the background composition gradient, and the mean local shear for instance. The specific mean-field instability discovered by Radko (2003) is called the γ-instability, and is triggered by a positive feedback loop between the heat and composition fluxes (whose ratio is traditionally called γ), and the background stratification (see reviews by, e.g. Radko 2013;Garaud 2018).
In subsequent investigations, Mirouh et al. (2012) ran and analyzed a wide range of DNS, to model the turbulent heat and composition fluxes associated with the basic weak form of ODDC. With that information, they were able to construct a quantitative model for the mean-field γ-instability that can be used to predict when layering is expected. A salient result of their analysis is that layering is always expected in ODDC-unstable regions located at the edge of the convective cores in intermediate-mass stars, suggesting that mixing in these regions is very efficient indeed.
With this in mind, Wood et al. (2013) used DNS of layered convection to quantify transport in that regime. They found that both heat and composition fluxes increase with the mean layer height in a way that is strongly reminiscent of (but not entirely identical to) Rayleigh-Bénard convection (i.e. thermal convection between plane parallel plates). But they also found that the layers have a propensity to merge with one another, in such a way that the ODDC-unstable region always eventually evolves into one that is fully convective. This suggests that the ODDC-unstable region is only expected to exist for a relatively short period of the star's life, and is rapidly replaced by a standard convection zone.
Using their results, Moore & Garaud (2016) modeled the evolution of ODDC-unstable regions in intermediate mass stars using MESA (Paxton et al. 2011), and found that mixing is so efficient that these regions can be approximated as being fully convective instead for simplicity. In other words, ignoring ODDC altogether and using the Schwarzschild criterion to identify convectively mixed regions satisfactorily predicts the star's evolution. A similar conclusion was recently reached by Anders et al. (2022), albeit for different reasons.
However, it is important to note that almost all results on ODDC so far have been obtained in the hydrodynamic limit (although see Hughes & Brummell 2021, for a related magnetized instability), but core convection in rotating stars can drive a substantial dynamo (Augustson et al. 2016). The magnetic field thus generated can diffuse out of the core and would likely influence the initial development of the ODDC instability, and/or its subsequent transition into layered convection.
In this work, we are therefore interested in quantifying the effects of magnetic fields on ODDC. For simplicity, we begin by looking at a very simple model setup in which an ODDC-unstable region is threaded by a uniform, vertical magnetic field. Section 2 describes the governing equations for the problem. Section 3 presents a linear stability analysis of ODDC in the presence of this field, demonstrating that the fastest-growing modes of instability are unaffected. Section 4 then presents nonlinear DNS of the system, where in this case the magnetic field is seen to slow down, and sometimes even suppress, the transition to layered convection. We interpret our findings in the light of Radko's γ-instability theory in Section 5, and conclude in Section 6 with a discussion of the impacts of our findings on stellar modeling.

MATHEMATICAL MODEL AND GOVERNING EQUATIONS
We consider a Cartesian domain in a small stellar region with dimensions (L x , L y , L z ) where gravity defines the vertical direction e z . We assume that the vertical extent of the domain is smaller than a pressure scale height, which allows us to use the Boussinesq approximation (Spiegel & Veronis 1960). Within this approximation, we assume that the temperature T and composition C can be decomposed as a linear background plus perturbations, namely where x, y, and z are the spatial coordinates and t is time, T m and C m are constants that represent the mean values of temperature and composition in the region considered, and dT 0 /dz and dC 0 /dz are negative constant gradients of temperature and composition, with L z |dT 0 /dz| T m and L z |dC 0 /dz| C m . Consistent with the use of the Boussinesq approximation, we assume that the density perturbationsρ satisfy a linear equation of state: where ρ m is the mean density of the domain, and α = − 1 ρm ∂ρ ∂T and β = 1 ρm ∂ρ ∂C are constants of thermal expansion and compositional contraction respectively.
With these assumptions, and neglecting magnetic buoyancy effects, the governing equations are: where u = (u x , u y , u z ) is the velocity of the fluid, B = (B x , B y , B z ) is the magnetic field, g is gravity, µ 0 is the vaccuum permeability, andp is the pressure perturbation away from hydrostatic equilibrium. The diffusion coefficients, namely the kinematic viscosity ν, the thermal diffusivity κ T , the compositional diffusivity κ C and the magnetic diffusivity η are all assumed to be constant. Finally, dT ad /dz = −g/c p is the vertical adiabatic temperature gradient, where c p is the specific heat at constant pressure of the fluid. To be considered Schwarzschild-unstable and Ledoux-stable the region satisfies These equations must be complemented by boundary conditions. Two options are possible: to run "global" simulations in bounded domains with boundary conditions enforced on the dynamical fields (e.g. Anders et al. 2022), or to run "local" simulations where background gradients are fixed and periodic boundary conditions are enforced for the perturbed fieldsT ,C, u, and B (e.g. Rosenblum et al. 2011;Mirouh et al. 2012;Wood et al. 2013). The regions in stellar interiors where ODDC might occur are far from boundaries; thus, to minimize their effect on our results, we perform local simulations with periodic boundary conditions. This numerical setup is meaningful and valid as long as two conditions are met. First, the domain size must be larger than the dominant eddy length scales in the saturated state to ensure that the boundaries are not unphysically influencing the turbulence. Second, we note that the periodic model implicitly maintains large-scale temperature and composition jumps across the height of the domain, letting the turbulence develop accordingly. This is only a good approximation to the true physical problem if the timescale over which the turbulence would in reality modify (and reduce) both jumps is long compared with the timescale for the development and nonlinear saturation of the instability. This was shown to be the case for the related fingering instability by Zemskova et al. (2014), and we assume that it holds here as well for ODDC. This assumption is very likely justified in the case where ODDC remains in its weak, small-scale, wave-like form, but could be invalidated in the layered case (where convective fluxes transport both heat and composition very rapidly). Studying the transition to layered convection in a different model setup is the subject of a forthcoming paper.
We non-dimensionalize the system using the expected width of diffusive structures d given by (Stern 1960) [ where N T is the local buoyancy frequency associated only with the temperature field. The corresponding choices for the other units are: Finally, we assume the existence of a constant background field of amplitude B 0 , so and we use this amplitude as the unit for the magnetic field strength. With this, the non-dimensional system of governing equations is: where hatted quantities, as well as time and space, are now non-dimensional. Several parameters appear, namely the so-called inverse density ratio which measures the stratification of the region, three diffusivity ratios which only depend on the properties of the fluid, and a parameter controlling the magnetic field strength, which is the square of the ratio of the Alfvén velocity to the anticipated ODDC velocity κ T /d.
Note that in the hydrodynamic case, a region is linearly unstable to ODDC provided linearly stable if R −1 0 > R −1 c , and unstable to thermosolutal convection if R −1 0 < 1 (Baines & Gill 1969). In stellar interiors, the diffusivity ratios are always quite small (cf. Garaud et al. 2015), so R −1 c is very large. Also note that an estimate for H B near the core of intermediate-mass stars is showing that a magnetic field with amplitude of about 10 4 G can in principle have a substantial impact on ODDC.
The system of equations (9)-(13) is identical to the one used by Harrington & Garaud (2019) to study magnetized fingering convection, a related double-diffusive instability, except for the negative signs in front of the termsû z and R −1 0û z in the temperature and composition equations respectively (see their equations 14-18).

LINEAR STABILITY ANALYSIS
We analyze the stability of the system (9)-(13) to infinitesimal perturbations following standard procedures. First, we explicitly writeB = e z +b, whereb are magnetic field perturbations. The background state around which we linearize is a steady state withû =b =T =Ĉ =p = 0. We assume the perturbations to be sufficiently small so that nonlinearities can be ignored. Finally, we seek solutions of the formq = q exp(ik x x + ik y y + ik z z +λt) whereλ is the growth rate of an unstable mode with wavenumberk = (k x ,k y ,k z ). With these steps and assumed ansatz, the governing equations become: wherek 2 =k 2 x +k 2 y +k 2 z =k 2 h +k 2 z andk h is a horizontal wavenumber. The linearized induction equation can be simplified as while the temperature and composition equations can be combined into Substituting (24) and (25) into the linearized momentum equation, taking the dot product of the result withk, and using incompressibility, yields an expression for p . After substituting it back into the z-component of the momentum equation, we finally obtain a quartic equation forλ: (λ +k 2 )(λ + τk 2 ) λ + Prk 2 + H Bk This quartic has several notable properties. We see that setting H B = 0 recovers the standard cubic associated with hydrodynamic ODDC, see, e.g. Rosenblum et al. (2011). In addition, the only difference with the quartic obtained by Harrington & Garaud (2019) describing the linear stability of magnetized fingering convection is the sign of the right-hand side of this equation. As such, many of their results hold here as well. For instance, the magnetic field has no effect on the stability of modes withk z = 0 (i.e. modes that are invariant along both gravity and the background magnetic field). This is expected, since magnetic fields only interact with flows perpendicular to the field lines. However, since thek z = 0 modes are also the most rapidly growing modes when H B = 0, we conclude that the presence of a vertical magnetic field does not change the stability of the most rapidly growing modes. 2 Writing (26) in the formλ 4 + α 3λ 3 + α 2λ 2 + α 1λ + α 0 = 0 gives the following expressions for the coefficients α j : We solve this quartic forλ for given input parameters (R −1 0 , H B , Pr, τ and D B ), and given values of the wavenumberk. Out of the four complex roots, we select the one with the maximum real part. From here on, we only ever discuss the properties of this particular root. We present the results of the linear stability analysis in Figure 1. Each row corresponds to a different magnetic field strength, namely H B = 0, 0.01, 1, 100, from top to bottom. From left to right in each column, we present (a) the real part ofλ as a function ofk x andk z (b) the imaginary part ofλ as a function ofk x andk z (selecting for simplicity the root with positive imaginary part), (c) the real part ofλ as a function ofk z , atk x = 0.417, and (d) the imaginary part ofλ as a function ofk z , atk x = 0.417. Note that the color scales on the first and second columns are quite different, and thatk x = 0.417 is the fastest growing mode atk z = 0 for the given parameters. In the third column, we compare Re(λ) to a purely diffusive solution of the form −c(k 2 x +k 2 z ), where c is the minimum of Pr, τ , and D B . In the last column, we compare |Im(λ)| to the oscillation frequency of a corresponding gravity wave, namely and to the frequency of a corresponding Alfvén waves, namely In general, and consistent with the discussion above, we see in the left-most column that the magnetic field partially or fully stabilizes every mode except those withk z = 0. More specifically, we see that at a given (non-zero) value ofk x , increasing the magnetic field strength reduces the range ink z of unstable modes, and reduces the growth rates of all modes exceptk z = 0.
In the second column, we see that as the value of H B increases, the oscillation frequency of the unstable modes becomes gradually influenced by the presence of the field. Indeed, in the first two rows of the last column, the field is zero or weak, and we find that the oscillation frequency of thê k z = 0 modes are related to (but not equal to) the oscillation frequency of a pure gravity wave ω g . For larger values of H B (third and fourth row, where H B ≥ 1), the oscillation frequency of the modes approaches the corresponding Alfvén frequency ω A instead, implying that the magnetic field now dominates the dynamics of the instability.
A cursory examination of the decay rate of the decaying modes in the third column reveals that the latter has a parabolic dependence onk z for large enoughk z . This can be explained by noting that for large wavenumber, the quartic (26) asymptotically tends to: This can be factored as (λ − Prk 2 )(λ − τk 2 )(λ − D Bk 2 )(λ −k 2 ) = 0, which has rootsλ = −ck 2 , where c can be 1, Pr, τ or D B . The slowest decaying mode is therefore obtained by taking c = min(1, P r, τ, D B ). The parabolaλ = −ck 2 is shown in green in the third column, and we see that the actual solution of the quartic (26) (blue curve) tends to the diffusive solution (green curve) at large values ofk z . In stars, we expect c = τ , as we almost always have the ordering of parameters τ < Pr < D B < 1.
To summarize, we have found that a vertical magnetic field has no effect on the growth rate and nature of the fastest-growing, z−invariant mode of instability in this model setup. It does, however, affect both the growth rate and oscillation frequency of inclined modes, that presumably play a role in the saturation of the instability.

NUMERICAL SIMULATIONS
To study the behavior of magnetized ODDC beyond the initial growth of the instability, and determine whether magnetic fields affect the processes by which convective layers form (see Section 1), we now turn to DNS. Equations (9)-(13) are evolved with time using the PADDIM code developed by Harrington & Garaud (2019). PADDIM itself is based on the original pseudo-spectral PADDI code (Traxler et al. 2011b;Stellmach et al. 2011) used in our previous non-magnetic studies (Rosenblum et al. 2011;Mirouh et al. 2012;Wood et al. 2013).
We select a parameter regime where layers are known to form in the absence of magnetic fields, namely Pr = τ = 0.3, and R −1 0 = 1.2 (Mirouh et al. 2012), and run three simulations with increasing background magnetic field strength: H B = 0 (hydrodynamic reference case), H B = 0.03, and H B = 0.1. The magnetic diffusion coefficient is fixed and equal to D B = 0.3. Initial conditions in each simulation are:B x =B y = 0,B z = 1 (for the magnetic simulations only), and both temperature and composition fields are initialized with small-amplitude white noise. In each case, the domain size used has size 100d × 100d × 300d, and the resolution is 192 × 192 × 576 equivalent grid points.
H B increases from 0 to 100 from the top row to the bottom row. The first column shows Re(λ) as a function ofk z andk x . The second column similarly shows |Im(λ)|. The third and fourth columns show Re(λ) and |Im(λ)|, respectively, at fixedk x = 0.417 (the fastest growing mode atk z = 0). In the third column, Re(λ) (blue solid line) is compared with the diffusive solution −ck 2 (green dashed line), see the main text for details. In the fourth column, Im(λ) (blue solid line) is compared with the oscillation frequencies of corresponding gravity waves ω g (brown dashed curve) and corresponding Alfvén waves ω A (green dashed curve).
Typical snapshots from the simulation with H B = 0.03 are shown in Figure 2, illustrating the evolution of the system. Note that this would correspond to a weak-field case according to linear theory. The left column shows the composition fieldĈ, the middle column shows the vertical velocity fieldû z , and the right column shows the vertical magnetic fieldB z . The top row shows the early evolution of the instability and the emergence of the oscillating vertically-invariant modes. These modes quickly saturate in the second row into a state that is dominated by weakly nonlinear waves, which is what is commonly called homogeneous oscillatory double-diffusive convection. Later in the simulation (third row) a stack of convective layers emerge, and then successively merge (fourth row), until a single layer remains (not shown). The amplitude of the flow velocity and of the magnetic field perturbations is much larger in the layered phase than in the homogeneous phase, consistent with the findings of Rosenblum et al. (2011) and Wood et al. (2013). Layers are much more visible in the composition field than in the vertical velocity and magnetic field. This is because the density contrast across the interface is weak at the parameters selected, and is therefore not a strong barrier to vertical fluid motions. Instead, the interface is very turbulent and regularly pierced by strong upflows or downflows (cf. Wood et al. 2013).
The left panel shows the early-time behavior, while the right panel shows the entire evolution. The horizontal black lines on the right panel show the kinetic energy density of the non-magnetic simulation measured, from bottom to top, in the four-layer phase, the three-layer phase, and the two-layer phase, to help compare it to that of the magnetized simulations.
A more quantitative way of looking at the evolution of the flow, which allows us to compare the outcome of simulations run at different magnetic field strengths, is presented in Figure 3. This figure shows the temporal evolution of the kinetic energy density KE and magnetic energy density M E, defined as where the angled brackets denote a volume average. The left-side panel focuses on the early-time behavior of the energies, while the right-side panel shows the entire time evolution.
For H B = 0, the simulation behaves as described by Mirouh et al. (2012). The primary ODDC instability develops rapidly, as can be seen by the exponential growth of the kinetic energy, and saturates at a low amplitude around t 250. The kinetic energy grows again between t 250 and t 500, as a result of the excitation of larger-scale gravity waves that then saturate at a fairly constant amplitude between t 500 and t 1000. During that time, a stack of seven layers forms as the result of an underlying mean-field γ-instability (see next section for details). This is accompanied by a step-like increase in the kinetic energy (around t 1000), which corresponds to the onset of overturning convection in each layer. The seven layers very quickly merge into five, which then more slowly merge into four, three and then two layers around times t 1200, t 1500 and t 2000, respectively. Had the simulation been carried out for longer, we anticipate that the last two layers would also eventually merge. Each successive merger is accompanied by another strong step-wise increase in the kinetic energy. This is consistent with the findings of Wood et al. (2013), who noted that the efficiency of layered convection increases significantly with the average layer height.
We find that the presence of a vertical magnetic field with H B = 0.03 and 0.1, for the chosen parameters (P r = τ = D B = 0.3, R −1 0 = 1.2) has substantial quantitative effects on both the initial saturated state of the instability and on the development and evolution of layered convection, despite being relatively weak from the perspective of linear theory. More specifically, we see that the magnetic field can reduce the kinetic energy in the flow prior to layer formation substantially (e.g. comparing energy levels in the phase prior to layer formation between the H B = 0 case, where KE 2.4 and the H B = 0.1 case, where KE 1.8). Crucially, we also find that layers take much longer to appear in the magnetic cases than for H B = 0, and that the first layered configuration has fewer layers (see Section 5.4). For H B = 0.03 for instance, layered convection starts around t 4300 with four layers which rapidly merge into three layers around t 4500. By t 5000, they have merged down to two layers and by the end of the simulation are in the process of merging into a single layer. In the case of H B = 0.1, the system stays in a state of weak wave-like convection for a very long time before layers start to emerge. Layered convection does not start until around t 8300, and again begins in a four-layer phase. It is interesting to note that the rate at which the layers merge, however, appears to be relatively independent of H B .
Finally, we also see that the kinetic energy in the layered phase depends on the magnetic field strength. Indeed, the black horizontal lines in Figure 3 indicate, from top to bottom, the average kinetic energy of the two, three and four layered phases in the non-magnetic (H B = 0) case, and help guide the eye for comparison with the magnetic cases. We see that the kinetic energy in each of the layered phases are lower for H B = 0.03 and H B = 0.1 than in the non-magnetic case, and that the effect becomes stronger as the magnetic field increases. This suggests that a strong enough magnetic field can substantially reduce the efficiency of layered convection. Meanwhile, we also see that the magnetic energy increases not only with H B (as expected), but also with each merger during the layered phase, demonstrating that some of the kinetic energy from the convective motions is being converted into magnetic energy.
In summary, our DNS demonstrate that the presence of a vertical magnetic field can both slow down (or suppress, in some cases shown in the next sections) the formation of layers, and reduce the kinetic energy of the flow in layered convection. Given the likely importance of layer formation and mergers in ODDC near the core of intermediate-mass stars (see Section 1 and Moore & Garaud 2016), it is therefore crucial to establish whether these statements still hold at stellar parameter values (which would have significant consequences for stellar evolution). To do so, we first need to determine how the layer-forming γ-instability is affected by the magnetic field, both qualitatively and quantitatively. The next section investigates this question in detail. The reader merely interested in the answer may skip to Section 6.

Theory
As mentioned in Section 1, in the absence of magnetic fields, the spontaneous emergence of convective layers from homogeneous ODDC has been attributed to the so-called γ-instability (Radko 2003;Rosenblum et al. 2011). We now demonstrate that a large-scale uniform magnetic field does not directly influence the γ-instability mechanism, so the latter is still expected to operate as before.
While most standard fluid instabilities are instabilities of the original Navier-Stokes equations, and describe the evolution of perturbations from an initial laminar state, the γ-instability is an instability of the Reynolds-averaged equations (where in this case, the average is simply a horizontal average) which model the evolution of large-scale spatial modulations to a pre-existing turbulent state. It is caused by a positive feedback loop between the vertical turbulent temperature and compositional fluxes (that cause an evolution of the temperature and composition gradients) and the local stratification (that controls these turbulent fluxes), which is schematically illustrated, for instance, in Fig. 3 of Garaud (2018). The physical mechanism driving the γ-instability is quite subtle. While we describe the mathematical details below, readers interested in a more schematic illustration are referred to Fig. 3 of Garaud (2018).
This feedback loop can be studied as follows (Radko 2003;Mirouh et al. 2012). First, we take the horizontal average of (11) and (12), which results in: whereT is the horizontal average of the fluctuationsT , and similarly forC. We have introduced the turbulent temperature and composition fluxes F T =û zT and F C =û zĈ , as well as the total temperature and composition fluxes Equations (32) and (33) describe the first part of the feedback loop, namely how spatial modulations of the fluxes cause a temporal evolution of the temperature and composition stratification. Note how the magnetic field does not appear in their derivation explicitly -in fact, the same equations apply in the hydrodynamic case. However, it does so implicitly since the turbulence, and therefore the turbulent fluxes, are affected by the field. In a strictly homogeneous turbulent system,T C 0, and the temperature and composition gradients are constant and equal to the non-dimensional background values of −1 and −R −1 0 , respectively. In that case the temperature and composition fluxes F tot T and F tot C are independent of height (z) in the domain. This is the state illustrated in the second row of Figure 2.
A crucial ingredient for the second part of the feedback loop in Radko's γ-instability theory is the fact that, for a given fluid (i.e. at fixed P r, τ and D B ) and a given background magnetic field (i.e. fixed H B ), the intensity and properties of double-diffusive turbulence depend on the local temperature and composition gradients only via the local inverse density ratio (Radko 2003), written by Mirouh et al. (2012) for ODDC as This quantity is the ratio of the local density gradient due to composition stratification, to the local density gradient due to temperature stratification, and is equal to R −1 0 in the homogeneous state. To model the γ-instability mathematically, Radko (2003) then introduced two important nondimensional quantities: the thermal Nusselt number, which is the ratio of the total temperature flux to the diffusive temperature flux and the (inverse) flux ratio which is the ratio of the total composition flux to the total temperature flux. The key assumptions discussed above can be expressed mathematically by requiring that both N u T and γ −1 tot should only be functions of R −1 (at fixed P r, τ , D B and H B ). With this assumption, we have If the functions N u T (R −1 ) and γ −1 tot (R −1 ) are known, then equations (32), (33), (36), and (39) form a closed system describing the evolution of large-scale spatial inhomogeneitiesT (z, t) andC(z, t). These equations can then be linearized around the previously introduced homogeneous turbulent state (which hasT =C = 0, and R −1 = R −1 0 ), in the limit whereT andC are very small, to demonstrate that perturbations of the kindT (z, t) ∼ e ikz+Λt andC(z, t) ∼ e ikz+Λt grow or decay exponentially with time, with a growth rate where σ satisfies the quadratic equation and where the coefficients a and b depend on the properties of the homogeneous state as with A detailed derivation of this result is presented in Rosenblum et al. (2011) and Mirouh et al. (2012). Crucially, we see from the steps outlined above that the presence of a magnetic field does not explicitly appear in this theory. However it does so implicitly by influencing the functions N u T (R −1 ) and γ −1 tot (R −1 ), and therefore the coefficients a and b of the quadratic equation (41).
Equation (41) immediately shows that a necessary condition for instability, i.e the existence of solutions with σ > 0, is i.e. that the inverse flux ratio should be a decreasing function of the inverse density ratio (Radko 2003;Mirouh et al. 2012). This γ-instability theory, and more specifically equation (41) was found to correctly predict the growth rate of large-scale perturbations in both fingering convection in the ocean (Radko 2003;Stellmach et al. 2011) and ODDC in stars (Rosenblum et al. 2011;Mirouh et al. 2012). We now investigate whether this remains true for magnetized ODDC.

Comparison with numerical experiments 1: layering vs. no layering
To test the γ-instability theory derived above, we follow the same steps as in Mirouh et al. (2012). The first step consists in measuring the functions N u T (R −1 ) and γ −1 tot (R −1 ) in homogeneous ODDC turbulence at fixed values of P r, τ , D B and H B . To do so, we run and analyze a number of smalldomain DNS, each of which is performed in a triply-periodic cube of size 100d × 100d × 100d with a resolution of 384 × 384 × 384 equivalent grid points. Each simulation is initialized as those of Section 4, with all fields set to zero except forT andĈ which are seeded with small random noise, andB z which is set to 1. The equations are evolved either until the first set of layers form, or until we have acquired enough data to be certain that they do not. As in Mirouh et al. (2012), we use the fact that when the flow is in a statistically stationary state, to compute We then take a time average of N u T and γ −1 tot in the statistically stationary homogeneous ODDC phase.
Note that our determination of the appropriate time interval for this average deviates somewhat from the method used by Mirouh et al. (2012) (see their Appendix B), mostly for reasons of simplicity (see below), but also because their method turns out to be somewhat inconsistent with the assumptions of the γ-instability theory outlined above. In what follows, we use the same steps as Mirouh et al. (2012) to determine the averaging start time, but set the averaging end time to be either when the horizontally-averaged density profile begins to deviate substantially from the background linear profile (when layers form), or at the end of the simulation (when no layers form). Mirouh et al. (2012), by contrast, stopped the averaging either when fully convective layers appear, or when large-scale gravity waves begin to dominate the simulation, whichever happens the soonest. As a result, their estimates for the convective fluxes are systematically higher than ours in very low R −1 0 simulations where layers form early, because their average includes times during which the layers have almost (but not completely) overturned, while we stop before this happens. Meanwhile, their estimates for the convective fluxes are often lower than ours at moderate and high R −1 0 , because the large-scale gravity waves that eventually appear in that limit (see Mirouh et al. 2012;Moll et al. 2016), which they discard but we keep, cause a non-negligible amount of transport.  Table 5, which they had extracted using their method. The difference between the solid and dashed red curves illustrates the fairly substantial impact of choosing a different interval for the time averages, as mentioned above. That impact is, however, smaller than the impact of adding a magnetic field. Indeed, and consistent with the results briefly discussed in Section 4, we see that even a weak magnetic field (e.g. H B = 0.03, green curves) can sometimes substantially reduce the turbulent heat flux, decreasing N u T . This in turn changes the shape of the γ −1 tot (R −1 0 ) curve, which gradually tends towards the diffusive flux ratio γ −1 dif f (R −1 0 ) = τ R −1 0 (solid black lines). The theory developed in the previous section, combined with these results, therefore shows that a large-scale magnetic field has an indirect impact on layer formation, by affecting N u T and γ −1 tot , which changes the values of the coefficients of the quadratic (41), and in turn, the γ-instability growth rate.
A rapid, qualitative way of testing the predictions of the γ-instability theory is to verify that all simulations that were run at parameters for which dγ −1 tot /dR −1 0 < 0 (cf. equation 44) eventually transition into layered convection. This is indeed almost always the case. In Figure 4 (bottom row), simulations that ultimately become layered have an additional open square surrounding the original symbol. We see that the square is present whenever γ −1 tot is a rapidly decreasing function of R −1 0 , as in Mirouh et al. (2012). We also see that no layers ever form when γ −1 tot increases with R −1 0 , consistent with the theory. The only simulations in which the data disagrees with the theory are those for which γ −1 tot decreases very slowly with R −1 0 (near the minimum of the curve), where layering was not observed even though the γ-instability ought to be active. This is possibly because the growth rate of the γ-instability is too low in that limit, and other effects that are not accounted for (such as the presence of large-scale gravity waves in the system) further damp it (see, e.g. Traxler et al. 2011a, for related effects in oceanic fingering convection).
Notwithstanding this minor discrepancy, we now also see that a sufficiently strong magnetic field can, in some cases, completely suppress layer formation. Indeed, consider the case with P r = τ = D B = 0.1. At these parameters, Figure 4 (bottom right panel) shows that layers form at R −1 0 = 1.4 in the non-magnetic case, and in the H B = 0.03 case, but not for H B = 0.1. At the same time it provides a tentative explanation why, namely by moving the minimum of the γ −1 tot (R −1 0 ) curve slightly to the left, which eventually stabilizes the system to the γ-instability at these parameters.

Comparison with numerical experiments 2: growth rate of the layering modes in small domains
Even when layers eventually form, we have seen in Section 4 that a weak magnetic field can significantly delay the onset of layered convection. The same was found to be true in all of the small domain simulations discussed in Section 5.2. The theory presented in Section 5.1 suggests a possible explanation for this delay, namely that the magnetic field reduces the turbulent temperature and composition fluxes, which in turn reduces the growth rate of the γ-instability. We now test this idea more quantitatively.
We begin by comparing the γ-instability theory predictions against data from small-domain simulations. This turns out to be easier than starting with the large-domain simulations of Section 4, because the latter have many modes growing simultaneously (layering modes and sometimes largescale gravity waves). Their presence can obfuscate the dynamics of the γ-instability, as discussed below. We therefore focus on four available small-domain layer-forming simulations, whose parameters are presented in Table 1. The table also shows the extracted values of N u T , γ −1 tot , A N u and A γ for these simulations. The derivative terms A N u and A γ are computed using second-order finite differences using simulations at values of R −1 0 on both sides of the target one. With this information, we can then evaluate the quadratic coefficients a and b using (42), and solve (41) for σ (also shown in Table 1). Finally, to compute the growth rate Λ n of a particular layering "mode" with n layers, we use the fact that We immediately see from Table 1 that σ decreases substantially when H B increases from 0.03 to 0.1 with all other parameters fixed. This confirms our above hypothesis that the magnetic field can reduce the growth rate of layers.
For a more direct test of the γ-instability theory against the data, we now compare the growth rates of the layering modes observed in the simulations to those predicted by the theory. We extract the time-dependent amplitude of these modes from the DNS by computing the Fourier expansion of  Table 1. Values of N u T , γ −1 tot , A N u , A γ extracted from DNS at R −1 0 = 1.15, P r = τ = D B = 0.3 (top two lines) and R −1 0 = 1.3, P r = τ = D B = 0.1 (bottom two lines). The derivatives A N u and A γ are computed using finite differencing with data at neighboring values of R −1 0 .
the horizontally-averaged density perturbations, namelȳ The quantity |ρ n | 2 (t), called the density spectral power hereafter, is plotted in Figure 5 for the four simulations presented in Table 1, for modes leading to n = {1, 2} layers in each case. Modes with n > 2 are not usually found to grow in small domain simulations. The γ-instability theory predicts that |ρ n | 2 (t) ∝ e 2Λnt , with Λ n given by (48) for a given Fourier mode (equivalently, number of layers) n. We compare these predictions (colored solid lines) to the data for each mode in Figure 5. We see that in all cases, the model is appropriate for the n = 2 mode at very early times, but overestimates its growth rate by a factor of about two at later times (after t ∼ 600). Predictions made with half the growth rate (dashed lines) appear to fit the data better then. The situation is not as clear for the n = 1 mode (in some cases, the above statements hold, and in some others, they do not), which is perhaps not too surprising because the latter is growing intrinsically slowly, in a fairly turbulent environment.
The fact that some layering modes grow slower than expected at later times can already be seen in some of the hydrodynamic simulations of Mirouh et al. (2012) but was not discussed in that paper. However, we now see that this is a relatively systematic effect. Inspection of the total density profile before, during, and after the time where the mode growth rate starts decreasing reveals that this corresponds to the point where the total density profile is no longer linear. At this point, the linearization procedure that leads to the derivation of the quadratic growth rate equation (41) is no longer valid, and it is therefore not surprising to see that the model no longer fits the data at the quantitative level. We do see, however, that the mode continues to grow, albeit at a smaller rate.

Comparison with numerical experiments 3: growth rate of the layering modes in large domains
To conclude this analysis, we now compare the predictions of the γ-instability theory with the layering mode data for the three large-domain simulations presented in Section 4, which have P r = τ = D B = 0.3, R −1 0 = 1.2. The results are presented in Figure 6. Each row of Figure 6 corresponds to a different magnetic field strength, from top to bottom H B = 0, 0.03 and 0.1. The density spectral power of the first n L layering modes are shown in each case, where n L is the number of layers of the mode that dominates at early times. As discussed in Section 4, n L = 7 for H B = 0 (pink curve), while n L = 4 for H B = 0.03 and H B = 0.1 (red curve). Modes with n > n L (not shown) are not found to grow in the simulation. It is worth noting that it is not clear why n L is smaller in the magnetized cases than in the hydrodynamic case, and whether this is a systematic result, or a coincidence. We discuss this issue in Section 6.  Table 1 (see text for detail). In each panel, the colored curve show |ρ n (t)| 2 for n = 2 (blue) and n = 1 (green), and the solid line of the same color shows the predicted exponential growth of the mode according to the γ-instability theory. The dashed line of the same color shows the same with half the growth rate, which appears to fit the data much better overall. Table 2 shows N u T , γ −1 tot , A N u , A γ and σ extracted from the small domain simulations, this time for the parameter values used here. All the quantities are computed as in Table 1. The predicted growth rates, calculated as in the previous Section for each dominant mode (i.e. a mode whose amplitude is much larger than the others), are shown in Figure 6 as solid lines. Dashed lines show how these modes would grow at half the predicted growth rates. As before, we find that most modes grow at the predicted rate at early times (before t ∼ 600), but then later continue at about half the predicted rate. The strongest field case (H B = 0.1) is a little different, however, as discussed below.
For the non-magnetic case, we see that the 7-layer mode clearly dominates at early times until t 1000. It grows roughly at the predicted rate until t 600, at which point it begins to grow at about half the predicted rate. Around t 900, we saw in Figure 3 that the kinetic energy of the flow increases substantially. Inspection of the total density profile at that time (see right-hand panel, solid black line) shows that it is no longer linear, but instead, clearly exhibits the presence of a 7-layer  Table 2. Values of N u T , γ −1 tot , A N u , A γ extracted from DNS at R −1 0 = 1.2, P r = τ = D B = 0.3. The derivatives A N u and A γ are computed using finite differencing with data at nearby values of R −1 0 .
mode, with some of the layers already being fully convective (i.e. with a density that increases with height). It is therefore not surprising to see that the mode stops growing shortly after t = 900. Note that the amplitude above which a single mode with n layers causes an inversion in the total density profileρ tot (z, t) = (1 − R −1 0 )z +ρ(z, t) was given by Rosenblum et al. (2011) to be This threshold, computed for n = n L , is shown as a horizontal black line in each panel of Figure 6. We see that in the non-magnetic calculation, the time at which the kinetic energy of the flow begins to increase (t = 900, vertical black line) corresponds to the time at which the 7-layer mode amplitude reaches |ρ n | crit . Beyond that point, and as discussed in Section 4, convective layers appear and rapidly begin to merge. We see, accordingly, that the dominant mode changes with time. A second density profile is shown at t = 2200 (black dashed line), showing two convective layers separated by thin interfaces.
In the weak magnetic field case (middle row) the situation is overall similar -a dominant layering mode grows from the γ-instability, and layers eventually appear. This mode has fewer layers (n L = 4) than in the non-magnetic case, but its growth rate continues to be reasonably well predicted using half the theoretical growth rate after t = 600. By contrast with the non-magnetic case, however, it does not remain dominant until it overturns, but instead, appears to somehow aid the growth of larger-scale layering modes (n = 3 and n = 2) around t = 1800. Inspection of the density profile at that time reveals the presence of a single, shallow convective layer in the lower half of the domain, that may have been created earlier than expected by a large-scale gravity wave breaking in a region of the domain whose stratification was already weakened by the presence of the layering modes. This convective layer remains in place for the rest of the simulation, and therefore couples the n = 2, n = 3 and n = 4 layering modes. All three modes then continue to grow at a slower rate, until the convective layers are fully established around t = 3800 (see right-side panel with two layers, one centered around z = 80 and another centered around z = 0.) Finally, the situation for the stronger field case H B = 0.1 is once again a little different. The n L = 4 mode is dominant at early times, but seems to couple with both n = 3 and n = 2 modes around t = 2000. By contrast with the H B = 0.03 simulation, however, there is no evidence for a convective layer at that point, and the total density profile remains close to being linear until t 6000. All three modes nevertheless continue to grow at a rate that is consistent with half the predicted growth rate of the n = 2 mode. Overturning convection is again triggered a little earlier than the time at which the amplitude of any of the modes reaches the critical threshold |ρ n | crit .  Figure 6. Left: Density spectral power |ρ n | 2 as a function of time for various layering modes for the simulations discussed in Section 4 with P r = τ = D B = 0.3, and R −1 0 = 1.2. The magnetic field increases from top to bottom. The matching colored solid lines show the predicted growth of the layering mode according to the γ-instability theory, and the dashed colored lines shows the same with half that growth rate. The horizontal black line shows the overturning threshold for the mode with n L layers, with n L = 7 in the non-magnetic case (top), and n L = 4 in the other cases (middle and bottom). The solid and dashed vertical lines mark the times at which density profiles are shown (right). The solid vertical line also marks the time at which the kinetic energy first increases in Figure 3. Right: Horizontally averaged total density profiles shown at selected times in the same simulations.
In the previous section, we have demonstrated that the γ-instability theory of Radko (2003) continues to be a good model for layer formation in magnetized ODDC, at least qualitatively, correctly predicting whether layers form or not. It can sometimes overestimate the growth rate of layer-forming modes by a factor of order unity, but this is not too much of a concern. Indeed, being driven by turbulent mixing processes, the timescale for layer formation is much shorter than any evolutionary timescale (regardless of the magnetic field strength). This suggests that layers would appear almost instantaneously from the perspective of stellar evolution whenever the γ-instability is excited. We have also seen that layers rapidly merge once they form, ultimately leading (in stars) to a fully-mixed convective layer. This merger process does not appear to be affected by the magnetic field, at least for the parameters achievable in the DNS.
These results, when combined, suggest that there are two possible outcomes for ODDC-unstable regions in stars: either the γ-instability is excited (i.e. has a positive growth rate), in which case the region rapidly becomes fully convective, or the γ-instability is not excited (i.e. has a negative growth rate), in which case the region remains in a state of weakly turbulent ODDC. This conclusion is not new (Moore & Garaud 2016), but we now confirm that it remains true for magnetized ODDC.
In Section 5.1, we recalled that a sufficient criterion for the γ-instability to occur is that the inverse flux ratio γ −1 tot (i.e. the ratio of the total composition flux to the total temperature flux caused by the basic ODDC instability) should be a decreasing function of the inverse density ratio R −1 0 . This criterion still applies in magnetized ODDC. But in Section 5.3 we also found that even a relatively weak magnetic field can change the shape of the γ −1 tot (R −1 0 ) curve and move the position of its minimum, thereby shrinking the region of parameter space unstable to layering. As a result, systems whose stratification is unstable to the γ-instability in the absence of magnetic fields can be stabilized when the field exceeds a certain threshold.
Unfortunately, the threshold magnetic field that was relevant to our DNS is not directly applicable to stars. This is because stellar fluids generally have much smaller diffusivity ratios P r, τ , and D B than what can be achieved numerically, which means that we cannot directly compute the γ −1 tot (R −1 0 ) curve relevant for stars (at least, not with currently available supercomputing power). Future theoretical work will therefore need to identify the mechanism responsible for the saturation of magnetized ODDC, to better predict the dependence of the turbulent fluxes on all input parameters, especially the field strength H B . With that information, we might then be able to predict the shape of γ −1 tot (R −1 0 ) at stellar values of P r, τ and D B , for varying H B (similar to the model of Mirouh et al. 2012, for the non-magnetic case). The position of the minimum of that curve, and its dependence on the magnetic field strength, then determines whether a particular ODDC-unstable region in the star, with a given stratification characterized by R −1 0 , is unstable to layering. In the event the new model reveals that ODDC in stars can indeed be stabilized against the γ-instability by a sufficiently strong (but still realistic) magnetic field, this could lead to interesting observational predictions. Indeed, Moore & Garaud (2016) showed the ODDC-unstable region surrounding the Ledoux-sized core of intermediate-mass stars rapidly becomes fully convective (as a result of the γ-instability) in the non-magnetic case. These stars therefore have a larger-thanexpected convective core whose size is appropriately computed using the Schwarzschild criterion instead. Intermediate-mass stars with a sufficiently strong magnetic field would, by contrast, remain in a state that has a smaller convective core, surrounded by a region of weak ODDC. Asteroseismic observations of Ledoux-sized cores, should they arise, would therefore point to the presence of a strong magnetic field.
Finally, it is worth noting that other authors have also argued that the boundaries of convective cores are best-described by the Schwarzschild criterion on grounds that are entirely distinct from the existence of the ODDC instability. As demonstrated by Anders et al. (2022), convective entrainment gradually pushes the location of the boundary predicted by the Ledoux criterion outwards until it agrees with the Schwarzschild criterion, on a timescale that is fast compared to stellar evolutionary timescales (see also Gabriel et al. 2014;Paxton et al. 2018Paxton et al. , 2019, where issues stemming from miscalculations of convective boundaries are discussed, and where the appropriateness of the Schwarzschild criterion is also argued). Thus, there are several distinct physical arguments for using the Schwarzschild criterion over the Ledoux criterion in stellar evolution models, in the absence of magnetic fields. In this paper, we have proposed that there may be magnetic fields of sufficient strength to stop the weak form of ODDC from spontaneously evolving into standard convection in these regions that are Schwarzschild-unstable but Ledoux-stable. This does not address whether the other arguments for using the Schwarzschild criterion over the Ledoux criterion still hold in MHDto address this question, further studies on convective entrainment in MHD are necessary.
A. S, A. F. and P. G. acknowledge funding by NSF AST 1908338. A. S. acknowledges funding from the Other Worlds Laboratory at UC Santa Cruz. Most simulations were run on the Lux supercomputer at UC Santa Cruz, funded by the NSF MRI grant AST-1828315. This work also used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. We thank Evan Anders and Adam Jermyn for useful discussions.