Magnetized fingering convection in stars

TL;DR

This study tests and extends a model of magnetized fingering convection in stars, revealing its limitations and improving its accuracy for stellar conditions, which impacts understanding of stellar mixing processes.

Contribution

The paper provides an extended model for magnetized fingering convection that better predicts chemical mixing rates across a wider parameter space in stellar environments.

Findings

The original model does not always accurately predict mixing rates when magnetic diffusivity is high.

An extended model improves prediction accuracy, especially near marginal stability conditions.

Implications for stellar structure and evolution are discussed based on the new model.

Abstract

Fingering convection (also known as thermohaline convection) is a process that drives the vertical transport of chemical elements in regions of stellar radiative zones where the mean molecular weight increases with radius. Recently, Harrington & Garaud (2019) used three-dimensional direct numerical simulations to show that a vertical magnetic field can dramatically enhance the rate of chemical mixing by fingering convection. Furthermore, they proposed a so-called "parasitic saturation" theory to model this process. Here, we test their model over a broad range of parameter space, using a suite of direct numerical simulations of magnetized fingering convection varying the magnetic Prandtl number, magnetic field strength, and composition gradient. We find that the rate of chemical mixing measured in the simulations is not always predicted accurately by their existing model, in particular when the magnetic diffusivity is large. We then present an extension of the Harrington & Garaud (2019) model which resolves this issue. When applied to stellar parameters, it recovers the results of Harrington & Garaud (2019) except in the limit where fingering convection becomes marginally stable, where the new model is preferred. We discuss the implications of our findings for stellar structure and evolution.

Tables

Table 1. Table 1: Rms error E between the model and the data, computed using ( 45 ), for different values of C 2 subscript 𝐶 2 C_{2} . The associated value of C 1 subscript 𝐶 1 C_{1} (see text for details) is also shown. Chosen fiducial parameters are marked in bold.

$C_2$	$C_1$	E
1.0	34.03	0.670
0.9	5.14	0.632
0.8	3.29	0.601
0.66	2.06	0.503
0.5	1.23	0.339
0.33	0.62	0.308
0.25	0.39	0.453
0.2	0.27	0.590
0.1	0.08	0.988

Equations

$$1<R_0<{\tau }^{-1},$$

$$R_0=\frac{|N_T^2|}{|N_C^2|}$$

$$D_{\mathrm{mix}}=C_t{\kappa }_T/R_0,$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\begin{gathered} {\rho }_m\left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}\right)=-\nabla p+{\rho }_m\nu {\nabla }^2\mathbf{u} \\ +\frac{1}{{\mu }_0}\left(\nabla \times \mathbf{B}\right)\times \mathbf{B}+{\rho }_m\left(-\alpha T+\beta C\right)\mathbf{g}, \end{gathered}$$

$$\frac{\partial T}{\partial t}+\mathbf{u}\cdot \nabla T+u_z\left(\frac{d{T}_0}{dz}-\frac{d{T}_{\mathrm{ad}}}{dz}\right)={\kappa }_T{\nabla }^{2}T,$$

$$\frac{\partial C}{\partial t}+\mathbf{u}\cdot \nabla C+u_z\frac{d{C}_0}{dz}={\kappa }_C{\nabla }^{2}C,$$

$$\frac{\partial \mathbf{B}}{\partial t}=\nabla \times \left(\mathbf{u}\times \mathbf{B}\right)+\eta {\nabla }^2\mathbf{B},$$

$$\nabla \cdot \mathbf{B}=0.$$

$$d={\left(\frac{{\kappa }_T\nu }{\alpha g\left(\frac{d{T}_0}{dz}-\frac{d{T}_{\mathrm{ad}}}{dz}\right)}\right)}^{1/4}={\left(\frac{{\kappa }_T\nu }{N_T^2}\right)}^{1/4},$$

$$[t]=\frac{d^2}{{\kappa }_T},[v]=\frac{{\kappa }_T}{d},$$

$$[T]=d\left(\frac{d{T}_0}{dz}-\frac{d{T}_{\mathrm{ad}}}{dz}\right),[C]=\frac{\alpha }{\beta }[T],$$

$$[B]=B_0,$$

$$\begin{gathered} \frac{\partial \hat{\mathbf{u}}}{\partial t}+\hat{\mathbf{u}}\cdot \nabla \hat{\mathbf{u}}=-\nabla \hat{p}+\mathrm{Pr}{\nabla }^2\hat{\mathbf{u}} \\ +H_B\left(\nabla \times \hat{\mathbf{B}}\right)\times \hat{\mathbf{B}}+\mathrm{Pr}\left(\hat{T}-\hat{C}\right){\hat{\mathbf{e}}}_z, \end{gathered}$$

$$\frac{\partial \hat{T}}{\partial t}+\hat{\mathbf{u}}\cdot \nabla \hat{T}+{\hat{u}}_z={\nabla }^2\hat{T},$$

$$\frac{\partial \hat{C}}{\partial t}+\hat{\mathbf{u}}\cdot \nabla \hat{C}+\frac{{\hat{u}}_z}{R_0}=\tau {\nabla }^2\hat{C},$$

$$\frac{\partial \hat{\mathbf{B}}}{\partial t}=\nabla \times \left(\hat{\mathbf{u}}\times \hat{\mathbf{B}}\right)+D_B{\nabla }^2\hat{\mathbf{B}},$$

$$\nabla \cdot \hat{\mathbf{u}}=0,$$

$$\nabla \cdot \hat{\mathbf{B}}=0,$$

$$\mathrm{Pr}$$

$$R_0$$

$$\mathrm{Pm}=\frac{\nu }{\eta }=\frac{\mathrm{Pr}}{D_B}.$$

$$\hat{q}={\hat{q}}_0\exp [\hat{\lambda }t+i(\hat{\mathbf{k}}\cdot \mathbf{x})],$$

$$F_T=⟨u_{z}T⟩,$$

$$F_C=⟨u_{z}C⟩,$$

$$F_C=-D_C\frac{d{C}_0}{dz}\leftrightarrow D_C=-\frac{⟨u_{z}C⟩}{\frac{d{C}_0}{dz}}.$$

$${\hat{D}}_C=\frac{D_C}{{\kappa }_T}=-R_0{\hat{F}}_C=-R_0⟨{\hat{u}}_z\hat{C}⟩.$$

$${\hat{F}}_C=C_1{\hat{w}}_f{\hat{C}}_f=-C_1\frac{{\hat{w}}_f^2}{R_0({\hat{\lambda }}_f+\tau {\hat{l}}_f^2)},$$

$${\hat{\lambda }}_f=C_2\hat{\sigma },$$

$${\hat{w}}_f=(C_{2}K{)}^{-1}\frac{{\hat{\lambda }}_f}{{\hat{l}}_f},$$

$${\hat{F}}_C=-\left(\frac{C_1}{C_2^2{K}^2}\right)\frac{{\hat{\lambda }}_f^2}{{\hat{l}}_f^2{R}_0({\hat{\lambda }}_f+\tau {\hat{l}}_f^2)},$$

Full text

††software: Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), SciPy (Virtanen et al. 2020)

1 Introduction

1.1 Astrophysical motivation

One of the original governing paradigms of stellar evolution models (Eddington 1926) is that stellar radiation zones are essentially laminar, because their strong stable density stratification quenches vertical turbulent motions. That paradigm has shifted over the last few decades, however, with the realization that there are a number of fluid instabilities that can develop and drive turbulence despite the stratification (see the recent review by Garaud 2021, and references therein). One such instability is the fingering instability (also called thermohaline instability), which can take place in the presence of what is often referred to as an "inverse" composition gradient, i.e., when the mean molecular weight of the star increases with radius. More precisely, a necessary condition for fingering convection (Stern 1960; Baines & Gill 1969; Ulrich 1972) is:

[TABLE]

where the diffusivity ratio $\tau=\kappa_{C}/\kappa_{T}$ is the ratio of the microscopic diffusivity $\kappa_{C}$ of the chemical species which most contributes to the inverse composition gradient, to the thermal diffusivity $\kappa_{T}$, and where the density ratio

[TABLE]

is the ratio of the stable temperature stratification (measured by $N_{T}^{2}>0$, the square of the Brunt-Väisälä frequency associated with the temperature stratification only) to the unstable composition stratification (measured by $N_{C}^{2}<0$, the corresponding quantity for composition). In stellar interiors, the thermal diffusivity is significantly greater than the compositional diffusivity so $\tau^{-1}\gg 1$, and thus even very slight inverse composition gradients may drive fingering convection.

A wide variety of stars are believed to undergo fingering convection, with consequences that are possibly observable, or may otherwise impact their evolution or the interpretation of observations. In low-mass red giant branch (RGB) stars, for example, extra mixing is required to explain abundance observations near the luminosity bump (e.g., Gratton et al. 2000; Shetrone et al. 2019) and is largely thought to be due to fingering convection (first noted by Charbonnel & Zahn 2007, see also Fraser et al. 2022b and references therein). The surface abundances of metal-rich white dwarfs (WDs) are often thought to indicate active accretion (Jura 2003; Farihi et al. 2009; Koester et al. 2014). This creates an unstable composition gradient near the surface that can drive fingering convection. The enhanced turbulent mixing rate must be taken into account in order to infer the accretion rate from observations (Deal et al. 2013; Bauer & Bildsten 2018, 2019; Wachlin et al. 2022). Also in WDs, a destabilizing composition gradient is thought to form when cooling drives crystalization in the core (Stevenson 1980; Mochkovitch 1983). The potential Rayleigh-Taylor instability of this gradient has been put forth as a possible source of dynamo activity to explain observed WD magnetic fields (Isern et al. 2017; Ginzburg et al. 2022); recent work by Fuentes et al. (2023) and Montgomery & Dunlap (2023), however, has demonstrated that this configuration is likely fingering-unstable, with important consequences to the possibility of a crystallization-driven dynamo. Fingering convection is also believed to occur in the envelopes of accretor stars in massive stellar binaries, and to significantly affect their interior composition profiles (Renzo & Götberg 2021). Similarly, fingering convection is thought to be an important process in the formation of carbon-enhanced metal-poor stars by accretion of carbon-rich material from a massive, evolved companion (Stancliffe et al. 2007).

Quantifying the precise impact of fingering convection on stellar evolution requires a model for mixing by the turbulence it creates. Unfortunately, despite the enormous progress made in supercomputing in the past 30 years, it is still not possible to perform three-dimensional direct numerical simulations (DNS) of fingering convection at the parameters appropriate for stellar interiors. Instead, one must rely on simple ad-hoc models, tested (whenever possible) using DNS at parameters that are achievable numerically. We now describe recent progress on that front.

1.2 From traditional to modern models of fingering convection

Early prescriptions for turbulent mixing by fingering convection based on dimensional analysis date back as far as Ulrich (1972) and Kippenhahn et al. (1980). Both studies argued that the turbulent compositional flux can be modeled as a diffusive process, with a diffusion coefficient $D_{\mathrm{mix}}$ that can be written as

[TABLE]

where $C_{t}$ is a free parameter. Plausible values for $C_{t}$ put forth in these original calculations vary by orders of magnitude, from $C_{t}=O(1000)$ (Ulrich 1972) to $C_{t}=O(10)$ (Kippenhahn et al. 1980). Charbonnel & Zahn (2007) found $C_{t}=O(1000)$ to be necessary in order for their 1D stellar evolution models of RGB stars to achieve mixing levels consistent with the observations of Gratton et al. (2000). This was later confirmed by a similar study of Denissenkov (2010) (for a review of efforts to match observations using this model, see Salaris & Cassisi 2017).

Since 2010, progress in high-performance computing has enabled us to estimate compositional fluxes in fingering convection using three-dimensional DNS in numerically tractable parameter regimes and test these models (Garaud 2018). In particular, Traxler et al. (2011) and Brown et al. (2013) found that the fluxes predicted by Ulrich (1972) overestimate the results of their hydrodynamic DNS by very large amounts. As a result, fingering convection alone does not seem to be a strong enough source of mixing to explain observations of RGB star abundances according to the results of Charbonnel & Zahn (2007). Furthermore, the functional dependence of the turbulent diffusivity on $R_{0}$ and $\tau$ predicted by Eq. (3) is not consistent with DNS results, prompting a search for new, better models.

Brown et al. (2013), following similar work in the oceanographic context by Radko & Smith (2012) (which itself built upon the pioneering work by Holyer 1984), developed a "parasitic saturation" model that agrees much more favorably with their simulations. Generally speaking, parasitic models (described in greater detail in Sec. 3.2 below) assume that a given primary instability saturates nonlinearly because of the development of secondary "parasitic instabilities" that feed on their energy (for applications in a variety of contexts beyond fingering convection, see, e.g., Goodman & Xu 1994; Pessah & Goodman 2009; Pessah 2010; Latter et al. 2009; Longaretti & Lesur 2010; Barker et al. 2019, 2020). In the case of fingering convection, Radko & Smith (2012) and Brown et al. (2013) demonstrated that shear instabilities can cause the saturation of the fingering instability when their respective growth rates are similar. Applying the theory yields estimates for the vertical velocity of the fingers at saturation, which can then be used to estimate the vertical heat and compositional transport by fingering convection.

The parasitic model of Brown et al. (2013) correctly predicts the turbulent compositional fluxes measured in hydrodynamic DNS of fingering convection for a wide range of input parameters (see, e.g., Figure 2 of Garaud 2018), suggesting that one could confidently use it at stellar parameters as well. Doing so, however, continues to predict compositional fluxes that are much smaller than those required to explain abundance observations in RGB stars (Charbonnel & Zahn 2007), and potentially challenges the notion that this instability is the answer to this particular observational conundrum. And yet, the fact that extra mixing is needed precisely at the time when an inverse composition gradient forms in the star (see counter-arguments by Tayar & Joyce 2022, however), and that the extra mixing required increases with the magnitude of the composition gradient (Fraser et al. 2022b), as one would expect in fingering convection, makes it very difficult to abandon the idea of fingering convection altogether.

1.3 Magnetized fingering convection

Harrington & Garaud (2019) (HG19 hereafter) recently demonstrated using DNS that a uniform, vertical magnetic field (that is aligned with the direction of gravity) can increase transport in fingering convection by suppressing the parasitic shear instabilities and letting finger velocities grow to larger amplitudes before saturation. They also extended the Brown et al. (2013) model by accounting for the magnetic field in the calculation of the growth rates of the parasitic shear instabilities. They found that this new magnetohydrodynamic (MHD) parasitic model correctly accounts for the increase in compositional flux with magnetic field strength observed in their DNS. The HG19 model predicts that magnetic field strengths as low as $O(100\text{G})$ could increase the efficiency of fingering convection in RGB stars by two orders of magnitude, thus potentially resolving the RGB stars abundance conundrum discussed above.

Despite these promising results, it is important to note that the HG19 work was very preliminary and only explored a very limited range of numerically accessible input parameters: they fixed $R_{0}=1.45$ and $\mathrm{Pr}=\tau=0.1$, where $\mathrm{Pr}$ is the Prandtl number (the ratio of the kinematic viscosity to the thermal diffusivity), and used a magnetic Prandtl number $\mathrm{Pm}=1$, where $\mathrm{Pm}$ is the ratio of kinematic viscosity to magnetic diffusivity. These choices are, in hindsight, somewhat problematic. Indeed, the radiative zones of many stars usually have $\mathrm{Pm}<1$. In addition, fingering regions likely span a wide range of density ratios $R_{0}$ (see, e.g., Fraser et al. 2022b). Thus, the HG19 model merits comparison with DNS over a much wider range of parameters than what has been covered so far.

1.4 Structure of the paper

In this paper we extend the work of HG19 in two ways. First, we perform DNS over a broader range of density ratios $R_{0}$ and at $\mathrm{Pm}<1$ (using the same uniform, vertical background magnetic field – we discuss some caveats pertaining to this choice in Sec. 6). We find that there are regimes where the HG19 model fails because viscosity and resistivity, which were neglected by HG19, begin to dominate the dynamics of the parasitic instability. Second, we extend the HG19 model in an attempt to account for these effects, and compare the results to our DNS. We first add viscosity and resistivity in the computation of the parasitic instability properties, but find that this is unable to predict the fluxes obtained in the DNS. Then, following Radko & Smith (2012), we also add the temperature and composition fields to the parasitic instability model, and find that it now matches the data with high fidelity.

The paper is outlined as follows: we introduce our governing equations, non-dimensionalization, and review the basic properties of the fingering instability in Sec. 2, followed by a summary of previous work in Sec. 3. In Sec. 4, we present the results of our DNS and compare them against the HG19 model. We describe our extension to the HG19 model and compare it to our DNS in Sec. 5. Our conclusions are presented in Sec. 6.

2 Problem setup

2.1 Governing equations

We study fingering convection in MHD using the same model as HG19, with the key features outlined here for clarity. We assume that the typical length scales of turbulent fluctuations are much shorter than the temperature, density and pressure scale heights, and that typical velocity fluctuations are much smaller than the local sound speed. These assumptions allow us to use the Boussinesq approximation for compressible gases (Spiegel & Veronis 1960). Similarly, we assume that the relevant length scales are much smaller than the radius of the star to justify the use of a Cartesian grid $(x,y,z)$, with $z$ taken as the vertical (i.e., radial) direction. The governing equations are then

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

In these equations $\mathbf{u}=(u_{x},u_{y},u_{z})$ is the velocity field, $\mathbf{B}=(B_{x},B_{y},B_{z})$ is the magnetic field, $p$ is pressure perturbation away from hydrostatic equilibrium (which is assumed for the background stratification), and $T$ and $C$ are the temperature and composition perturbations away from assumed linear background profiles with gradients $dT_{0}/dz$ and $dC_{0}/dz$, respectively. All other quantities are model constants (consistent with the use of the Boussinesq approximation). The mean density of the domain is $\rho_{m}$, and $\alpha$ and $\beta$ are the coefficients of thermal expansion and compositional contraction, respectively. The adiabatic temperature gradient is $dT_{\mathrm{ad}}/dz=-g/c_{p}$, where $g=|\mathbf{g}|$ is the local gravity and $c_{p}$ is the specific heat at constant pressure. In stellar radiative zones, $dT_{\mathrm{ad}}/dz<dT_{0}/dz<0$. The kinematic viscosity, thermal diffusivity, compositional diffusivity, and magnetic diffusivity are given by $\nu$, $\kappa_{T}$, $\kappa_{C}$, and $\eta$, respectively, and are also taken to be constant. Finally, $\mu_{0}$ is the permeability of free space. Note that magnetic buoyancy is not included in this model.

The dynamical fields of interest, $\mathbf{u}$, $\mathbf{B}$, $T$, and $C$, are assumed to satisfy periodic boundary conditions along each axis, a convenient property for direct numerical simulations because it enables the use of Fourier series expansions in pseudospectral codes, and avoids unphysical boundary layers that inevitably occur near impermeable boundaries.

We then non-dimensionalize the equations as in HG19. Length scales are measured in units of the characteristic width of fingers (Stern 1960),

[TABLE]

where $N_{T}$ is the local Brunt–Väisälä frequency defined from the temperature stratification. The other units for time, velocity, temperature, composition, and magnetic field are, respectively,

[TABLE]

[TABLE]

and

[TABLE]

where $B_{0}$ is the amplitude of an assumed mean vertical field (see Sec. 1 and below).

The non-dimensional governing equations are (cf. HG19)

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where hats over dynamical fields indicate non-dimensional quantities, and $\hat{\mathbf{e}}_{z}$ is the unit vector in the $z$ direction. From here on, unless otherwise stated, space and time coordinates are implicitly assumed to be non-dimensionalized in this way. The system is now specified by five dimensionless parameters:

[TABLE]

where $\mathrm{Pr}$ is the Prandtl number, $\tau$ and $D_{B}$ are the compositional and resistive diffusivity ratios, respectively, $R_{0}$ is the density ratio, and $H_{B}$ is the coefficient of the Lorentz force, which can be understood as the squared inverse of the Alfvén Mach number based on the unit magnetic field and the unit flow speed. In stars, the parameters $\mathrm{Pr}$, $\tau$, and $D_{B}$ are all extremely small (see, e.g., Garaud et al. 2015). The density ratio characterizes the stabilizing influence of the stable temperature gradient relative to the destabilizing composition gradient, and, in radiation zones, can take any value from 1 to infinity depending on local stellar conditions. The threshold $R_{0}=1$ is equivalent to the Ledoux criterion for convection, and $R_{0}>1$ applies to regions that are stable to convection but potentially unstable to fingering convection, as described below in Sec. 2.2. Throughout this paper, we will often refer to the magnetic Prandtl number, defined as

[TABLE]

This ratio is typically $O(10^{-2})-O(1)$ in the radiative zones of most stars.

2.2 Linear stability analysis

A thorough discussion of the linear stability properties of non-magnetic fingering convection at low $\mathrm{Pr}$ (the regime relevant to stellar interiors) can be found in Brown et al. (2013), with the magnetized case first explored in Harrington (2018) and HG19. The most important result is that for a given diffusivity ratio $\tau$, fingering convection takes place when $1<R_{0}<1/\tau$, with the system becoming less unstable as $R_{0}$ increases (Stern 1960; Baines & Gill 1969).

More specifically, eigenmode solutions to the linearized system for small-amplitude perturbations take the form

[TABLE]

where $q$ represents each of the dynamical fields, $\hat{\mathbf{k}}=(\hat{k}_{x},\hat{k}_{y},\hat{k}_{z})$ is the mode’s wavenumber, and $\hat{\lambda}$ is its growth rate. The fastest-growing modes are so-called “elevator modes", with a velocity that is in the vertical direction only, and where all perturbations are invariant in $z$ (so $\hat{k}_{z}=0$). Harrington (2018) and HG19 showed that the same is true for magnetized fingering convection with a uniform background magnetic field regardless of its orientation. When the magnetic field is vertical, which is the case considered here, the field lines are parallel to the fluid flow within each elevator, and therefore do not interact with the latter. As a result, the properties of the fastest-growing modes (growth rate and wavenumber) are identical in MHD and in the hydrodynamic limit. The magnetic field, however, stabilizes fingering modes with $\hat{k}_{z}\neq 0$ and also plays a crucial role in the nonlinear saturation of the elevators/fingers by stabilizing them against secondary shear instabilities (see HG19). This in turn controls the typical velocity of the fluid within the fingers after saturation, and sets the transport rate of heat and chemical species by the fingering convection.

2.3 Turbulent diffusion by fingering convection

The rate at which a turbulent fluid transports temperature and composition is quantified by the thermal and compositional fluxes, defined (dimensionally) as

[TABLE]

and

[TABLE]

respectively, where $\langle.\rangle$ denotes a volume average as well as a time average taken once the turbulence has reached a statistically stationary state (see Sec. 4.1). Note that both fluxes are negative in fingering convection (Radko 2013). In stellar evolution, it is customary to assume that the turbulent compositional flux can be modeled as a turbulent diffusion process, with a turbulent diffusivity $D_{C}$ defined (dimensionally) such that

[TABLE]

In the units of Section 2.1, this becomes

[TABLE]

In other words, modeling the turbulent diffusion coefficient $D_{C}$ for magnetized fingering convection in stars requires estimating the nondimensional flux $\hat{F}_{C}$ as a function of the local properties of the star, which control the non-dimensional input parameters of the model, namely $\mathrm{Pr}$, $\tau$, $R_{0}$, $H_{B}$, and $D_{B}$. In what follows, we first summarize recent preliminary results in that direction obtained by HG19, who used a combination of numerical simulations and semi-analytical modelling to predict $\hat{F}_{C}$ for a narrow range of parameter space. We then expand their work to a much wider range of parameter space.

3 Summary of previous work

3.1 Previous numerical experiments

The first DNS of fingering convection in the presence of a magnetic field, solving equations (14)-(19), were presented by Harrington (2018) and HG19. The majority of their simulations assumed a uniform vertical background magnetic field, as we do here, with an amplitude selected so that $H_{B}$ spans the interval $[10^{-2},10^{2}]$. In that preliminary work, all the other parameters were held constant, with the diffusivity ratios fixed to be $\mathrm{Pr}=\tau=0.1$ and $\mathrm{Pm}=1$, and the density ratio set to $R_{0}=1.45$. As discussed in Sec. 1, HG19 found that increasing the amplitude of the background magnetic field can significantly enhance the magnitude of the turbulent vertical fluxes of both heat and composition, a result they explained by noting that the field stabilizes the secondary shear instability that normally saturates the primary fingering instability, letting the fingers grow to larger amplitude before saturation. In order to quantify this idea, they proposed a new parasitic saturation model for fingering convection in MHD that is built upon the one successfully presented by Brown et al. (2013) in the hydrodynamic limit.

3.2 Parasitic saturation models

As described in Sec. 1, parasitic saturation models assume that the primary fingering instability (which is dominated by the fastest-growing elevator modes) saturates when secondary shear instabilities developing between the upflowing and downflowing elevators take over. If we let $\hat{w}_{f}$ and $\hat{C}_{f}$ be the amplitude of the vertical velocity and composition perturbation in the fastest-growing elevator mode at saturation, then the corresponding turbulent compositional flux can reasonably be assumed to be proportional to the product $\hat{w}_{f}\hat{C}_{f}$ (Radko & Smith 2012; Brown et al. 2013), i.e.,

[TABLE]

where $C_{1}$ is a model parameter that is assumed to be a universal constant, $\hat{\lambda}_{f}$ and $\hat{l}_{f}$ are the nondimensional growth rate and horizontal wavenumber of the fastest-growing elevator mode, and where in the second equality, $\hat{C}_{f}$ has been expressed in terms of $\hat{w}_{f}$ by using the linearized equations described in Sec. 2.2 (see Brown et al. 2013, Sec. 4.2 for details). Note that $\hat{F}_{C}<0$, as required; in addition, note that this expression is only valid in the fingering range, i.e. for $1<R_{0}<\tau^{-1}$.

The parasitic saturation condition is expressed mathematically as a rough balance between the growth rate of the fingering elevator mode $\hat{\lambda}_{f}$ and the growth rate of the fastest-growing parasitic mode $\hat{\sigma}$, namely

[TABLE]

where $C_{2}$ is another model parameter assumed to be constant. Because shear instabilities almost always grow faster when the shear is stronger, $\hat{\sigma}$ is (usually) a monotonically increasing function of $\hat{w}_{f}$. Eq. (28) can then be inverted to obtain $\hat{w}_{f}$ as a function of $\hat{\lambda}_{f}$, which then yields $\hat{F}_{C}$ as a function of known quantities only via Eq. (27).

In the hydrodynamic case, Brown et al. (2013) ignored the effect of viscosity on the development of the shear instability, in which case the computation of $\hat{\sigma}$ can be done semi-analytically to yield $\hat{\sigma}=K\hat{w}_{f}\hat{l}_{f}$, where $K\simeq 0.26$ is a known constant that is obtained numerically by solving the linear stability problem. Combining this result with Eq. (28) yields

[TABLE]

and provides an estimate of the vertical velocity of the elevator modes at saturation. Sengupta & Garaud (2018) compared this expression with the measured $\hat{u}_{z,\mathrm{rms}}$ in DNS to estimate the prefactor as $(C_{2}K)^{-1}\simeq 2\pi$, and hence $C_{2}\simeq 0.6$. Finally, substituting this into Eq. (27) yields

[TABLE]

which was fitted by Brown et al. (2013) against the available hydrodynamic data for $\hat{F}_{C}$ to estimate $C_{1}(C_{2}K)^{-2}\simeq 49$, and thus $C_{1}\simeq 1.24$. This expression, as demonstrated by Brown et al. (2013) and Garaud (2018), provides a very good estimate of the compositional flux $\hat{F}_{C}$ over most of parameter space.

In the magnetic case, the shear instability growth rate $\hat{\sigma}$ also depends on the magnetic field strength, via $H_{B}$, and (possibly) on the magnetic diffusivity, via $D_{B}$. Following Brown et al. (2013), HG19 chose to consider the so-called "ideal MHD" limit where now both viscosity and resistivity are neglected in the calculation of $\hat{\sigma}$. They showed that, in this limit, the growth rate of the dominant parasitic mode, normalized to the elevator mode amplitude and wavenumber, namely

[TABLE]

is purely a function of

[TABLE]

which is the squared ratio of the Alfvén velocity to the saturated elevator mode velocity. Moreover, they found that this function is approximated well by the expression

[TABLE]

(and corresponds to decaying modes otherwise), which is equivalent to

[TABLE]

This shows that $\hat{\sigma}$ is a monotonically increasing function of $\hat{w}_{f}$ and a monotonically decreasing function of $H_{B}$, and confirms that a vertical magnetic field can stabilize the parasitic shear instability. As a result, the finger velocity $\hat{w}_{f}$ is allowed to grow to larger amplitude in the presence of a magnetic field before the shear instability growth rate $\hat{\sigma}$ reaches the saturation condition Eq. (28). This in turn leads to an increase in the magnitude of $\hat{F}_{C}$ according to Eq. (27).

HG19 compared their model predictions for the compositional flux $\hat{F}_{C}$ (obtained by solving Eqs. (28) and (34) numerically to obtain $\hat{w_{f}}$ as a function of $\hat{\lambda}_{f}$, and then substituting the result into Eq. (27)) to the values of $\hat{F}_{C}$ extracted from their DNS. They found that the two agree very well, at least for the limited set of parameters considered ($\mathrm{Pr}=\tau=D_{B}=0.1$, $R_{0}=1.45$). In what follows, we now extend their work to a much wider range of parameter space, and present new DNS of magnetized fingering convection that explore the effects of increasing the density ratio $R_{0}$ (to study less-unstable systems) and reducing the magnetic Prandtl number $\mathrm{Pm}$ (which reduces the dynamical connection between the field and the flow). Both extensions are important and relevant for stars, where $R_{0}$ necessarily becomes large near the boundary between a fingering-unstable region and a stable radiation zone (Fraser et al. 2022b), and $\mathrm{Pm}$ is typically smaller than 1 in stellar radiative zones.

4 Simulation results

4.1 Methodology

Our simulations are performed using the pseudospectral code PADDI (Stellmach et al. 2011), extended to solve equations (14)-(19) as detailed in Harrington (2018) and used by HG19. Simulations are initialized with a uniform vertical magnetic field of unit strength, $\hat{\mathbf{B}}=(0,0,1)$, and all other dynamical fields set to zero. Low-amplitude white noise perturbations are added to $\hat{T}$ and $\hat{C}$ to seed the instability.

For all new simulations presented here, the horizontal dimensions of the simulation domain are taken as $L_{x}=L_{y}=100$, which is adequate to accommodate many finger widths. However, the domain height $L_{z}$ varies across simulations to ensure that it is always substantially larger than the natural finger height (as determined by visual inspection of simulation snapshots). In the majority of cases with extreme finger heights (larger $H_{B}$ and $R_{0}$, as described in Sec. 4.2 below), we performed additional simulations doubling $L_{z}$ to ensure that the saturated fluxes only change minimally in the process. We use domain heights as low as $L_{z}=100$ and as large as $L_{z}=800$ in the most extreme cases. Further details on the effects of $L_{z}$ in such simulations are given in Appendix A of Traxler et al. (2011). Cubic simulations with $L_{z}=100$ use a numerical resolution of $128^{3}$ Fourier modes (with standard dealiasing according to the 3/2 rule, so that nonlinearities are calculated in physical space on a $192^{3}$ grid)111All $R_{0}=1.45$ MHD simulations presented in this paper are taken from HG19, who used $192^{3}$ Fourier modes, with nonlinearities calculated in physical space on a $288^{3}$ grid., and simulations with taller domains use correspondingly larger numbers of Fourier modes in $z$ to ensure a fixed resolution. We check that this resolution is appropriate for all simulations by inspecting fluctuating fields and verifying that grid-scale oscillations (characteristic of the Gibbs phenomenon) are minimal.

As expected from the linear theory of Sec. 2.2, we find in all cases that the fingering instability undergoes an exponential growth phase that eventually saturates into a quasi-stationary state of homogeneous fingering convection (see, e.g., HG19, Fig. 2). We extract the compositional flux by taking a time-average of $\hat{F}_{C}$ over this quasi-stationary state, and a measure of the error is taken to be the rms variability of the flux about this average; note, however, that in all figures presented below, the error bars are smaller than the symbol sizes themselves, and are thus entirely obscured.

4.2 Trends with increased $R_{0}$

We start by discussing the effect of increasing the density ratio $R_{0}$, while holding the other physical parameters fixed as $\mathrm{Pr}=\tau=0.1$ and $\mathrm{Pm}=1$. Figure 1 shows how $|\hat{F}_{C}|$ varies with $R_{0}$ for the hydrodynamic case (black points) and in MHD for two different values of $H_{B}$ (colored symbols), with the parasitic saturation models described in Sec. 3.2 shown as dashed lines for comparison.

For fixed $H_{B}$, we see that $|\hat{F}_{C}|$ decreases as $R_{0}\to 1/\tau$. This is consistent with expectations from linear theory (Baines & Gill 1969, HG19), which shows that the fingering instability is gradually stabilized as $R_{0}$ increases. We also see that $|\hat{F}_{C}|$ increases with $H_{B}$ at fixed $R_{0}$, consistent with the findings of HG19. Our new simulations reveal that this effect is much more pronounced at larger $R_{0}$, where $|\hat{F}_{C}|$ can increase by more than one order of magnitude compared with the hydrodynamic case, even for moderate $H_{B}$. This can be explained at least qualitatively by the theory proposed by HG19 (see the dashed lines in Fig. 1), through the suppression of the parasitic mode by the magnetic field, although the latter clearly overestimates the effect.

The enhanced effect of magnetic fields on the turbulence at large $R_{0}$ relative to small $R_{0}$ can also be seen qualitatively in Fig. 2, which shows snapshots of $\hat{u}_{z}$ in the $x$-$z$ plane for the MHD simulations presented in Fig. 1. At $R_{0}=1.45$, we see that the velocity fields in the $H_{B}=0.01$ and $H_{B}=0.1$ simulations have similar structures. At larger $R_{0}$, however, the fingers are significantly longer and have higher flow speeds at larger $H_{B}$ than at smaller $H_{B}$, consistent with the findings of Fig. 1 (see above).

4.3 Trends with decreased $\mathrm{Pm}$

In stellar interiors, the magnetic diffusivity is typically somewhat larger than the kinematic viscosity, so the magnetic Prandtl number is roughly $\mathrm{Pm}\sim 0.01-1$. Meanwhile HG19 only presented simulations for $\mathrm{Pm}=1$. Yet, it is well-known that MHD turbulence in $\mathrm{Pm}<1$ fluids behaves very differently than at $\mathrm{Pm}=1$ (see, e.g., Christensen & Aubert 2006; Warnecke et al. 2023, for prominent examples in the stellar dynamo context). Thus, we now investigate the impact of using a smaller value of $\mathrm{Pm}$.

The impact of reducing $\mathrm{Pm}$ can be seen qualitatively in Fig. 3, which shows snapshots of $\hat{u}_{z}$ as in Fig. 2, with $\mathrm{Pr}=\tau=0.1$, $R_{0}=7$, and, in panels (a) - (c), $H_{B}=0.1$ with $\mathrm{Pm}$ decreasing from $1$ to $0.01$. Panel (d) shows a hydrodynamic simulation for comparison. We see that reducing $\mathrm{Pm}$ leads to shorter, lower-amplitude fingers, similar to the effect of reducing $H_{B}$ shown in Fig. 2. At $\mathrm{Pm}=0.01$, the dynamics are nearly indistinguishable from the hydrodynamic case.

To look at this effect more quantitatively, we show in the left panel of Fig. 4 the turbulent compositional flux against the density ratio $R_{0}$ (as in Fig. 1) for simulations with $\mathrm{Pr}=\tau=0.1$, for both hydrodynamic simulations (black dots) and MHD simulations (colored crosses) with fixed $H_{B}=0.1$ and $\mathrm{Pm}$ varying from $1$ to $0.01$. Consistent with Fig. 3, we see that $|\hat{F}_{C}|$ decreases as $\mathrm{Pm}$ decreases, with the fluxes in $\mathrm{Pm}=0.01$ simulations nearly identical to the hydrodynamic case, except for the largest value of $R_{0}$. The predicted fluxes from both hydrodynamic (black, Brown et al. (2013)) and MHD (grey, HG19) parasitic saturation models are also shown as dashed curves for comparison with the data. Because HG19 ignored viscosity and resistivity in their calculation of $\hat{\sigma}$ (see Sec. 3.2), their model does not depend on $\mathrm{Pm}$, which is why there is only one corresponding dashed curve. It is therefore clear that their model, by construction, cannot capture the strong dependence of $\hat{F}_{C}$ on $\mathrm{Pm}$ observed in our DNS data, especially at large $R_{0}$.

One may wonder why the HG19 model can successfully explain the properties of fingering convection at their selected parameter values ($\mathrm{Pr}=\tau=0.1$, $\mathrm{Pm}=1$ and $R_{0}=1.45$), while dramatically failing in our new higher $R_{0}$, lower $\mathrm{Pm}$ simulations. One plausible explanation for this breakdown in the model is provided in the right panel of Fig. 4, where the data points represent one possible measure of the magnetic Reynolds number $\mathrm{Rm}$ associated with the vertical velocity of the fingers, for each of the simulations presented on the left panel. We define this magnetic Reynolds number in terms of the rms vertical velocity $\hat{u}_{z,\mathrm{rms}}$ and the wavenumber of the fastest-growing elevator modes $\hat{l}_{f}$, as

[TABLE]

We find that while $\mathrm{Rm}$ is greater than one in the $\mathrm{Pm}=1$ simulations, the $\mathrm{Pm}=0.1$ and $\mathrm{Pm}=0.01$ simulations all have $\mathrm{Rm}\lesssim 1$. Thus, it is no longer appropriate to ignore the effects of resistivity (and in some cases, even of viscosity) in the calculation of the shear instability growth rate $\hat{\sigma}$, which may explain why the HG19 model fails at these parameters. In what follows, we now extend the HG19 model to account for additional effects including dissipation and buoyancy in the computation of the parasitic modes.

5 Extending the parasitic saturation model

5.1 The extended model

Here we extend the HG19 model to take into account diffusive effects (including viscosity and resistivity) as well as temperature and composition fluctuations on the development of the parasitic instability. As we will see, including diffusion alone is insufficient to bring the model into agreement with simulations, and one must additionally include temperature ($T$) and composition ($C$) fluctuations as was done by Radko & Smith (2012).

As in Brown et al. (2013) and HG19, we use the frozen-in approximation to assume that the fastest-growing elevator mode has a two-dimensional, sinusoidal, steady velocity field of the form $\hat{\mathbf{u}}_{0}=\hat{w}_{f}\cos(\hat{l}_{f}x)\hat{\mathbf{e}}_{z}$ (the choice of using $x$ as the oscillating direction is made here without loss of generality). It is embedded in a uniform, streamwise magnetic field $\hat{\mathbf{B}}_{0}=\hat{\mathbf{e}}_{z}$. For sufficiently large values of $\hat{w}_{f}$, the elevator mode becomes unstable to a parasitic instability whose growth rate $\hat{\sigma}$ can be computed by performing a standard linear stability calculation. Ignoring buoyancy fluctuations, this calculation would be the same as the stability calculation that was performed in Fraser et al. (2022a), albeit with a different non-dimensionalization.

Fraser et al. (2022a) found that 2D perturbations are always the most unstable in their system. We thus assume that this is also true here, and consequently, that the fastest-growing parasitic modes can be described using a streamfunction formulation where the perturbed flow and field are given by $\mathbf{\hat{u}}=\nabla\times\left(\hat{\psi}\hat{\mathbf{e}}_{y}\right)$ and $\mathbf{\hat{b}}=\nabla\times\left(\hat{A}\hat{\mathbf{e}}_{y}\right)$, respectively. We split the streamfunction $\hat{\psi}$ into elevator and perturbation as follows: $\hat{\psi}=\hat{\psi}_{E}+\hat{\psi}_{P}$ (where the subscript $E$ represents the elevator mode and $P$ the perturbation), and similarly let $\hat{T}=\hat{T}_{E}+\hat{T}_{P}$ and $\hat{C}=\hat{C}_{E}+\hat{C}_{P}$. Since the elevator mode does not affect the background field, we have $\hat{A}=\hat{A}_{P}$. As the elevator mode is invariant in $z$ and sinusoidal in $x$, we express it as:

[TABLE]

The elevator mode amplitude ($\hat{E}_{\psi}$) is related to $\hat{w}_{f}$ by the equation $2\hat{E}_{\psi}\hat{l}_{f}=\hat{w}_{f}$, where $\hat{l}_{f}$ is the horizontal wavenumber of the fastest-growing fingering mode for given values of $\mathrm{Pr}$, $\tau$ and $R_{0}$ ($\hat{l}_{f}$ can easily be computed from linear theory, see section 2.2). The $\hat{T}$ and $\hat{C}$ components of the elevator mode are similarly given by (cf. Eq [3.10] of Radko & Smith 2012)

[TABLE]

and

[TABLE]

where $\hat{E}_{T}$ and $\hat{E}_{C}$ are related to $\hat{E}_{\psi}$ (and therefore $\hat{w}_{f}$) via

[TABLE]

On their frozen-in spatially sinusoidal background state, Fraser et al. (2022a) found that the fastest growing linear perturbations have the same spatial period as the elevator mode. Here we assume the same, thus perturbations take the form

[TABLE]

where $\hat{k}_{z}$ is the vertical wavenumber of the perturbation, and $\hat{\sigma}$ is its growth rate. In practice, we truncate the series to a finite sum from $-N$ to $N$ so that we can numerically compute $\hat{\sigma}$ for a given $\hat{k}_{z}$. Using (40) in the governing equations (14)-(19), linearizing the system (assuming perturbations have small amplitudes compared with the background state), and projecting the results onto each Fourier mode yields the following set of linear equations:

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where we have defined $\hat{k}_{m}^{2}=m^{2}\hat{l}_{f}^{2}+\hat{k}_{z}^{2}$. For each given values of $\hat{k}_{z}$ and $\hat{w}_{f}$ and input parameters ($\mathrm{Pr}$, $\tau$, $R_{0}$, $H_{B}$, $D_{B}$), Eqs. (41) - (44) define an eigenvalue problem with possibly multiple eigenvalues $\hat{\sigma}$. For each $\hat{k}_{z}$, we only consider the eigenvalue with the largest real part and call this $\hat{\sigma}(\hat{k}_{z},\hat{w}_{f};\mathrm{Pr},\tau,R_{0},H_{B},D_{B})$. We then maximize the result over all $\hat{k}_{z}$ to obtain the fastest-growing parasitic mode growth rate as a function of the elevator mode velocity and input parameters, $\hat{\sigma}_{\rm max}(\hat{w}_{f};\mathrm{Pr},\tau,R_{0},H_{B},D_{B})$. Finally, knowing the elevator mode’s growth rate $\hat{\lambda}_{f}$ (which is again a known function of $\mathrm{Pr}$, $\tau$ and $R_{0}$) we find the smallest value of $\hat{w}_{f}$ for which Eq. (28) holds (corresponding to the largest amplitude the elevator mode can reach before being disrupted by the parasitic mode), and use the result to calculate $\hat{F}_{C}$ using Eq. (27). For validation, we have checked (in the hydrodynamic case) that this procedure exactly recovers Figure 7 of Radko & Smith (2012) with the parameters $\mathrm{Pr}=7$, $\tau=0.01$ (setting $C_{2}=1/4.3$ and $C_{1}=0.5$ to match their model parameters).

As in the HG19 model, the entire procedure depends on two parameters $C_{1}$ and $C_{2}$ which appear in equations (27) and (28) respectively. For now, we choose $C_{1}=0.62$ and $C_{2}=0.33$ as our fiducial parameters (we revisit this matter in Sec. 5.2).

The model predictions are compared with the data for a range of $H_{B}$, $\mathrm{Pm}$, and $R_{0}$ values and with $\mathrm{Pr}=\tau=0.1$ in Fig. 5, where the top-left panel shows the dissipative model without $T$ and $C$ fluctuations included (cf. Fraser et al. 2022a), and the top-right panel includes the full system described by Eqs. (41)-(44). In both panels, dissipative models are given by solid lines and DNS data by crosses. The number of Fourier modes used for this figure is 9 for the results without $\hat{T}$ and $\hat{C}$ (corresponding to $N=4$) and 21 when $\hat{T}$ and $\hat{C}$ are included (corresponding to $N=10$). These numbers were chosen so that adding additional modes did not change the results significantly. The top-left panel shows that adding dissipation does not, on its own, help resolve the discrepancy between model and DNS. Comparing to Fig. 4, we see that adding viscosity significantly worsens the fit between the model and the data at large $R_{0}$ even in the hydrodynamic case. This is because the predicted value of $\hat{w}_{f}$ would normally decreases as $R_{0}$ approaches $1/\tau$ in the non-dissipative model (and in the DNS data), but must exceed a certain critical value for the parasitic modes to overcome viscous damping in the dissipative model, and is hence significantly overestimated. For this reason, the model prediction for $|\hat{F}_{C}|$ largely overestimates the DNS results in that limit.

In the magnetic case the problem is similar (with $|\hat{F}_{C}|$ being vastly over-estimated at large $R_{0}$), although the dissipative model is somewhat better at capturing the data trends than the non-dissipative HG19 model. In particular, we see that the dissipative model correctly predicts that $|\hat{F}_{C}|$ should decrease as $\mathrm{Pm}$ decreases, but the predicted sensitivity of $\hat{F}_{C}$ to changes in $\mathrm{Pm}$ is much weaker than what is seen in DNS, particularly at large $R_{0}$.

The results of the full model with $T$ and $C$ fluctuations included is given in the top-right panel of Fig. 5. Focusing first on the hydrodynamic case (black), we see that including temperature and composition in the parasitic mode calculation greatly improves the accuracy of the prediction for $\hat{F}_{C}$ over that of the dissipative model without them, to the extent that it now recovers the quality of the original non-dissipative model of Brown et al. (2013)—see Fig. 4. But more importantly, we also see that the new model is far better at capturing the dependence of $\hat{F}_{C}$ on $H_{B}$, $R_{0}$, and $\mathrm{Pm}$ in the magnetic case. In particular, we now correctly predict that the turbulent compositional flux at $H_{B}=0.1$ tends to its corresponding hydrodynamic value when $\mathrm{Pm}\rightarrow 0$, which is expected as the magnetic field has an increasingly small effect on the flow in this limit.

In hindsight, we can see why including the $\hat{T}$ and $\hat{C}$ fields might be conceptually important as follows. If we ignore $\hat{T}$ and $\hat{C}$ in the dynamics of the parasitic mode, shear is only possible source of instability, and is easily stabilized by viscosity or magnetic fields when the elevator mode amplitude $\hat{w}_{f}$ is small. For this reason, $\hat{w}_{f}$ must always exceed a certain threshold for the parasitic mode to grow, which leads us to overestimate $|\hat{F}_{C}|$ in some regions of parameter space. However, when $\hat{T}$ and $\hat{C}$ are included, low-$\hat{k}_{z}$ parasitic modes can draw energy from the unstable background compositional stratification just like the basic fingering instability does, and can therefore be unstable even at very low (or even zero) $\hat{w}_{f}$. This seemed to have already been known by Radko & Smith (2012) who had included the $\hat{T}$ and $\hat{C}$ fields in their original study of the parasitic modes and their role in saturating the hydrodynamic fingering instability in the geophysical context.

5.2 Model parameters

As mentioned before, the model presented has two free parameters, which are assumed to be universal constants. The first, $C_{1}$, appears in Eq. (27) and controls the level of chemical mixing driven by elevator modes of a given amplitude $\hat{w}_{f}$. Adjusting $C_{1}$ simply corresponds to uniformly shifting the curves shown in Fig. 5 vertically up or down. The second parameter, $C_{2}$, represents the assumed ratio between the growth rate of the elevator mode and its parasitic perturbation at saturation (see equation Eq. 28). Its effect on the prediction for $\hat{F}_{C}$ is more subtle, as shown in Fig. 5 (bottom row).

Two new values of $C_{2}$ are tested: $C_{2}=0.8$ (lower left) and $C_{2}=0.2$ (lower right). In both cases, $C_{1}$ is adjusted a posteriori to obtain the best fit with the data, see Table 1 for detail. This figure shows that, generally speaking, a larger value of $C_{2}$ causes the model to be more sensitive to increases in $H_{B}$ and $\mathrm{Pm}$. This can be understood fairly intuitively as follows. For small values of $C_{2}$, the parasitic mode growth rate $\hat{\sigma}$ needs to reach a correspondingly larger value before saturation can occur, which therefore requires a larger amplitude $\hat{w}_{f}$. In turn, this implies that the effective magnetic Reynolds number of the parasitic mode is larger, reducing the effects of the magnetic diffusivity. Similarly, the elevator mode kinetic energy is higher, so the ambient magnetic energy is less significant in comparison, and the field has less impact on the parasitic perturbations. For larger values of $C_{2}$ the opposite is true. Comparing the three figures, we find that $C_{2}=0.33$ is roughly optimal because it has the right sensitivity to both $H_{B}$ and $\mathrm{Pm}$.

In practice, we have chosen the values $C_{1}=0.62$ and $C_{2}=0.33$ as our fiducial parameters using the following procedure. We define the rms error between the model and the data as follows:

[TABLE]

where the sum ranges over the available DNS data (which has $\mathrm{Pr}=\tau=0.1$) and $N_{p}$ is the number of data points. Note that we use the logarithm of the flux to ensure that all data points carry a similar weight. For a given value of $C_{2}$, we first compute the model flux $\hat{F}_{C,\text{model}}^{*}$ assuming that $C^{*}_{1}=1$. Then the multiplicative factor $C_{1}$ that best fits the data for that given $C_{2}$ satisfies

[TABLE]

Finally, we compute E using (45). The results are shown in Table 1 for a range of $C_{2}$. We can see among all of the values of $C_{2}$ tested, the choice $C_{2}=0.33,C_{1}=0.62$ minimizes the rms error E.

Finally, note that although $C_{2}$ and $C_{1}$ were tuned using DNS data for the compositional flux $\hat{F}_{C}$ only, the model is equally accurate in predicting the thermal flux $\hat{F}_{T}$ with the same model parameters, where $\hat{F}_{T}$ is calculated from $\hat{w}_{f}$ according to

[TABLE]

just as Eq. (27) is used to calculate $\hat{F}_{C}$ (Radko & Smith 2012; Brown et al. 2013). This is shown in Figure 6, which compares the model and the data over a wide range of input parameters for both $\hat{F}_{C}$ and $\hat{F}_{T}$, using the tuned parameters ($C_{2}=0.33,C_{1}=0.62$). We see that the model is accurate at predicting both fluxes to within a factor of two (i.e. within the grey lines) for almost all cases.

6 Summary and Conclusions

6.1 Summary

We have performed a suite of DNS of magnetized fingering convection at $\mathrm{Pr}=\tau=0.1$ with an initially uniform, vertical magnetic field, and compared the measured turbulent compositional fluxes to predictions from the parasitic saturation model put forth by HG19. Our work extends HG19’s simulation study (which used the same $\mathrm{Pr}$ and $\tau$, and fixed the density ratio at $R_{0}=1.45$, and the magnetic Prandtl number at $\mathrm{Pm}=1$) to higher $R_{0}$ and lower $\mathrm{Pm}$, which are both relevant to stellar interior conditions.

As in HG19, we found that the large-scale vertical magnetic field always enhances vertical transport by fingering convection. This is especially true at large $R_{0}$, where it elongates fingers and increases transport even more than at small $R_{0}$. However, while predictions of the simple parasitic saturation model proposed by HG19 are qualitatively consistent with this trend, we have found that they significantly over-estimate the magnitude of the turbulent fluxes measured in our DNS at large $R_{0}$. Lowering $\mathrm{Pm}$ has the opposite effect: in our simulations with $\mathrm{Pm}<1$, the effect of magnetic fields on the fingers becomes less pronounced and the turbulent transport is correspondingly lower than when $\mathrm{Pm}=1$. This is because the resistivity causes the magnetic field to decouple from the fluid, at which point the turbulent flow behaves essentially as in the corresponding hydrodynamic simulations. Because the HG19 model neglects viscosity and resistivity in the calculation of the parasitic shear instability growth rate $\hat{\sigma}$, its predictions do not depend on $\mathrm{Pm}$, so it is unable, by construction, to capture this trend.

By computing the magnetic Reynolds number for the fingers, $\mathrm{Rm}$, we found that the HG19 model fails to quantitatively predict the correct turbulent fluxes whenever $\mathrm{Rm}$ is order unity or smaller. Thus, magnetic diffusion is important to the dynamics of the fingers and their parasitic shear instabilities and cannot be neglected (at least in the parameter regimes covered by the DNS).

In an effort to remedy the problem, we extended the HG19 model to include (a) viscosity and resistivity, and (b) temperature and composition fluctuations in the calculation of the growth rates of the parasitic modes. While these effects were unnecessary to explain the data in Brown et al. (2013) and HG19 and thus deemed negligible, we found that both of these additions are necessary to reproduce the data when considering the wider range of parameter space explored in this paper. Accounting for these effects in the parasitic saturation model now provides an excellent match with the turbulent thermal and compositional flux data for all available simulations (see Fig. 6). We therefore conclude that the model presented in Section 5.1 can reliably be used to estimate mixing by magnetized fingering convection in the presence of a uniform, vertical magnetic field for a very wide range of possible input parameters ($R_{0},\mathrm{Pr},\tau,D_{B}$ and $H_{B}$).

6.2 Model caveats

There are several caveats one must bear in mind, however. First, a uniform, vertical magnetic field is a very special geometric configuration. While the uniform field assumption is justified by the extremely small-scale nature of the fingers compared to the pressure scale height (and thus presumably the gradient scale lengths of whatever magnetic fields exist in these stellar radiation zones; see Harrington & Garaud 2019, Table 1), fingers are in general more likely to encounter oblique fields than vertical ones. HG19 demonstrated that the fastest-growing mode of the fingering instability is unchanged if the field is inclined, but that the nonlinear saturation process is affected. More specifically, they saw in several examples that the turbulent mixing in inclined field cases can greatly exceed that of the corresponding vertical field configuration (see Harrington 2018, for examples with a horizontal field). If this result holds in general, it implies that the parasitic saturation model proposed in this paper only provides a lower bound to the magnetized fingering flux amplitudes. Future work will therefore be needed to quantify the effects of field inclination on fingering convection.

Second, great care must be taken in the application of the parasitic saturation model to ensure that, for any given set of input parameters, the predictions respect the constraints of the Boussinesq approximation. More specifically, one must check that the dimensional elevator mode velocity at saturation, namely $w_{f}=\hat{w}_{f}\kappa_{T}/d$, always remains much smaller than the local sound speed, and that the characteristic vertical scale of the fingers remains much smaller than the local pressure, density or temperature scale height. Both of these self-consistency checks can easily be done a posteriori. Note that a rough estimate of finger height is $d\hat{k}_{z,{\rm max}}^{-1}$, where $\hat{k}_{z,{\rm max}}^{-1}$ is the nondimensional vertical wavenumber of the fastest-growing parasitic mode at saturation.

Finally, we note that the computation of the fingering fluxes in our complete magnetized parasitic saturation model (including buoyancy and diffusion) is computationally much more onerous than in Brown et al. (2013) or HG19. As such, it is not entirely practical for stellar evolution calculations that need to be performed for a large number of stars over their entire lifetimes. However, it can in principle be implemented in such codes (work is currently under way on an implementation in MESA, Paxton et al. 2013), with the caveat that it will greatly slow down the computations. Possible ways of accelerating the flux calculations are discussed in Section 6.4.

6.3 Implications for stellar abundance observations

With these caveats in mind, we can now apply the complete parasitic saturation model described in Section 5.1 to two situations where fingering convection is known to be important in stars, as introduced in Section 1: just above the hydrogen-burning shell of RGB stars, and just below the surface of white dwarfs undergoing accretion from a companion or a debris disk. Using the stellar models for a typical 1$M_{\odot}$ RGB star just before the luminosity bump, and a typical 0.59$M_{\odot}$, $T_{\rm eff}=11,150$K DA white dwarf (both models are presented in Garaud et al. 2015), we can obtain very rough estimates222Note that these quantities can vary by one or two orders of magnitude in the fingering regions of RGB stars and white dwarfs, so the numbers provided here are indicative only. for the fingering region viscosity, thermal diffusivity, and magnetic diffusivity. The former are computed as described in Garaud et al. (2015), and the latter is computed using the standard Spitzer (1962) formula for non-degenerate gases, which is valid in these regions. We find that for RGB stars typical parameter values are $\mathrm{Pr}=O(10^{-6})$, $\tau=O(10^{-7})$, $\mathrm{Pm}=O(0.1)$, while for the near-surface layers of white dwarfs, $\mathrm{Pr}\simeq\tau=O(10^{-3})$, $\mathrm{Pm}=O(1)$. Figure 7 shows the corresponding model predictions for $D_{C}/\kappa_{C}$, where $D_{C}$ is the turbulent mixing coefficient for composition resulting from the fingering convection (see equation 25). The left panel shows the RGB star model, while the right panel shows the white dwarf model. In both cases, we present two values of $H_{B}$ corresponding to a 100G field (green dots) and a 1000G field (blue dots), respectively (see caption for detail). We also show, for reference, the model predictions from HG19 as dotted lines of the same colors, and the hydrodynamic model of Brown et al. (2013) as a dashed line. We have confirmed a posteriori that in all cases, the model predictions are well within the region validity of the Boussinesq approximation (see caveats described in Sec. 6.2).

We immediately see that for these parameter values, a magnetic field of 100G-1000G is sufficient to increase mixing by fingering convection by many orders of magnitude above the levels expected for the purely hydrodynamic instability alone, especially in the RGB star. This finding is consistent with the conclusions of HG19, and demonstrates that even a relatively weak magnetic field can have an enormous impact on vertical transport by this instability. As a result, we conclude that fingering convection is indeed sufficient to explain the rapid heavy-element surface abundance changes observed in RGB stars around the luminosity bump (Gratton et al. 2000), as originally proposed by Charbonnel & Zahn (2007), as long as a weak magnetic field is also present. In fact, the tight correlation between $D_{C}$ and $H_{B}$ seen in Fig. 7 (where $D_{C}\propto H_{B}$ depends quadratically on the field strength over much of the fingering range) might be used to constrain the interior field strength in the vicinity of the fingering region using the surface abundance observations (cf. Fraser et al. 2022b).

With regards to the white dwarf results, we also find that vertical transport by the fingering instability is greatly enhanced by the field, especially at larger density ratios. This implies that the downward transport by fingering convection of the heavy elements accreted onto the surface from a debris disk (cf. Deal et al. 2013; Bauer & Bildsten 2019) will be much faster in strongly magnetized white dwarfs than in non-magnetized or weakly magnetic ones. As a result, one may expect to see an inverse correlation between the magnetic field strength and the surface abundance of heavy elements in white dwarfs undergoing active accretion. Equivalently, this implies that a much larger accretion rate than previously thought is required to explain surface abundances of heavy elements in magnetic white dwarfs than in non-magnetic ones, as already noted by HG19.

6.4 Implications for stellar evolution models

Figure 7 also shows that the HG19 model (dotted lines) and the new complete parasitic saturation model (large dots) are consistent with one another (except at very large $R_{0}$, see below) within the margin of error of the model333Indeed, the factor of about 2 difference in the low and moderate $R_{0}$ regime is not unexpected given that the new model is itself only accurate within a factor of 2 of the DNS data.. This can easily be explained a posteriori by computing the estimated magnetic Reynolds number of the fingers $\mathrm{Rm}$ (using equation (35) and replacing $\hat{u}_{z,{\rm rms}}$ with $\hat{w}_{f}$) and noting that it is much larger than one at these parameter values, except when $R_{0}\rightarrow\tau^{-1}$ . Hence, the major discrepancies reported in Section 4 between the DNS data and the HG19 model can be understood, in retrospect, as an artefact of using much larger values of $\mathrm{Pr}$ and $\tau$ in the DNS than in real stars. When $\mathrm{Pr}$ and $\tau$ are much smaller, the predicted magnetic Reynolds number $\mathrm{Rm}$ of the fingers is correspondingly much larger, and HG19 generally applies except at very large $R_{0}$.

The complete parasitic model and the HG19 model do differ significantly in the limit $R_{0}\rightarrow\tau^{-1}$, however, despite the fact that $\mathrm{Rm}$ remains large for most of that range (for instance, for the $B=1000$G case for the RGB parameters at $R_{0}=8\times 10^{6}$, where the HG19 model and the new model begin to disagree significantly, the new model predicts $\mathrm{Rm}\approx 700$). Instead, in that limit the parasitic modes take the form of buoyancy-driven modes whose growth rate is amplified by the shear, which were absent in HG19 because their model neglected buoyancy. As a result, HG19 failed to capture the decrease in the amount of turbulent mixing that must inevitably occur as the system becomes stable to fingering convection. The complete parasitic model, on the other hand, correctly predicts that $D_{C}\rightarrow 0$ in that limit. As such, we strongly advocate for its use over HG19 if at all possible.

Yet, as mentioned earlier, the added computational cost is quite substantial, and a simplified parametrization of the new model would be preferable. To create it, one may leverage the knowledge that HG19 is adequate when $\mathrm{Rm}\gg 1$ and $R_{0}\ll\tau^{-1}$, and use it when these conditions are satisfied. The complete model would then only be needed to compute the turbulent fluxes when $\mathrm{Rm}=O(1)$ or less, and/or $R_{0}\rightarrow\tau^{-1}$. Work towards implementing such a model in MESA is currently in progress.

6.5 Other model predictions

There are several other interesting implications of the new model that are worth mentioning as well. First, because the magnetic Reynolds number $\mathrm{Rm}$ is much greater than one in most of the fingering range at stellar values of $\mathrm{Pr}$, $\tau$ and $\mathrm{Pm}$, it is quite plausible that the fluid motions associated with the saturated fingering convection could drive a small-scale dynamo. Whether that small-scale dynamo would then substantially change the model predictions, however, remains to be determined.

Another interesting prospect is that because the magnetic field can greatly enhance the turbulent temperature and compositional fluxes in fingering regions of stars, they may now become unstable to large-scale mean-field instabilities such as the collective instability and the layering instability (see the review by Garaud 2018). Garaud et al. (2015) had shown that the turbulent fingering fluxes are too weak to trigger these instabilities in the hydrodynamic case. Our results suggest that they may actually be active in the magnetized case (although will necessarily also be affected by the large-scale field). This will be the subject of future work.

In the meantime, this study has once more illustrated the power of parasitic saturation models in predicting the nonlinear saturation of a particular instability and the resulting turbulent mixing it causes. The key, as we have demonstrated here, is to take great care in ensuring that all physical processes are included, so as to capture all possible parasitic modes in the calculation – or otherwise run the risk of vastly overestimating or underestimating the level of saturation. Based on their natural similarities, we anticipate that a similar parasitic model could adequately predict the saturation of the magnetized GSF instability, for instance, as well as many other instabilities where saturation occurs because of parasitic modes.

Acknowledgements

AF and PG acknowledge support from National Science Foundation grant AST-1908338. AF acknowledges additional support from NASA HTMS grant 80NSSC20K1280, and from the George Ellery Hale Postdoctoral Fellowship in Solar, Stellar and Space Physics at University of Colorado. Simulations performed using the Extreme Science and Engineering Discovery Environment (XSEDE) Expanse supercomputer at the San Diego Supercomputer Center through allocation TG-PHY210050, and the Lux supercomputer at UC Santa Cruz, funded by NSF MRI grant AST-1828315. We gratefully acknowledge helpful discussions with Matteo Cantiello and Daniel Lecoanet at the KITP program "Probes of transport in stars" (supported in part by the National Science Foundation under Grant No. NSF PHY-1748958), helpful discussions with Justin Walker, David Hughes, Timour Radko, Stephan Stellmach, and Rafa Fuentes, and helpful discussions with Evan Bauer and Rich Townsend who, while working on a separate project, independently verified some of the calculations and figures in this work.