The interaction between shear and fingering (thermohaline) convection

TL;DR

This study investigates how shear flows influence thermohaline convection in stars, revealing that shear reduces mixing efficiency and proposing an analytical model to incorporate these effects into stellar evolution simulations.

Contribution

It is the first to analyze the impact of shear on fingering convection using simulations and to extend existing models to include shear effects.

Findings

Moderate shear decreases fingering convection efficiency.

A simple analytical model explains the effect of shear on mixing.

Turbulent viscosity correlates with compositional mixing.

Abstract

Fingering convection is a turbulent mixing process that can occur in stellar radiative regions whenever the mean molecular weight increases with radius. In some cases, it can have a significant observable impact on stellar structure and evolution. The efficiency of mixing by fingering convection as a standalone process has been studied by Brown et al. (2013), but other processes such as rotation, magnetic fields and shear can affect it. In this paper, we present a first study of the effect of shear on fingering (thermohaline) convection in astrophysics. Using Direct Numerical Simulations we find that a moderate amount of shear (that is not intrinsically shear-unstable) always decreases the mixing efficiency of fingering convection, as a result of the tilt it imparts to the fingering structures. We propose a simple analytical extension of the Brown et al. (2013) model in the presence of shear that satisfactorily explains the numerically-derived turbulent compositional mixing coefficient for moderate shearing rates, and can trivially be implemented in stellar evolution codes. We also measure from the numerical simulations a turbulent viscosity, and find that the latter is strongly tied to the turbulent compositional mixing coefficient. Observational implications and caveats of the model are discussed.

Tables

Table 1. Table 1: Summary of the parameters and main results for simulations with Pr = τ = 0.3 Pr 𝜏 0.3 {\rm Pr}=\tau=0.3 . All measurements are taken as time averages once the system has reached a statistically stationary state. For each pair of parameters ( R 0 , Ri ) subscript 𝑅 0 Ri (R_{0},{\rm Ri}) , the third and fourth columns show the time-average of u ^ rms ​ ( t ) subscript ^ 𝑢 rms 𝑡 \hat{u}_{\rm rms}(t) and w ^ rms ​ ( t ) subscript ^ 𝑤 rms 𝑡 \hat{w}_{\rm rms}(t) (see equation 37 ), and the errors represent a standard deviation around the mean. The fifth column represents the time-average of the maximum value of u ^ ¯ ¯ ^ 𝑢 \overline{\hat{u}} , see text for detail, and the error represents a standard deviation around the mean. The sixth column represents the amplitude of the shear S ^ m subscript ^ 𝑆 𝑚 \hat{S}_{m} , see text for detail; the errorbar is the error on the fitting parameter. The seventh column represents the volume average of | ∇ C ^ | 2 superscript ∇ ^ 𝐶 2 |\nabla\hat{C}|^{2} , and the error is the standard deviation about the mean. Finally, the eighth column represents the turbulent viscosity estimated by fitting the relationship between u ^ ​ w ^ ¯ ¯ ^ 𝑢 ^ 𝑤 \overline{\hat{u}\hat{w}} and S ^ ^ 𝑆 \hat{S} , see text for detail. If the relationship is linear, the errorbar is the error on the fitting parameter. If the relationship is nonlinear, the fit of the mean profile is reported first, and the bracket contain the minimum and maximum values of ν ^ turb subscript ^ 𝜈 turb \hat{\nu}_{\rm turb} measured in the center and wings of the profile (see text for detail).

$R_0$	$\mathrm{Ri}$	${\hat{u}}_{\mathrm{rms}}$	${\hat{w}}_{\mathrm{rms}}$	max$\left(\overline{\hat{u}}\right)$	${\hat{S}}_m$	$⟨{\left|\nabla \hat{C}\right|}^2⟩$	${\hat{\nu }}_{\mathrm{turb}}$
1.5	0.1	2.55 $\pm$ 0.01	1.36 $\pm$ 0.04	3.4 $\pm$ 0.05	0.107 $\pm$ 0.02	13.32 $\pm$ 0.24	2.46 $\pm$ 0.01
1.5	0.3	1.46 $\pm$ 0.006	1.37 $\pm$ 0.01	1.8 $\pm$ 0.04	0.056 $\pm$ 0.02	14.4 $\pm$ 0.29	2.72 $\pm$ 0.01
1.5	3	0.81 $\pm$ 0.01	1.39 $\pm$ 0.01	0.55 $\pm$ 0.04	0.017 $\pm$ 0.02	14.8 $\pm$ 0.28	2.8 $\pm$ 0.02
1.5	30	0.72 $\pm$ 0.007	1.39 $\pm$ 0.01	0.18 $\pm$ 0.04	0.006 $\pm$ 0.02	14.8 $\pm$ 0.30	2.23 $\pm$ 0.08
1.75	0.3	3.00 $\pm$ 0.001	0.809 $\pm$ 0.006	4.17 $\pm$ 0.06	0.13 $\pm$ 0.01	4.43 $\pm$ 0.08	1.18 $\pm$ 0.01
1.75	3	0.912 $\pm$ 0.005	0.983 $\pm$ 0.007	1.10 $\pm$ 0.03	0.035 $\pm$ 0.01	6.29 $\pm$ 0.11	1.49 $\pm$ 0.02
1.75	30	0.525 $\pm$ 0.004	0.99 $\pm$ 0.008	0.33 $\pm$ 0.03	0.01 $\pm$ 0.01	6.41 $\pm$ 0.11	1.48 $\pm$ 0.02
1.75	100	0.489 $\pm$ 0.004	0.99 $\pm$ 0.007	0.24 $\pm$ 0.04	0.007 $\pm$ 0.01	6.37 $\pm$ 0.18	1.52 $\pm$ 0.02
2	0.1	11.46 $\pm$ 0.08	0.27 $\pm$ 0.01	16.3 $\pm$ 0.1	0.51 $\pm$ 0.03	0.37 $\pm$ 0.01	0.40 $\pm$ 0.01
2	0.3	5.81 $\pm$ 0.04	0.30 $\pm$ 0.005	8.2 $\pm$ 0.1	0.25 $\pm$ 0.02	0.60 $\pm$ 0.02	0.54 [0.30,0.68]
2	1	2.80 $\pm$ 0.01	0.505 $\pm$ 0.004	3.9 $\pm$ 0.1	0.12 $\pm$ 0.01	1.61 $\pm$ 0.02	0.65 [0.42,0.92]
2	3	1.74 $\pm$ 0.006	0.64 $\pm$ 0.004	2.42 $\pm$ 0.04	0.076 $\pm$ 0.007	2.41 $\pm$ 0.03	0.58 $\pm$ 0.01
2	10	0.95 $\pm$ 0.01	0.70 $\pm$ 0.005	1.26 $\pm$ 0.03	0.04 $\pm$ 0.006	2.86 $\pm$ 0.04	0.67 $\pm$ 0.01
2	30	0.59 $\pm$ 0.007	0.72 $\pm$ 0.006	0.7 $\pm$ 0.02	0.02 $\pm$ 0.007	2.96 $\pm$ 0.05	0.68 $\pm$ 0.01
2.25	0.1	16.03 $\pm$ 0.03	0.18 $\pm$ 0.01	22.6 $\pm$ 0.17	0.71 $\pm$ 0.05	0.16 $\pm$ 0.01	0.24 $\pm$ 0.01
2.25	0.3	9.04 $\pm$ 0.01	0.17 $\pm$ 0.004	12.8 $\pm$ 0.08	0.40 $\pm$ 0.03	0.18 $\pm$ 0.008	0.27 [0.21,0.28]
2.25	3	2.53 $\pm$ 0.006	0.275 $\pm$ 0.003	3.6 $\pm$ 0.08	0.11 $\pm$ 0.01	0.47 $\pm$ 0.01	0.33 [0.12,0.63]
2.25	30	1.18 $\pm$ 0.004	0.43 $\pm$ 0.003	1.62 $\pm$ 0.04	0.05 $\pm$ 0.005	0.99 $\pm$ 0.01	0.13 [0.07,0.35]
2.25	300	0.46 $\pm$ 0.006	0.49 $\pm$ 0.004	0.57 $\pm$ 0.02	0.02 $\pm$ 0.004	1.25 $\pm$ 0.02	0.09 $\pm$ 0.01
2.5	0.1	19.0 $\pm$ 0.08	0.147 $\pm$ 0.015	26.8 $\pm$ 0.26	0.84 $\pm$ 0.07	0.09 $\pm$ 0.01	0.16 $\pm$ 0.01
2.5	0.3	11.37 $\pm$ 0.03	0.12 $\pm$ 0.005	16.0 $\pm$ 0.12	0.50 $\pm$ 0.03	0.09 $\pm$ 0.007	0.16 $\pm$ 0.01
2.5	3	3.38 $\pm$ 0.008	0.14 $\pm$ 0.003	4.77 $\pm$ 0.04	0.15 $\pm$ 0.007	0.125 $\pm$ 0.005	0.19 [0.06,0.23]
2.5	30	1.30 $\pm$ 0.002	0.23 $\pm$ 0.003	1.83 $\pm$ 0.05	0.06 $\pm$ 0.005	0.28 $\pm$ 0.006	0.1 [-0.05,0.36]
2.5	100	0.96 $\pm$ 0.001	0.27 $\pm$ 0.003	1.34 $\pm$ 0.04	0.04 $\pm$ 0.004	0.38 $\pm$ 0.007	0.03 [-0.09,0.26]
2.5	1000	0.57 $\pm$ 0.001	0.32 $\pm$ 0.003	0.78 $\pm$ 0.02	0.025 $\pm$ 0.003	0.51 $\pm$ 0.009	-0.12 [-0.17,0]
2.5	10000	0.21 $\pm$ 0.001	0.34 $\pm$ 0.003	0.21 $\pm$ 0.02	0.0065 $\pm$ 0.003	0.57 $\pm$ 0.01	-0.18 $\pm$ 0.01

Table 2. Table 2: As for Table 1 , but for Pr = τ = 0.03 Pr 𝜏 0.03 {\rm Pr}=\tau=0.03 . These simulations were run in a shorter domain size, with L x = 200 subscript 𝐿 𝑥 200 L_{x}=200 ( L y = 25 subscript 𝐿 𝑦 25 L_{y}=25 and L z = 200 subscript 𝐿 𝑧 200 L_{z}=200 remain unchanged), and have an effective resolution of ( 384 × 48 × 384 384 48 384 384\times 48\times 384 ) equivalent grid points. A statistically stationary state was achieved in all cases.

$R_0$	$\mathrm{Ri}$	${\hat{u}}_{\mathrm{rms}}$	${\hat{w}}_{\mathrm{rms}}$	max$\left(\overline{\hat{u}}\right)$	${\hat{S}}_m$	$⟨{\left|\nabla \hat{C}\right|}^2⟩$	${\hat{\nu }}_{\mathrm{turb}}$
10	1	0.66 $\pm$ 0.003	0.11 $\pm$ 0.001	0.95 $\pm$ 0.01	0.027 $\pm$ 0.005	0.144 $\pm$ 0.003	0.144 $\pm$ 0.01
10	3	0.328 $\pm$ 0.001	0.136 $\pm$ 0.001	0.45 $\pm$ 0.01	0.015 $\pm$ 0.005	0.207 $\pm$ 0.003	0.17 $\pm$ 0.01
20	3	1.299 $\pm$ 0.002	0.015 $\pm$ 0.0006	1.84 $\pm$ 0.02	0.057 $\pm$ 0.004	0.002 $\pm$ 0.0001	0.021 $\pm$ 0.002
20	10	0.683 $\pm$ 0.001	0.017 $\pm$ 0.0005	0.97 $\pm$ 0.01	0.029 $\pm$ 0.003	0.0023 $\pm$ 0.0001	0.024 [0.015,0.027]

Equations

$${\kappa }_{\mathrm{turb}}=\xi \frac{{\kappa }_T}{R_0-1},$$

$$T(x,y,z,t)$$

$$C(x,y,z,t)$$

$$\frac{{\rho }_0(z)}{{\rho }_m}=1-\alpha T_{0z}z+\beta C_{0z}z$$

$$\frac{\tilde{\rho }}{{\rho }_m}=-\alpha \tilde{T}+\beta \tilde{C},$$

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}$$

$$\nabla \cdot \mathbf{u}$$

$$\frac{\partial \tilde{T}}{\partial t}+\mathbf{u}\cdot \nabla \tilde{T}+w(T_{0z}-T_{\mathrm{ad},z})$$

$$\frac{\partial \tilde{C}}{\partial t}+\mathbf{u}\cdot \nabla \tilde{C}+w{C}_{0z}$$

$${\mathbf{u}}_L(x,y,z)=\frac{F_0}{{\rho }_m\nu k_s^2}\sin (k_{s}z){\mathbf{e}}_x.$$

$$N^2\equiv N_T^2+N_C^2=\alpha g(T_{0z}-T_{\mathrm{ad},z})-\beta g{C}_{0z},$$

$$\mathrm{Ri}=\frac{N^2}{S_L^2}=(\alpha g(T_{0z}-T_{\mathrm{ad},z})-\beta g{C}_{0z}){\left(\frac{{\rho }_m\nu k_s}{F_0}\right)}^2.$$

$$d={\left(\frac{{\kappa }_T\nu }{N_T^2}\right)}^{1/4},$$

$$\frac{\partial \hat{\mathbf{u}}}{\partial t}+\hat{\mathbf{u}}\cdot \nabla \hat{\mathbf{u}}$$

$$\nabla \cdot \hat{\mathbf{u}}$$

$$\frac{\partial \hat{T}}{\partial t}+\hat{\mathbf{u}}\cdot \nabla \hat{T}+\hat{w}$$

$$\frac{\partial \hat{C}}{\partial t}+\hat{\mathbf{u}}\cdot \nabla \hat{C}+R_0^{-1}\hat{w}$$

$$\mathrm{Pr}=\frac{\nu }{{\kappa }_T},\tau =\frac{{\kappa }_C}{{\kappa }_T},$$

$$R_0=\frac{\alpha (T_{0z}-T_{\mathrm{ad},z})}{\beta C_{0z}},{\hat{k}}_s=d{k}_s=\frac{2\pi d}{L_z},{\hat{F}}_0=\frac{F_0}{{\rho }_m}\frac{d^3}{{\kappa }_T^2}.$$

$${\hat{\mathbf{u}}}_L(x,y,z)={\hat{u}}_L\sin ({\hat{k}}_{s}z){\mathbf{e}}_x\text{mbox}where{\hat{u}}_L=\frac{{\hat{F}}_0}{\mathrm{Pr}{\hat{k}}_s^2},\text{mbox}and{\hat{S}}_L=\frac{{\hat{F}}_0}{\mathrm{Pr}{\hat{k}}_s},$$

$$\mathrm{Ri}=\frac{{\mathrm{Pr}}^3{\hat{k}}_s^2}{{\hat{F}}_0^2}(1-R_0^{-1})=\frac{\mathrm{Pr}(1-R_0^{-1})}{{\hat{S}}_L^2}=\frac{{\hat{N}}^2}{{\hat{S}}_L^2}.$$

$$\frac{\partial \tilde{\mathbf{u}}}{\partial t}+{\hat{\mathbf{u}}}_L\cdot \nabla \tilde{\mathbf{u}}+\tilde{\mathbf{u}}\cdot \nabla {\hat{\mathbf{u}}}_L$$

$$\nabla \cdot \tilde{\mathbf{u}}$$

$$\frac{\partial \tilde{T}}{\partial t}+{\hat{\mathbf{u}}}_L\cdot \nabla \tilde{T}+\tilde{w}$$

$$\frac{\partial \tilde{C}}{\partial t}+{\hat{\mathbf{u}}}_L\cdot \nabla \tilde{C}+R_0^{-1}\tilde{w}$$

$$\frac{\partial \tilde{w}}{\partial t}$$

$$\frac{\partial \tilde{T}}{\partial t}+\tilde{w}$$

$$\frac{\partial \tilde{C}}{\partial t}+R_0^{-1}\tilde{w}$$

$$\frac{\partial \tilde{u}}{\partial t}+\tilde{w}{\hat{u}}_L{\hat{k}}_s\cos ({\hat{k}}_{s}z)=\mathrm{Pr}{\nabla }^2\tilde{u},$$

$$\tilde{q}(x,y,z,t)=\sum _n{q}_n\exp (i{\hat{k}}_{x}x+in{\hat{k}}_{s}z+\hat{\lambda }t),$$

$$\hat{\lambda }{u}_n+\frac{{\hat{u}}_L{\hat{k}}_x}{2}(u_{n-1}-u_{n+1})+\frac{{\hat{k}}_s{\hat{u}}_L}{2}(w_{n-1}+w_{n+1})=-i{\hat{k}}_x{p}_n-\mathrm{Pr}(n^2{\hat{k}}_s^2+{\hat{k}}_x^2)u_n,$$

$$\hat{\lambda }{w}_n+\frac{{\hat{u}}_L{\hat{k}}_x}{2}(w_{n-1}-w_{n+1})=-in{\hat{k}}_s{p}_n+\mathrm{Pr}(T_n-C_n)-\mathrm{Pr}(n^2{\hat{k}}_s^2+{\hat{k}}_x^2)w_n,$$

$$\hat{\lambda }{T}_n+\frac{{\hat{u}}_L{\hat{k}}_x}{2}(T_{n-1}-T_{n+1})+w_n=-(n^2{\hat{k}}_s^2+{\hat{k}}_x^2)T_n,$$

$$\hat{\lambda }{C}_n+\frac{{\hat{u}}_L{\hat{k}}_x}{2}(C_{n-1}-C_{n+1})+R_0^{-1}{w}_n=-\tau (n^2{\hat{k}}_s^2+{\hat{k}}_x^2)C_n,$$

$${\hat{k}}_x{u}_n+n{\hat{k}}_s{w}_n=0,$$

$$\hat{\lambda }{w}_n+n{\hat{k}}_x{\hat{k}}_s^2\frac{{\hat{u}}_L}{2K_n^2}((n-2)w_{n-1}-(n+2)w_{n+1})+\frac{{\hat{u}}_L{\hat{k}}_x^3}{2K_n^2}(w_{n-1}-w_{n+1})=\frac{{\hat{k}}_x^2}{K_n^2}\mathrm{Pr}(T_n-C_n)-\mathrm{Pr}{K}_n^2{w}_n$$

$${\hat{u}}_{\mathrm{rms}}(t)=⟨{\hat{u}}^2{⟩}^{1/2}\text{mbox}and{\hat{w}}_{\mathrm{rms}}(t)=⟨{\hat{w}}^2{⟩}^{1/2},$$

$${\hat{F}}_C(t)=\left|⟨\hat{w}\hat{C}⟩\right|\simeq \tau R_0⟨|\nabla \hat{C}{|}^2⟩.$$

$$\frac{D_C}{{\kappa }_C}=\frac{{\kappa }_C+{\kappa }_{C,\mathrm{turb}}}{{\kappa }_C}=\frac{-{\kappa }_C{C}_{0z}+⟨wC⟩}{-{\kappa }_C{C}_{0z}}=1-\frac{R_0}{\tau }⟨\hat{w}\hat{C}⟩=1+\frac{R_0}{\tau }{\hat{F}}_C,$$

Full text

hydrodynamics — instabilities — turbulence — stars:evolution

1 Introduction

The double-diffusive fingering instability (also called thermohaline instability) can be a significant source of turbulent mixing in radiative regions of stars where the mean molecular weight $\mu$ increases with radius (Ulrich, 1972). These so-called inverse $\mu$ gradients can arise in different situations (see the review by Garaud, 2018). For instance, they form whenever a star accretes high $\mu$ material on its surface, either from planetary infall (Vauclair, 2004; Garaud, 2011; Théado & Vauclair, 2012), planetary debris accretion and especially on White Dwarfs (Deal et al., 2013; Wachlin et al., 2017; Bauer & Bildsten, 2018), or mass transfer in binaries (Marks & Sarna, 1998; Chen & Han, 2004; Stancliffe et al., 2007). They can also form as a result of the radiative levitation of high-opacity atomic species such as iron and nickel (Théado et al., 2009; Zemskova et al., 2014; Deal et al., 2016), or from nuclear shell-burning in reactions that decrease the mean molecular weight (Ulrich, 1972; Charbonnel & Zahn, 2007; Denissenkov & Merryfield, 2011), such as 3He + 3He $\rightarrow$ p + p + 4He.

Double-diffusive instabilities in general were first discussed in the oceanographic context by Stommel et al. (1956) and Stern (1960). They occur in fluids whose overall density profile is stable to standard convection (i.e. Ledoux-stable in the stellar context), but where the density itself depends on multiple components that diffuse at different rates. When the rapidly-diffusive component (here, temperature) is stably stratified while the slowly diffusive component (here, the chemical composition) is unstably stratified, as it would be in stellar radiative zones with an inverse $\mu$ gradient, the fingering instability can take place. Fluid parcels displaced from high $\mu$ regions downward rapidly equilibrate thermally with their surroundings (due to the fast thermal diffusion) but retain most of their excess $\mu$ (due to the slow compositional diffusion). As a result, they become denser than the environment and sink further. The nonlinear development of the instability in stellar interiors takes the form of roughly spherical parcels of higher entropy, higher $\mu$ fluid moving down, and lower entropy, lower $\mu$ fluid moving up (Traxler et al., 2011a).

Ulrich (1972) and Kippenhahn et al. (1980) were the first to propose a parametrization for mixing by the fingering instability in the context of stellar astrophysics. While their approaches differ substantially from one another, the mixing coefficient they arrive at is the same aside from a constant multiplicative factor $\xi$:

[TABLE]

where $\kappa_{T}$ is the thermal diffusivity, $R_{0}=\delta(\nabla-\nabla_{\rm ad})/\phi\nabla_{\mu}$ is the density ratio, and where $\delta$, $\phi$, $\nabla$, $\nabla_{\rm ad}$ and $\nabla_{\mu}$ are defined as usual in stellar evolution (e.g. Kippenhahn et al., 1990). Ulrich (1972) argued that the constant $\xi$ would be of the order of 1000, while Kippenhahn et al. (1980) put forward a much smaller value of the order of 10. Nowadays this simple model is commonly used in stellar evolution computations, with the constant $\xi$ left as a tunable parameter (as in Paxton et al., 2011, for instance).

Recent developments in computational fluid dynamics have revitalized the field, and provide means of performing direct numerical simulations (DNSs) to study the nonlinear behavior of the fingering instability in detail. Early works along those lines (Denissenkov, 2010; Traxler et al., 2011a; Denissenkov & Merryfield, 2011) focussed on testing the validity of the turbulent mixing prescription (1), and unanimously ruled in favor of the smaller constant proposed by Kippenhahn et al. (1980). The more systematic numerical experiments of Traxler et al. (2011a) and Brown et al. (2013) however revealed an additional dependence of the mixing coefficient on the Schmidt number (i.e. the ratio of the kinematic viscosity $\nu$ to the compositional diffusivity $\kappa_{C}$), that had so far been ignored. To explain their results, Brown et al. (2013) (see also Radko & Smith, 2012) proposed that the exponential growth of the fingering instability stalls when parasitic shear instabilities between adjacent fingers cause their nonlinear saturation. They assumed that this happens when the growth rate of these parasitic instabilities is of the order of the growth rate of the basic fingering instability (see Section 5 for more detail). Using this idea, they were able to quantify the nonlinear saturation of the fingering instability and to better account for the variation of the turbulent mixing coefficient with input parameters measured in the DNSs.

More recently, Garaud et al. (2015b) looked into the geophysically-inspired notion that the basic fingering instability might excite larger-scale dynamics such as internal gravity waves or thermo-compositional staircases in stars, but ultimately showed that these processes do not seem to take place at stellar parameters, at least in the absence of any other physical process at play. This naturally raised the follow-up question of whether (and how) processes such as rotation, horizontal temperature and compositional gradients, magnetic fields or shear would affect the dynamics of fingering convection. While most of the subsequent research remains quite preliminary, interesting discoveries have been made. Medrano et al. (2014) (see also Holyer, 1983) for instance showed that horizontal gradients greatly widen the region of parameter space that is linearly unstable, to the extent that any inverse $\mu$ gradient, however weak, could drive so-called lateral intrusions. However, the turbulent mixing coefficient in the regime stable to fingering convection but unstable to lateral intrusions remains weak. Sengupta & Garaud (2018) considered fingering convection in a rotating fluid, with the assumption that the rotation axis is aligned with gravity (i.e. in the polar regions of a star). They found that rotation generally reduces the efficiency of turbulent mixing, except in some specific regions of parameter space where large-scale vortices aligned with the rotation axis form. The presence of these vortices appears to greatly enhance vertical mixing by fingering convection, by many orders of magnitude above what one would observe in the non-rotating case. Whether such vortices are actually present in stars, however, is unknown. Finally, Harrington & Garaud (2019) looked at the effect of large-scale magnetic fields. They found that when the field is aligned with gravity (and consequently with the fastest-growing fingering modes), it stabilizes the fingers against their parasitic instabilities and allows them to grow to much larger amplitudes before saturating. This greatly enhances transport by fingering convection compared with the non-magnetic case.

An important process that remains to be investigated is the effect of shear on fingering instabilities. Shear is always present in stars, albeit with a wide range of amplitudes, lengthscales and timescales. Global-scale slowly varying shear flows are observed in stars thanks to helio- and asteroseismology. For instance the solar tachocline is a thin rotational shear layer located at the base of the convection zone, whose radial shear is about 10% of the mean rotation rate of the Sun (Christensen-Dalsgaard & Schou, 1988; Brown et al., 1989; Charbonneau et al., 1999). The envelopes of subgiants are observed to rotate significantly slower than their cores, demonstrating the presence of (presumably) large-scale rotational shear as well (Beck et al., 2012; Deheuvels et al., 2014). Smaller-scale and more rapidly varying shear flows are also very likely in stars. They could arise for instance from the interaction between gravity waves and rotation (in a manner similar to the quasi-biennal oscillation in the Earth’s atmosphere), as discussed by Kumar et al. (1999) in the solar context. They could also be driven by tides in close binary stars (see the review by Ogilvie, 2014), as well as a large number of other possible mechanisms. Since the size of fingering structures is typically of the order of meters to hundreds of meters at stellar parameters, even a small-scale shear flow from the perspective of the entire star would appear to be large-scale from the perspective of the fingers themselves. We can therefore ask the question of how a “large-scale” shear affects fingering convection in stars, and in particular, its associated turbulent mixing properties.

The related question of how fingering convection transports momentum and influences shear flows is also of interest, especially in the context of the aforementioned observations of rotational shear in subgiant stars. Indeed, these observations are at odds with predictions from stellar evolution calculations in which angular momentum is conserved by each spherical mass shell as the core contracts and the envelope expands. Such calculations predict a much higher level of differential rotation than what is observed (Marques et al., 2013). To reconcile theory and observations requires efficient angular momentum transport between the core and the envelope as the star enters and evolves along the subgiant branch. Fingering convection has naturally been proposed as a possible mechanism for turbulent momentum transport in these objects, but little is known about its efficiency. The latter is usually assumed to be proportional to the turbulence compositional mixing coefficient, with the proportionality constant left as a tunable parameter again (cf. Paxton et al., 2011). Whether this prescription is reasonable or not remains to be determined.

In what follows we therefore undertake a preliminary set of numerical experiments to quantify mixing of both chemical species and momentum in sheared fingering convection in stellar interiors. To our knowledge, this has never been studied before, though related work on the effect of shear on geophysical fingering convection exists (Linden, 1974; Ruddick, 1985; Kunze, 1990, 1994; Wells et al., 2001; Smyth & Kimura, 2007; Kimura & Smyth, 2007; Radko et al., 2015, see below for a discussion of their results). Section 2 presents our mathematical model, and Section 3 analyses its linear stability properties. In Section 4 we then describe and analyse a series of DNSs spanning a wide region of parameter space and discuss the emergence of two distinct regimes: a fingering-dominated regime, and a shear-dominated regime. Section 5 presents a very simple geometric model to explain the various quantitative trends observed in our simulations, and Section 6 concludes by discussing astrophysical implications of our results, and future work.

2 Model setup

Given that our goal is to study the effect of shear on fingering instabilities in the context of stellar interiors where there are no solid boundaries, we must adopt a model setup that is also free of solid boundaries. To do so, it is common to use a computational domain that is periodic in all directions. In order to add shear to this triply-periodic system, there are two possibilities: the first is to use a shearing box setup, which assumes the existence of a fixed uniform shear throughout the domain, and the second is a body-forced setup, in which the shear flow is driven by a steady spatially-periodic force. We adopt the body-forced setup primarily because it enables us to study the two-way interaction between fingering and shear, and in particular, how the turbulent flux of momentum induced by fingering convection controls the evolution of the large-scale shear flow (which cannot be done in the shearing box model). As we demonstrate in Section 5, however, there are also drawbacks to this setup, so further work using a shearing box will be presented in a follow-up publication.

In all that follows, we model a small region of a star using a local Cartesian domain $(x,y,z)$ with gravity aligned with the vertical axis: ${\bf g}=-g{\bf e}_{z}$. We ignore the effects of rotation and magnetic fields for now, and refer the reader to the works of Sengupta & Garaud (2018) and Harrington & Garaud (2019) for preliminary results on the topic. We use the Boussinesq approximation for weakly compressible gases (Spiegel & Veronis, 1960), which is valid as long as the height of the computational domain $L_{z}$ is much smaller than any density or temperature scaleheight, and the typical turbulent velocities are much smaller than the local sound speed. Both conditions are satisfied deep in stellar interiors, primarily because double-diffusive turbulence at stellar parameters is both weak and small scale (Garaud, 2018). However, they may fail very close to the surface where both the pressure scaleheight and the sound speed drop significantly, and the double-diffusive lengthscale increases. As such, we restrict the scope of this study to double-diffusive regions that are well below the surface. Consistent with the Boussinesq approximation, we also take $g$ and all the diffusivities ($\nu$, $\kappa_{T}$ and $\kappa_{C}$) to be constant within the domain.

We then write the temperature and composition as the sum of a linear background stratification plus perturbations that are triply periodic:

[TABLE]

where $T_{0}(z)=T_{m}+T_{0z}z$ and $C_{0}(z)=C_{m}+C_{0z}z$, where $T_{m}$, $C_{m}$, $T_{0z}$ and $C_{0z}$ are constant. The compositional field $C$ can be viewed as the concentration of a particular species per unit mass, or as the mean molecular weight $\mu$. As required by the Boussinesq approximation, we linearize the equation of state around the mean state denoted by the subscript $m$ so the mean density profile $\rho_{0}(z)$ and perturbations $\tilde{\rho}$ are given by

[TABLE]

where $\rho_{m}=\rho(p_{m},T_{m},C_{m})$ is the mean density of the region, $p_{m}$ is the mean pressure, and where $\alpha=-\rho_{m}^{-1}\frac{\partial\rho}{\partial T}$, and $\beta=\rho_{m}^{-1}\frac{\partial\rho}{\partial C}$ are the corresponding partial derivatives of the equation of state at $T_{m}$, $C_{m}$, and $p_{m}$.

The Boussinesq equations governing the fluid evolution under the effect of a body-force ${\bf F}$ are:

[TABLE]

where ${\bf u}=(u,v,w)$, $T_{{\rm ad},z}=-g/c_{p}$ is the adiabatic temperature gradient, and $c_{p}$ is the specific heat at constant pressure.

In all that follows, the shear will be driven by a sinusoidal body force ${\bf F}=F_{0}\sin(k_{s}z){\bf e}_{x}$ where $k_{s}=2\pi/L_{z}$ is the wavenumber associated with the domain height. With this assumption, the laminar steady state solution of the governing equations is

[TABLE]

For this laminar flow, the maximum value of the shear is $|S_{L}|=\frac{F_{0}}{\rho_{m}\nu k_{s}}$, and is achieved at $z=0$, $z=L_{z}/2$, and $z=L_{z}$. The constant buoyancy frequency in this system is

[TABLE]

where $N_{T}^{2}=\alpha g(T_{0z}-T_{{\rm ad},z})$ and $N_{C}^{2}=-\beta gC_{0z}$ are the contributions from the temperature and compositional stratifications, respectively. With this information we can compute the Richardson number of the laminar flow as

[TABLE]

It is customary in studies of double-diffusive convection to non-dimensionalize the variables and equations using the anticipated finger scale,

[TABLE]

as the unit length, the thermal diffusion timescale across $d$ (namely $d^{2}/\kappa_{T}$) as the unit time, $d(T_{0z}-T_{{\rm ad},z})$ as the unit temperature and $d\alpha(T_{0z}-T_{{\rm ad},z})/\beta$ as the unit composition. The non-dimensional equations are then

[TABLE]

where all the hatted quantities are from here on non-dimensional111To simplify the notation, we have not added hats on the independent variables $x,y,z$ and $t$, or on the differential operators; their non-dimensionalization is implicit. and

[TABLE]

The Prandtl number ${\rm Pr}$ is typically very small in stellar interiors, ranging from values of the order of $10^{-2}$ in degenerate regions of White Dwarfs, to values of the order of $10^{-6}$ or even smaller in non-degenerate regions of Red Giants and Main Sequence stars. The diffusivity ratio $\tau$ is also very small in stars, and is typically a few times to an order of magnitude smaller than the Prandtl number. The density ratio $R_{0}$ is the ratio of the background density gradient due to the entropy (or equivalently potential temperature) stratification, to the background density gradient due to the compositional stratification. When no shear is present, a fingering instability developing on linear background profiles of $T$ and $C$ is known to take place for values of $R_{0}$ in the interval $[1,\tau^{-1}]$ (Baines & Gill, 1969) with $R_{0}=1$ corresponding to the Ledoux criterion (so $R_{0}<1$ regions are convectively unstable), and $R_{0}=\tau^{-1}$ corresponding to marginal stability to fingering convection, see Garaud (2018) for a review.

The non-dimensional laminar flow and shear amplitudes are, respectively,

[TABLE]

so the Richardson number defined in equation (12) can also be expressed as

[TABLE]

which defines $\hat{N}^{2}$. This expression is used to select the forcing amplitude $\hat{F}_{0}$ for each simulation, so that we can achieve a desired laminar Richardson number ${\rm Ri}$ for any given Prandtl number, density ratio and domain size.

The stability properties of sinusoidal shear flows in the limit of infinite density ratio (i.e. in the absence of compositional stratification) have been studied extensively. When both viscosity and thermal diffusion are neglected (i.e. by artificially removing the Laplacian in the momentum and temperature equations), the shear is linearly unstable provided ${\rm Ri}<1/4$, and linearly stable if ${\rm Ri}>1/4$. This is in accordance with the Howard-Miles theorem for non-diffusive shear flows (Miles, 1961; Howard, 1961). When viscosity is included, the critical Richardson number can increase to a value of 1, thanks to the emergence of a long-wavelength oscillatory viscous instability (Balmforth & Young, 2002, 2005; Garaud et al., 2015a). When thermal diffusion is also included then the stabilizing effects of the temperature stratification are reduced. In that case the critical Richardson number below which instability exists can increase dramatically (see the works of Zahn, 1974, 1992; Prat & Lignières, 2013, 2014; Garaud et al., 2015a; Garaud & Kulenthirarajah, 2016; Prat et al., 2016; Garaud et al., 2017), and the system further becomes nonlinearly unstable to small-scale perturbations. An empirical criterion for nonlinear instability found independently by Prat & Lignières (2014) and Garaud et al. (2017) suggests that turbulence can be sustained in a region with local shear $S$ and local buoyancy frequency $N_{T}$ provided $J_{T}{\rm Pr}<(J{\rm Pr})_{c}\simeq 0.007$, where $J_{T}=N_{T}^{2}/S^{2}$ is the gradient Richardson number associated with the temperature stratification.

By contrast, the effect of shear on a double-diffusively stratified system has been given less attention in general, and so far none in the astrophysical context. Theoretical and experimental studies in the geophysical context, i.e. with Prandtl number of order unity, have shown that shear can stabilize fluid motions that vary along the direction of the mean flow (hereafter called the streamwise direction), but does not affect fluid motions that are streamwise invariant. As a result, so-called salt sheets are expected instead of salt fingers, and have indeed been observed experimentally (Linden, 1974; Kunze, 1990, 1994; Wells et al., 2001; Smyth & Kimura, 2007; Kimura & Smyth, 2007; Smyth & Kimura, 2011; Radko et al., 2015). As summarized by Smyth & Kimura (2011), these salt-sheets are quite resilient to the presence of a moderate steady shear. The main effect of the shear is to decrease the vertical compositional transport despite the fact that kinetic energy is added into the system. Whether this result continues to hold at low Prandtl numbers appropriate for astrophysical environments is one of the questions we aim to address in this paper.

In the following section, we now evaluate the stability of the doubly-stratified sinusoidal shear flow $\hat{\bf u}_{L}$ described above (see equation 20) to infinitesimal perturbations, and then study its nonlinear evolution in Section 4.

3 Linear stability analysis

3.1 Method

Assuming that $\hat{\bf u}=\hat{\bf u}_{L}+\tilde{\bf u}$, where $|\tilde{\bf u}|\ll|\hat{\bf u}_{L}|$, and assuming that $\hat{T}\equiv\tilde{T}\ll 1$ (and similarly for the compositional perturbations), we can linearize the governing equations (14)-(17) as

[TABLE]

This is a set of linear partial differential equations with coefficients that are constant in time and in the horizontal coordinates, but that vary sinusoidally with height $z$ (via $\hat{\bf u}_{L}$).

If we consider perturbations that are invariant both in the streamwise and vertical directions (so $\partial/\partial x$ and $\partial/\partial z$ are null), then we must also have $\tilde{v}=0$ by continuity. The linear stability of these perturbations is identical to that of the corresponding unsheared ones since they satisfy

[TABLE]

The only difference is in the evolution of $\tilde{u}$, which would merely decay viscously in the unsheared case, but satisfies

[TABLE]

when the background shear is present. Since $\tilde{u}$ does not enter in the (linear) evolution of $\tilde{w}$, $\tilde{T}$ and $\tilde{C}$, it does not affect the linear stability of the problem. The solutions to the system of equations (26)-(29) are well known and describe the dynamics of the aforementioned salt sheets (see Section 2, and e.g. Linden (1974)). For Richardson numbers above the threshold for shear instabilities, they are the fastest growing unstable modes.

On the other hand, in the more general situation where $\partial/\partial x$ is non-zero, then the effect of the shear cannot be ignored. To compute the growth rate of such modes, we assume an ansatz of the form

[TABLE]

for $\tilde{q}\in\left\{\tilde{\bf u},\tilde{T},\tilde{C},\tilde{p}\right\}$. This now implicitly assumes that the perturbations are two-dimensional by being invariant in the spanwise ($y$) direction, and that they are periodic with the same period in the vertical direction as the background flow. These assumptions are made for simplicity. Substituting this ansatz into the system of equations (22)-(25), and then projecting the result on individual Fourier modes, we obtain a linear system of coupled eigenvalue equations for the coefficients $q_{n}$, namely

[TABLE]

for any integer index $n$. To solve this numerically, the Fourier expansion needs to be truncated, so $-N\leq n\leq N$. The variables $p_{n}$ and $u_{n}$ can easily be eliminated, so

[TABLE]

where $K_{n}^{2}=n^{2}\hat{k}_{s}^{2}+\hat{k}_{x}^{2}$.

With this, the system can be cast in the form of an eigenvalue problem $\bf A\bf v=\hat{\lambda}\bf v$ where ${\bf v}=(w_{-N},\dots,w_{N},T_{-N},\dots,T_{N},C_{-N},\dots,C_{N})$, and solved numerically using the LAPACK routine DGEEV. For each set of input parameters (${\rm Pr},\tau;{\rm Ri},R_{0};\hat{k}_{s},\hat{k}_{x}$), this procedure yields $3\times(2N+1)$ eigenvalues, from which we select only the one with the largest real part. Finally, we maximize this growth rate over all possible horizontal wavenumbers $\hat{k}_{x}$ to find the growth rate of the fastest-growing (or most slowly decaying) mode for a given physical system. Note that it is quite important to verify that the growth rates found using this method are independent of $N$, the truncation order of the system. This has been checked for all the results presented below.

3.2 Results

Figure 1 summarizes the linear stability analysis of sheared doubly-stratified systems in the fingering regime to spanwise-invariant periodic perturbations, and presents the results in the form of stability diagrams in the $({\rm Ri},R_{0})$ plane for fixed values of ${\rm Pr}$, $\tau$ and $\hat{k}_{s}$ (indicated on the bottom right corner). In each panel, the colors represent the logarithm (in base 10) of the growth rate of the fastest-growing mode found using the method described in the previous section, with the eigensystem (33)-(34) and (36) truncated at $N=100$ (we have verified that this sufficiently large for the results to be independent of $N$). A white color is used when all the modes have negative growth rates, implying that the flow is linearly stable. In the top two panels, the Prandtl number Pr and diffusivity ratio $\tau$ are both equal to 0.3. This is the value used in the majority of the numerical simulations presented in Section 4. Two values of $\hat{k}_{s}$ are presented, for comparison. In the third panel Pr and $\tau$ are both equal to $10^{-3}$, a value that is much lower (though not as low as what one may expect in non-degenerate regions of stellar interiors). These three panels illustrate a number of interesting trends that actually hold for all parameters.

Generally speaking, we find that sheared doubly-stratified systems in the fingering regime can be classified into three categories. For sufficiently small values of the Richardson number ${\rm Ri}$, we find that the system is linearly unstable to shearing modes. These are characterized by a large growth rate commensurate with the flow shearing rate $\hat{S}_{L}$, and a small wavenumber $\hat{k}_{x}$ of the order of the wavenumber of the imposed shear $\hat{k}_{s}$. This mode of instability is relatively independent of the density ratio $R_{0}$ associated with the thermocompositional stratification. For sufficiently large values of the Richardson number, on the other hand, the effect of the shear on the intrinsic fingering instability is negligible, and the modes of instability are fingering modes, characterized by a horizontal wavenumber of the order of unity (in the non-dimensionalization selected) and a growth rate that depends sensitively on $R_{0}$, Pr and $\tau$. We recover the well-known result that the system is only unstable in that limit when $R_{0}<\tau^{-1}$.

Finally, for intermediate values of the Richardson number we see that the system is stabilized (to the assumed infinitesimally small $y-$invariant 2D perturbations) for values of the density ratio smaller than $\tau^{-1}$, showing that moderate shear has a stabilizing effect on the system. The shape of the marginal stability curve depends on the values of Pr, $\tau$ and $\hat{k}_{s}$ selected. In particular we see that shear has a more strongly stabilizing effect at lower values of Pr and $\tau$ (for fixed $\hat{k}_{s}$). To understand qualitatively why this may be the case, recall that the intrinsic growth rate $\hat{\lambda}_{f}$ of the pure fingering instability (i.e in the absence of shear) is much smaller at low Prandtl number than at large Prandtl number (e.g. Brown et al., 2013). As such, the same shearing rate is expected to have a stronger effect on the fingers at low Pr than at high Pr.

In fact, we found empirically that the ratio $\hat{\lambda}_{f}L_{z}/\hat{S}_{L}d$ (where $d=1$ in the units selected) plays an important role in the linear instability of sheared fingering convection. To see this we have overlaid the level contour $\hat{\lambda}_{f}L_{z}/\hat{S}_{L}d=4$ on the growth rate plots (the constant 4 was fitted to the data). It appears that this contour is a good predictor for the edge of the unstable region (at least approximately) in all cases. To understand the significance of this parameter, note that $\hat{\lambda}_{f}L_{z}/\hat{S}_{L}d$ is equal to the ratio of the vertical to horizontal changes in finger velocities (in the linear regime) due to the combined effects of the exponential growth of vertical motions from the fingering instability, and the growth of horizontal motions due to advection by the shear. This suggests that shear can stabilize the fingering instability (at least from a linear perspective) as soon as the shearing rate $\hat{S}_{L}$ exceeds $(L_{z}/4d)\hat{\lambda}_{f}$.

It is worth remembering, however, that all of these linear stability results were obtained for 2D perturbations that are invariant in $y$, and that streamwise-invariant modes are by contrast unaffected by the shear. In addition, it is also important to note that in stars, ${\rm Ri}$ is unknown since it is the Richardson number associated with a hypothetical laminar shear, rather than the one associated with the actual shear in the system. For both reasons, this linear stability analysis is of limited practical interest, and we now proceed to present the results of nonlinear DNSs of the problem.

4 Numerical simulations

4.1 Methodology

In this section, we use DNSs to investigate the nonlinear evolution of the system until it reaches a statistically stationary state, as well as the properties of that state. We use the pseudo-spectral PADDI code (Traxler et al., 2011b; Stellmach et al., 2011), which solves equations (14)-(17) using the Patterson-Orszag algorithm. Linear terms in the equations are advanced implicitly in spectral space, while nonlinear terms are computed first in real space, then transformed back into spectral space before advancing them explicitly using a third order Adams-Bashforth backward differencing algorithm. The code has been extensively used in the astrophysical literature to study both double-diffusive instabilities (e.g. Traxler et al., 2011a; Rosenblum et al., 2011; Mirouh et al., 2012; Brown et al., 2013; Wood et al., 2013) and shear instabilities (e.g. Garaud et al., 2015a; Garaud & Kulenthirarajah, 2016; Gagnier & Garaud, 2018).

The computational domain is triply-periodic, with overall dimensions ($L_{x},L_{y},L_{z}$). Our choice for $L_{x}$, $L_{y}$ and $L_{z}$ is essentially dictated by computational constraints. Since in stars the scale of the shear is expected to be much larger than the scale of individual fingers, we must respect that hierarchy here as well. As shown by Traxler et al. (2011a) the vertical scale of fingers is close to their horizontal scale $2\pi/\hat{l}_{f}$ where $\hat{l}_{f}$ is the wavenumber of the fastest-growing fingering mode (in the absence of shear). Since $\hat{l}_{f}\simeq 0.5$ in our selected units across almost all parameter space, we have to satisfy the condition $L_{z}\gg 2\pi/\hat{l}_{f}\sim 10$. For all the simulations shown below we have picked $L_{z}=200$. When the forcing is strong (low Ri) and the density ratio is relatively high, the fingers are strongly stretched horizontally. In order to avoid being overly constrained by the domain size, most of our simulations are therefore run in long domains with $L_{x}=500$ except where specifically mentioned. With such large values of $L_{x}$ and $L_{z}$, we have then been forced to use a relatively narrow domain in the remaining direction, choosing $L_{y}=25$. As discussed by Garaud & Brummell (2015), the choice of $L_{y}$ in the fingering regime does not affect the results much as long as $L_{y}$ is larger than about 2 finger widths (hence our choice of $L_{y}$). Tables 1 and 2 summarize all the runs we have performed.

Except when specifically mentioned, all the simulations were run until a statistically stationary state was reached, either starting from the laminar state with small added perturbations in the compositional field, or starting from the end-state of another simulation at nearby parameters (e.g. gradually increasing or decreasing either the Richardson number or the density ratio, as appropriate). We have found that the initial conditions used have no influence on the nature of the statistically stationary state reached by the simulation. Owing to the high computational costs of low Prandtl number / low diffusivity ratio simulations (where the viscous and compositional boundary layers at the edges of the fingers are very thin, and yet must be fully resolved), most of the simulations were run with ${\rm Pr}=\tau=0.3$. While this is far from stellar values, these parameters are smaller than one, as they would be in stellar interiors. A few additional runs at ${\rm Pr}=\tau=0.03$ are also presented to verify the adequacy of any model or scaling law we may derive from analyzing the ${\rm Pr}=\tau=0.3$ data. The numerical resolution used for each run (written in terms of equivalent grid points used) varies depending on the parameters selected, and is reported in each case below. We checked that for all the runs shown, the simulations are fully resolved.

4.2 Results for ${\rm Pr}=\tau=0.3$.

We focus for now on a first set of simulations that have ${\rm Pr}=\tau=0.3$. In all cases the domain size is selected to have $L_{x}=500$, $L_{y}=25$ and $L_{z}=200$, and the simulations have been run until a statistically stationary state is reached. For cases where $R_{0}=1.5$, $1.75$ and $2$, the resolutions used is $(960\times 48\times 192)$ equivalent grid points, while for cases where $R_{0}=2.25$ and 2.5, the resolution used is $(480\times 24\times 96)$ equivalent grid points. A summary of all the runs performed in this manner is presented in Table 1.

Figure 2 shows the temporal evolution of a typical simulation in the shear-stable case, that was started from laminar initial conditions plus small-amplitude perturbations, and has a density ratio $R_{0}=2$ (which is in the middle of the fingering-unstable range for $R_{0}$ at $\tau=0.3$) and ${\rm Ri}=3$. Shown are the rms values of the streamwise and vertical velocities

[TABLE]

where $\langle\cdot\rangle$ denotes a volume average over the computational domain. Also shown is a proxy for the absolute value of the compositional flux, which filters out rapid time variations due to internal gravity waves (c.f. Malkus, 1954; Wood et al., 2013) when the system is in a statistically stationary state:

[TABLE]

Note that $\langle\hat{w}\hat{C}\rangle$ in fingering convection is always negative, and so $\langle\hat{w}\hat{C}\rangle=-\hat{F}_{C}$.

We see two phases of evolution in Figure 2. The early time behavior is dominated by the development of the fingering instability on the background shear flow (at this value of ${\rm Ri}$, the shear instability is not active). We see that $\hat{u}_{\rm rms}(t)$ is initially steady while $\hat{w}_{\rm rms}(t)$ grows exponentially, until the perturbations gain sufficient amplitude to affect the mean flow (around $t=30$). This marks the onset of the second phase, in which the mean streamwise flow and the fingering field interact nonlinearly with one another on a slower timescale. When this happens, the streamwise velocity begins to decrease, as a result of the small-scale fingering which acts as an effective turbulent viscosity (more on this below). Because the shear impedes mixing by fingering convection, proxies for the turbulence (namely $\hat{w}_{\rm rms}$ and $\hat{F}_{C}$) slowly grow with time as the shear decreases. This contest between the background shear and fingering convection eventually leads to a statistically stationary state.

Figure 3 shows snapshots of the compositional perturbations in the statistically stationary state for three different values of the density ratio spanning the fingering instability range (which, at the selected values of Pr and $\tau$ is $R_{0}\in[1,10/3]$), and three different values of the Richardson number. When the density ratio is close to the onset of convection ($R_{0}=1.5)$, the intrinsic growth rate of the fingering instability is large, and the shear does not affect the latter much even for small values of ${\rm Ri}$. By contrast, for a large density ratio closer to marginal stability for the fingering instability ($R_{0}=2.5)$, the shear clearly tilts the much more slowly growing fingers even for large Richardson number.

The effect of the shear on vertical compositional mixing by fingering convection can easily be quantified by measuring the ratio

[TABLE]

where $\hat{F}_{C}$ is the compositional flux defined in equation (38). This defines $D_{C}$ as the effective compositional diffusivity, which is the sum of the microscopic diffusivity $\kappa_{C}$ and the turbulent diffusivity $\kappa_{C,{\rm turb}}=-\langle wC\rangle/C_{0z}$. The ratio $D_{C}/\kappa_{C}$ is often called the compositional Nusselt number. Note that the first three expressions in (39) contain only dimensional quantities, while the last two contain non-dimensional quantities. Figure 4 presents the time average of $D_{C}/\kappa_{C}$, measured once the system has reached a statistically stationary state. We see that decreasing ${\rm Ri}$ (i.e. increasing the shear) always acts to reduce $D_{C}/\kappa_{C}$ (except at very low values of the Richardson number where shear instabilities are finally excited, not shown here). However, we also see that the effect of the shear on the fingers is almost negligible at low density ratios, but becomes very significant at higher density ratios, where the flux is reduced even at large values of ${\rm Ri}$. This is consistent with the visual inspection of the simulation snapshots shown in Figure 3. A simple model for the reduction of the compositional flux in the presence of shear is presented in Section 5.

Turbulent fingering convection does not only transport chemical elements, but also transports momentum. This could be seen for instance in the response of the mean flow to the development of the fingering instability in Figure 2. Figure 5 (top row) shows the horizontally averaged and time-averaged streamwise velocity $\overline{\hat{u}}$ as a function of $z$, for two of the simulations that were shown in Figure 3, namely one with $R_{0}=1.5$ and ${\rm Ri}=0.3$, where the fingering is strong, and one with $R_{0}=2.5$ and ${\rm Ri}=30$, where the fingering is weak. As usual, the time average is taken once the simulation is in a statistically stationary state. Also shown is the corresponding laminar mean flow $\hat{u}_{L}\sin(\hat{k}_{s}z)$ for the same simulation. In both cases, we see that the turbulent mean flow has a similar sinusoidal shape to that of the laminar flow (though perhaps slightly more triangular in the $R_{0}=2.5$ and ${\rm Ri}=30$ case), but has a different amplitude. This difference can be used to estimate the turbulent viscosity of the system. In a statistically stationary state, the horizontal average of the momentum equation is

[TABLE]

Assuming, as it is common to do so, that the Reynolds stress can be written as

[TABLE]

(where $\hat{\nu}_{\rm turb}$ is a non-dimensional turbulent viscosity) then the momentum balance reduces to

[TABLE]

whose solution is

[TABLE]

Comparing this expression with the definition of the laminar flow, we should therefore have

[TABLE]

as long as the assumption made in (41) applies. Equation (44) can be used to estimate $\hat{\nu}_{\rm turb}$ from the data (see below).

However, a direct inspection of the relationship between the mean shear and the Reynolds stress $\overline{\hat{u}\hat{w}}$ reveals that (41) does not always apply. The bottom panels of Figure 5 show the time average of the Reynolds stress $\overline{\hat{u}\hat{w}}(z)$ against $\hat{S}(z)\equiv\partial\overline{\hat{u}}/\partial z$, as a parametric plot with $z$ being the parameter. According to (41) these two quantities should be linearly related to one another, and we see that this is indeed the case for the simulation with $R_{0}=1.5$ and ${\rm Ri}=0.3$. However the $R_{0}=2.5$ and ${\rm Ri}=30$ simulation behaves very differently, and reveals a strongly nonlinear relationship between $\overline{\hat{u}\hat{w}}$ and $\hat{S}$. This nonlinear behavior is in fact common for all simulations with relatively low shear and high density ratios. In extreme cases, such as the $R_{0}=2.5$, ${\rm Ri}=30$ case shown here, fingering convection can even behave anti-diffusively, i.e. with $\overline{\hat{u}\hat{w}}$ locally increasing with $\hat{S}$. This recovers a relatively well-known result first obtained by Holyer (1984) and discussed in various other publications since (e.g. Paparella & Spiegel, 1999; Xie et al., 2019). Anti-diffusive tendencies only happen for relatively high ${\rm Ri}$ cases, however. For higher shearing rates, high density ratio simulations have a $\overline{\hat{u}\hat{w}}$ vs. $\hat{S}$ relationship that is once again both linear and decreasing.

These finding suggest two ways of estimating the turbulent diffusivity: one using the stress-strain relationship (41), and one using the mean flow, from (44). In the latter case, we measure the amplitude of $\overline{\hat{u}}(z)$ by fitting the profile to a sinusoidal function. We then solve (44) for $\hat{\nu}_{\rm turb}$ using that measured amplitude. In the former case, we compute the time average of the Reynolds stress profile $\overline{\hat{u}\hat{w}}(z)$ as well as the vertical derivative of the mean flow $\overline{\hat{u}}$ to compute $\hat{S}(z)$, as was done in Figure 5 (bottom panels). We then fit a linear function of the kind

[TABLE]

to this data. For cases where the $\overline{\hat{u}\hat{w}}$ vs. $\hat{S}$ relationship is linear, we set $b=0$. For cases where the $\overline{\hat{u}\hat{w}}$ vs. $\hat{S}$ relationship is not linear, three fits are carried out. The first one, with $b=0$, is a fit for the whole dataset, and reports a mean $\hat{\nu}_{\rm turb}$. Note that for a few simulations this measurement yields a negative value of $\hat{\nu}_{\rm turb}$. The second one, with $b=0$, is taken for low values of $\hat{S}$ only, and reports a minimum value of $\hat{\nu}_{\rm turb}$, that is in some cases negative as well. The third fit is made with $b\neq 0$ to relate $|\overline{\hat{u}\hat{w}}|$ to $|\hat{S}|$ for large values of $|\hat{S}|$ only, which captures both regions of large positive and negative shear at the same time. This final fit estimates a maximum value of $\hat{\nu}_{\rm turb}$. Examples of the results of such fits are shown in Figure 5 (bottom panels), for both the linear case and the nonlinear case. Table 1 presents the results for each simulation.

Figure 6 compares the two methods of estimating the turbulent diffusivity to one another, with the mean-flow based method (using equation 44) on the horizontal axis, and the stress-strain based method (using equation 45) on the vertical axis. The horizontal errorbar is computed from the measured error on the amplitude of $\overline{\hat{u}}(z)$. The vertical errorbar is computed either from the error on the linear fit for cases where the relationship between $\overline{\hat{u}\hat{w}}$ and $\hat{S}$ is linear (typically that error is smaller than the symbol size), or ranges from the minimum to the maximum value of $\hat{\nu}_{\rm turb}$, with the symbol placed at the mean $\hat{\nu}_{\rm turb}$, for cases where the relationship between $\overline{\hat{u}\hat{w}}$ and $\hat{S}$ is nonlinear (see above). For the low density ratio simulations, $R_{0}=1.5$ and $R_{0}=1.75$, the $\overline{\hat{u}\hat{w}}$ vs. $\hat{S}$ relationship is linear for all values of ${\rm Ri}$, and the two methods of computing $\hat{\nu}_{\rm turb}$ are consistent except for two cases. These are the two simulations with largest Richardson numbers, which have very large horizontal errorbars because the horizontal mean flow is very weak. For these runs the mean-flow based method is inaccurate. For simulations at intermediate ($R_{0}=2$) and high ($R_{0}=2.25$ and $R_{0}=2.5$) density ratio, we begin to see examples of nonlinear relationship between $\overline{\hat{u}\hat{w}}$ and $\hat{S}$ especially at intermediate values of ${\rm Ri}$. Note that again for very large ${\rm Ri}$, the data is difficult to interpret as both $\overline{\hat{u}\hat{w}}$ and $\hat{S}$ are small compared with the fingering noise. Despite the nonlinearity of the stress-strain relationship, however, we see that the values of $\hat{\nu}_{\rm turb}$ computed using the mean flow method compares well with that computed by fitting a linear function to the full $\overline{\hat{u}\hat{w}}$ vs. $\hat{S}$ data (regardless of whether it is actually linear or not). This suggests that the details of the stress-strain relationship are less important than the mean behavior in setting the mean flow amplitude.

Finally, as discussed above we see that for the largest density ratio simulations the mean value of $\hat{\nu}_{\rm turb}$ becomes slightly negative (i.e. the turbulence acts in an anti-diffusive way on average), and would on its own act to enhance the shear rather than suppress it. However, we also see that in all these anti-diffusive cases $|\hat{\nu}_{\rm turb}|<{\rm Pr}=0.3$, so the system overall still behaves diffusively overall (i.e. with ${\rm Pr}+\hat{\nu}_{\rm turb}>0$). In fact, so far we have found that $|\hat{\nu}_{\rm turb}|\ll{\rm Pr}$ in all such runs, so the anti-diffusive effect is in practice negligible. Whether this would continue to be the case in stellar interiors remains to be determined, but is quite likely since the turbulent velocities are expected to scale as $\sqrt{{\rm Pr}/(R_{0}-1)}$ (in the pure fingering case at least) so the Reynolds stresses should scale as ${\rm Pr}/(R_{0}-1)$ (see, e.g. Brown et al., 2013; Sengupta & Garaud, 2018). For large $R_{0}$ approaching marginal stability, which is the region of parameter space where anti-diffusive behavior may be expected (Xie et al., 2019), the Reynolds stress is therefore expected to be very small. For this reason, we now set aside the mathematically interesting but in practice probably irrelevant scenario of anti-diffusive fingering convection.

A more practical result can be obtained by comparing the turbulent viscosity to the turbulent compositional diffusivity computed from the simulations, which is related to the flux $\hat{F}_{C}$ defined earlier in equation (38) via

[TABLE]

Note that the second term in that expression is the ratio of two dimensional quantities, while the third is non-dimensional. Figure 7 shows the time-average of $\hat{\kappa}_{C,{\rm turb}}$ (during the statistically stationary state) against $\hat{\nu}_{\rm turb}$. We see that the data follows some interesting trends, and falls into two categories. Simulations that are clearly fingering-dominated, such as those with low density ratio ($R_{0}=1.5$ and $R_{0}=1.75$), or intermediate density ratio ($R_{0}=2$) but low shear, satisfy the relationship $\hat{\nu}_{\rm turb}\simeq 0.25\hat{\kappa}_{C,{\rm turb}}$ (red line). By contrast, simulations that are clearly shear-dominated, such as those with high density ratio and high shear, satisfy the relationship $\hat{\nu}_{\rm turb}\simeq\hat{\kappa}_{C,{\rm turb}}$ (blue line). For some values of the density ratio (especially $R_{0}=2$), the data spans both limits, and continuously moves from one to the other as the Richardson number decreases (i.e. moving from right to left on the figure). The only data points that do not fall on either lines (or in between them) correspond to simulations with high density ratio and low shearing rates. In almost all cases, these also correspond to parameters where anti-diffusive behavior is observed, so $\hat{\nu}_{\rm turb}$ cannot remain proportional to $\hat{\kappa}_{C,{\rm turb}}$ in that limit222By construction, $\hat{\kappa}_{C,{\rm turb}}$ has to be positive, since $\langle\hat{w}\hat{C}\rangle$ has to be negative, so it cannot remain proportional to $\hat{\nu}_{\rm turb}$ when $\hat{\nu}_{\rm turb}$ changes sign.. Since this behavior is unlikely to persist at stellar values of the Prandtl number, these data points a probably not significant for astrophysical purposes.

These results are important because they relate the momentum transport to the compositional transport. Hence if one can somehow be derived from observations (of, e.g. surface abundances or subsurface velocity profiles), the other can be indirectly inferred.

5 Model

We now propose a simple model to explain the trends observed in the previous section, which extends the work of Brown et al. (2013) in the presence of a large-scale shear. For pedagogical purposes, we first briefly recall the salient properties of the original Brown et al. (2013) model for mixing by fingering convection, and then add the shear.

Given that the saturation of the fingering instability occurs as a result of parasitic shear instabilities that develop between up-flowing and down-flowing fingers, a key ingredient of the Brown et al. (2013) model lies in assuming that this saturation occurs when the growth rate of the parasitic instabilities $\hat{\sigma}$ equals a universal constant times the growth rate of the original fingering instability $\hat{\lambda}_{f}$, i.e. when $\hat{\lambda}_{f}=C_{B}\hat{\sigma}$. Furthermore, from dimensional analysis (or exact computations), it can be shown that $\hat{\sigma}=\eta\hat{w}_{f}\hat{l}_{f}$, where $\eta$ is a known constant (whose value is irrelevant). Hence the model predicts a simple relationship between the vertical velocity within the fingers $\hat{w}_{f}$, their growth rate $\hat{\lambda}_{f}$, and their wavenumber $\hat{l}_{f}$, namely

[TABLE]

This idea was successfully verified by Sengupta & Garaud (2018). Next, using the fact that

[TABLE]

from equation (17) in the absence of shear, Brown et al. (2013) estimated the amplitude of the compositional perturbation $\hat{C}_{f}$ to be

[TABLE]

Finally, by dimensional analysis the compositional flux must be proportional to $\hat{w}_{f}\hat{C}_{f}$, which implies that

[TABLE]

where $K_{B}$ is, again, a universal constant. Note how all these universal constants now conveniently combine into a single one, that can be calibrated against the data. The model was successfully validated by Brown et al. (2013) against DNSs of unsheared fingering convection with the constant $K_{B}/(C_{B}\eta)^{2}\simeq 49$.

In the presence of shear, the model must be modified to account for the fact that the fingers become tilted. This changes both their intrinsic velocity (which acquires a horizontal component) and their wavenumber (which decreases as a result of the tilt), see Figure 8. To include this in the model, we begin with a thought experiment in which the fingers initially develop as they would without the shear (i.e. they grow at the same rate $\hat{\lambda}_{f}$ and have the same wavenumber $\hat{l}_{f}$ as the unsheared fingers), but are then subject to a homogeneous shear flow with constant shearing rate $\hat{S}$, whose effect is to tilt them away from the vertical.

With this assumption, a fluid parcel that would be flowing purely vertically in the absence of shear acquires a horizontal velocity relative to the mean flow given by

[TABLE]

where $dz$ is the distance it has travelled, $\hat{w}_{tf}$ is its vertical velocity, and $dt$ is the travel time. The subscript “$tf$” refers to the tilted finger (by contrast with the subscript “$f$” which refers to the non-sheared vertical finger). We then assume that all parcels travel on average for a time $dt$ before losing coherence through nonlinear effects, and that this time is inversely proportional to the growth rate of the fingers $\hat{\lambda}_{f}$ (the reason for this assumption stems from the work of Brown et al. (2013) described above). With this, we predict that

[TABLE]

where $\chi$ is a universal constant of order unity. Note that up-flowing parcels thus acquire a positive horizontal velocity relative to the mean flow in this model, while downflowing parcels acquire a negative horizontal velocity relative to the mean flow. This relationship (and its underlying assumptions) can in fact be verified from the data. To do so, we compute an estimate of the typical vertical finger velocity $\hat{w}_{tf}$ in equation (52) using the rms vertical velocity within the domain, namely $\hat{w}_{\rm rms}$ defined in equation (37). We also estimate the typical horizontal velocity of parcels $\hat{u}_{tf}$ relative to the background from

[TABLE]

Finally, the shear $\hat{S}$ is estimated by fitting a sinusoidal function to $\partial\bar{u}/\partial z$, and reporting its amplitude (which corresponds to the maximum value of the shear in the flow) $\hat{S}_{m}$. Figure 9 shows $\hat{u}^{\prime}_{\rm rms}/\hat{w}_{\rm rms}$ as a function of $\hat{S}_{m}/\hat{\lambda}_{f}$ for all the simulations of Table 1, and we see that the data collapses onto a single universal curve, with two distinct regimes. For low shearing rates (relative to $\hat{\lambda}_{f}$), $\hat{u}^{\prime}_{\rm rms}/\hat{w}_{\rm rms}\simeq 0.5$ is independent of the shear. This corresponds to the fingering-dominated limit, where any horizontal motion of the fluid is simply due to nonlinear saturation of the fingering instability rather than the tilting of the fingers. In that limit, $\hat{u}^{\prime}_{\rm rms}$ is not a good predictor for $\hat{u}_{f}$, which is expected. On the other hand for larger shearing rates (when $\hat{S}_{m}/\hat{\lambda}_{f}$ is larger than one), then we see that $\hat{u}^{\prime}_{\rm rms}/\hat{w}_{\rm rms}\simeq\hat{S}_{m}/3\hat{\lambda}_{f}$, which confirms that $\hat{u}^{\prime}_{\rm rms}$ is now indeed originating from the shear itself, and that equation (52) applies with $\chi\simeq 1/3$. Note that the value of $\chi$ found here is specific to the sinusoidal shear setup used, because both $\hat{u}^{\prime}_{\rm rms}$ and $\hat{w}_{\rm rms}$ are derived from whole-domain averages, while the shear is not constant within the domain. In other words, using a different model setup such as a shearing box could yield a different value of $\chi$, though we expect it to remain close to unity.

With this information, we can now compute the total velocity $\hat{v}_{tf}$ within a tilted finger (as illustrated in Figure 8):

[TABLE]

As in the unsheared case described earlier, we now assume that saturation of the instability occurs when the growth rate of parasitic shear instabilities between two adjacent fingers is of the order of the growth rate of the fingering instability itself. In the tilted finger $\hat{\sigma}=\eta\hat{v}_{tf}\hat{l}_{tf}$ where $\hat{l}_{tf}$ is the wavenumber of tilted finger, and $\eta$ is the same universal constant as before. Using simple geometrical arguments (see Figure 8), it is relatively easy to show that

[TABLE]

so we can expect saturation when

[TABLE]

which implies

[TABLE]

using equation (47). With this, we can compute the expected compositional perturbation within the tilted finger $\hat{C}_{tf}$, using the modal equation for the evolution of the compositional perturbations:

[TABLE]

Note that we have used $\hat{l}_{f}$ instead of $\hat{l}_{tf}$ in this expression, which is done for simplicity (see below); the results are not strongly affected by this simplification, since the compositional diffusion term is usually very small in stars anyway ($\tau\ll 1$). Finally, we find that the magnitude of the compositional flux in tilted fingers is given by

[TABLE]

using equations (50) and (57). Note how the ratio $\hat{F}_{C}/\hat{F}_{C}(0)$ only depends on $\hat{S}/\hat{\lambda}_{f}$ in this model, a result that is reminiscent of the linear stability analysis performed in Section 3.2.

This prediction can be compared with the data from our DNSs, but the comparison is complicated by the fact that the shear is not constant in the domain, which has a number of consequences described below. If we ignore the problem for now, and simply compare the ratio $\hat{F}_{C}/\hat{F}_{C}(0)$ of the volume-averaged compositional flux measured in sheared simulations to the corresponding flux from unsheared simulations, with $\hat{S}_{m}/\hat{\lambda}_{f}$ (where $\hat{S}_{m}$ is the maximum amplitude of the shear measured from the mean flow, see above), we obtain Figure 10a. In all cases, $\hat{F}_{C}=\tau R_{0}\langle|\nabla\hat{C}|^{2}\rangle$ (see equation 38) is computed once the system has achieved a statistically stationary state, while $\hat{F}_{C}(0)$ is similarly measured from a simulation at the same values of $R_{0}$, Pr, and $\tau$, but without shear. We see that the data collapses onto a single universal curve, demonstrating that the relative effect of shear depends only on $\hat{S}_{m}/\hat{\lambda}_{f}$, as in the model (59). The solid black line is the prediction from (59) with $\chi=1/3$, and we see that it fits the data very well up to $\hat{S}_{m}/\hat{\lambda}_{f}\simeq 3$. Note that no additional constant was fitted to create Figure 10a, since the value of $\chi=1/3$ had already been independently fitted earlier to the data relating $\hat{u}^{\prime}_{\rm rms}/\hat{w}_{\rm rms}$ to $\hat{S}_{m}/\hat{\lambda}_{f}\simeq 3$ in Figure 9.

Beyond $\hat{S}_{m}/\hat{\lambda}_{f}\simeq 3$, we see that the model seriously underestimates transport, and goes to 0 as $(\hat{S}_{m}/\hat{\lambda}_{f})^{-4}$, while the data reaches a plateau with $\hat{F}_{C}/\hat{F}_{C}(0)\simeq 0.15$. The origin of this discrepancy remains to be determined and may (or may not) depend on the model setup used in this paper. Several possibilities come to mind, but none of them satisfactorily solve the problem.

The most obvious explanation for this discrepancy is the fact that we neglected the variation of $\hat{S}$ with height in the domain when producing Figure 10a, and effectively assumed that $\hat{S}=\hat{S}_{m}$ everywhere. This overestimates the amount of shear, and therefore underestimates the averaged turbulent compositional flux (which is larger in regions of weaker shear). To compute a “better” estimate of the compositional flux integrated over the computational domain in our sinusoidally forced simulations, we can approximate $\hat{S}(z)\simeq\hat{S}_{m}\cos(\hat{k}_{s}z)$, and evaluate the integral

[TABLE]

We computed the integral numerically as a function of $\hat{S}_{m}/\hat{\lambda}_{f}$ (with $\chi=1/3$) and the results are shown in the green dashed line on Figure 10a. We see that, as expected from the discussion above, $\hat{F}_{C}$ is larger than in the case where we assumed $\hat{S}=\hat{S}_{m}$ everywhere, but still goes to 0 as $\hat{S}_{m}/\hat{\lambda}_{f}$ increases, and fails to explain the plateau in the data. This is not entirely surprising, however. In order to properly account for the effect of a varying shear, we would also need to account for the fact that the vertically varying compositional flux implies a vertically varying background compositional gradient, so the local density ratio is no longer constant in the domain. As a result, the growth rate of the fingers will also depend on $z$, and the model given in (59) needs to be further modified to include that effect. For these reasons, it appears that a better way of testing (59) in the limit of strong shear would be to use simulations in a shearing box model, where the background shear is held constant by construction; this will be the subject of future work.

Other explanations for the enhanced transport at high shearing rates are also possible. For instance, the shear itself could become unstable to diffusive shear instabilities on the finger scale. While we have checked that standard shear instabilities are never excited in the simulations presented because the gradient Richardson number $J=\hat{N}^{2}/\hat{S}_{m}^{2}={\rm Pr}(1-R_{0}^{-1})/{\hat{S}_{m}^{2}}$ is always significantly larger than one, diffusive shear instabilities could take place, as explained in Section 2 (see, e.g. Garaud & Kulenthirarajah, 2016; Garaud et al., 2017). These are difficult to identify unambiguously, however, since they are subcritical in nature, and it is not clear how to distinguish them beyond doubt from other dynamics. The fact that $\hat{\nu}_{\rm turb}\simeq\hat{\kappa}_{C,{\rm turb}}$ in that limit (see Figure 7 and associated text at the end of Section 4) could be a clue to their presence, since this scaling seems to be an intrinsic property of diffusive shear-induced turbulence (Prat et al., 2016; Garaud et al., 2017). However, if that was indeed the case, it is not clear why the compositional flux $\hat{F}_{C}$ should be equal to a constant fraction of $\hat{F}_{C}(0)$ in that limit. Simulations in a different model setup should first be performed to determine whether they also have $\hat{F}_{C}/\hat{F}_{C}(0)\simeq 0.15$ for $\hat{S}_{m}/\hat{\lambda}_{f}>3$, before proceeding with constructing theories for this regime.

Finally note that for even larger values of $\hat{S}_{m}/\hat{\lambda}_{f}$ (not shown here), $\hat{F}_{C}/\hat{F}_{C}(0)$ increases again once the shear instability takes over the fingering instability. This happens when the gradient Richardson number $J$ drops below a critical threshold, at which point the nature of the perturbations changes entirely. Studying this regime is much more computationally demanding, however, because it requires a domain that is at a minimum the same size as we presently use in the $x-$ and $z-$ directions, but much wider in the $y-$direction, to capture the three-dimensionality of the much larger-scale shear-induced eddies. This may be discussed in a future publication.

Having compared our model for compositional transport by sheared fingering convection against data from simulations that had reached a statistically stationnary state, we may now wonder whether that model also applies when the system is slowly evolving in time, a situation that is more likely to be realized in stellar interiors. Figure 10b shows that it does. The red line shows the temporal evolution of $\hat{F}_{C}(\hat{S}(t))/\hat{F}_{C}(0)$ against $\hat{S}_{m}(t)/\hat{\lambda}_{f}$, in a run with $R_{0}=1.5$ and ${\rm Ri}=0.01$ that was initialized from the statistically stationary state reached by the $R_{0}=1.5$, ${\rm Ri}=0.1$ simulation described in Table 1 (with the same resolution and domain size). The increase in the forcing associated with the reduction of ${\rm Ri}$ causes the mean flow to gradually accelerate, and the mean shear to increase with time. As a result, the compositional flux gradually decreases. We see that the flux law (59) derived above still provides a good estimate for transport in this system at each point in time.

Similarly, to see whether the model applies for other values of the Prandtl number and diffusivity ratio, we have run 4 additional simulations for ${\rm Pr}=\tau=0.03$ (see Table 2). Two of these have a density ratio of $R_{0}=10$, which is relatively low at these parameters, and two have a density ratio of $R_{0}=20$, which is closer to marginal stability. Despite the large differences in the governing parameters, we see that these simulations lie close to the existing data at higher ${\rm Pr}$ and $\tau$, for all values of $\hat{S}_{m}/\hat{\lambda}_{f}$ (smaller and greater than 3). This suggests that, at smaller values of ${\rm Pr}$ and $\tau$, (59) continues to hold for low shearing rates, and that $\hat{F}_{C}/\hat{F}_{C}(0)\simeq 0.15$ continues to be true as well for larger shearing rates. We also see in Figure 7 that $\hat{\nu}_{\rm turb}\simeq 0.25\hat{\kappa}_{C,{\rm turb}}$ for the fingering-dominated $R_{0}=10$ simulations and that $\hat{\nu}_{\rm turb}\simeq\hat{\kappa}_{C,{\rm turb}}$ for the shear-dominated $R_{0}=20$ simulations, so these results also appear to be robust regardless of Pr and $\tau$.

6 Implications for stellar interiors, and future work

Despite the uncertainties remaining in the theory described above, we have established that the effect of shear on fingering convection depends on a single quantity, namely the ratio of the shearing rate to the intrinsic (unsheared) finger growth rate $\hat{S}/\hat{\lambda}_{f}$, or, in terms of dimensional quantities, $S/\lambda_{f}$. We have also established that moderate shear reduces the efficiency of fingering, and that the effect becomes pronounced when $S/\lambda_{f}$ is of order unity, or equivalently, when the shearing rate approaches the fingering growth rate $\lambda_{f}$. The formula for the compositional flux given in equation (59) can be re-written in terms of a turbulent mixing coefficient as

[TABLE]

where $\chi=1/3$ for sinusoidal shear flows (though could be a little different in a different setup), and $\kappa_{C,{\rm turb}}(S=0)$ is the turbulent mixing coefficient in the absence of shear, which can be computed for instance using the model of Brown et al. (2013). The formula appears to be valid up to $S/\lambda_{f}\simeq 3$, beyond which it underestimates $\kappa_{C,{\rm turb}}$. For $S/\lambda_{f}>3$, we find that $\kappa_{C,{\rm turb}}(S)\simeq 0.15\kappa_{C,{\rm turb}}(S=0)$ instead, although this result needs to be confirmed using a different model setup (such as a shearing box, for instance).

We have also studied momentum transport by fingering convection, and have found that the turbulent viscosity and the turbulent mixing coefficient are indeed (almost) proportional to one another, with $\nu_{\rm turb}\simeq 0.25\kappa_{C,{\rm turb}}$ when $S/\lambda_{f}<3$, gradually increasing to $\nu_{\rm turb}\simeq\kappa_{C,{\rm turb}}$ when $S/\lambda_{f}>3$. This higher-shear limit remains to be confirmed using a different model setup, however, as discussed in Section 5.

Note that all of these results have been obtained for moderate shear, prior to the onset of standard shear instabilities (i.e. when $J=N^{2}/S^{2}$ is larger than one). For larger shearing rates, we expect shear instabilities to develop, although this limit has not been investigated yet.

Using this information, we can conjecture on the conditions where shear might substantially reduce the fingering-induced mixing efficiency in stellar interiors. To estimate the dimensional finger growth rate, we use the asymptotic formula of Brown et al. (2013) (see their Appendix B.2) for $1\ll R_{0}\ll\tau^{-1}$ (which is the more likely scenario in stars where $\tau$ is asymptotically small) namely

[TABLE]

where $N^{2}_{T}=\alpha g(T_{0z}-T_{{\rm ad},z})$ is the Brunt-Väisälä frequency associated with the temperature stratification only (see Section 2). With this estimate, we predict that shear may become relevant when $\lambda_{f}\sim S$, which is

[TABLE]

In other words, fingering becomes affected by shear when the gradient Richardson number associated with the temperature stratification only is of the order of (or lower than) the density ratio $R_{0}$. Such conditions, while not necessarily ubiquitous, are nevertheless possible depending on the current phase of evolution and conditions within the star.

In Red Giant Branch stars for instance, a fingering region induced by nuclear reactions that transform 3He into 4He is located between the hydrogen burning shell and the convective envelope (Ulrich, 1972), and is thought to contribute to the transport of both Lithium and CNO cycle by-products between the two regions (Charbonnel & Zahn, 2007). Using the estimates from Denissenkov (2010) for typical parameters appropriate of that region, we have $N_{T}^{2}\simeq 2\times 10^{-4}$s*-2*, and $R_{0}\simeq 2\times 10^{3}$. The local shear on the other hand is unknown, the best available constraints coming from asteroseismic estimates of the core and envelope angular velocities $\Omega_{\rm core}/2\pi\simeq 10^{-6}$Hz and $\Omega_{\rm env}/2\pi\simeq 10^{-7}$Hz (Deheuvels et al., 2014). If we assume the transition takes place smoothly over the entire core (from the center to the base of the convective envelope), then $S\simeq\Delta\Omega\simeq 2\pi\times 10^{-6}$s*-1*, which would imply that $J_{T}\simeq 5\times 10^{6}$. Since this is much larger than $R_{0}$, we conclude that a smooth shear profile in these objects would not have any effect on the fingering efficiency. On the other hand if the transition between the core and envelope rotation rates occurs over a distance $\Delta$ that is a fraction of the core size $r_{\rm core}$, then

[TABLE]

and so the shear could begin to affect the fingering efficiency if $\Delta/r_{\rm core}\simeq 0.01$ or less333unless $R_{0}$ is much larger than estimated by Denissenkov (2010), which seems unlikely. As a result, we predict that the turbulent mixing coefficient $\kappa_{C,{\rm turb}}$ in Red Giant Branch stars (ignoring magnetic fields) can be estimated using the model proposed by Brown et al. (2013) in the absence of shear, and further find that the turbulent viscosity can be computed as $\nu_{\rm turb}\simeq 0.25\kappa_{C,{\rm turb}}$. If magnetic fields are present, however, then the model of Harrington & Garaud (2019) for $\kappa_{C,{\rm turb}}$ should be used instead, although no model presently exists for $\nu_{\rm turb}$ in that case.

In White Dwarfs undergoing accretion of material from a surrounding debris disk, a fingering region is thought to be located just beneath the surface (c.f. Deal et al., 2013; Bauer & Bildsten, 2018, for instance), and could participate in draining the excess metallicity of accreted material inward in addition to the effect of gravitational settling. Assuming a steady-state balance between the rate of debris accretion and the inward subsurface flux of heavy elements, one can in theory predict their expected surface abundances. Conversedly, the observed abundances can in principle be used to infer the accretion rate of debris assuming the subsurface flux is known. Unfortunately this inference is complicated by the fact that the fingering flux depends on many different processes, as discussed in Section 1, which have all been neglected to date. To determine whether shear in particular could substantially reduce the efficiency of fingering convection in White Dwarfs, one would need to know how both the gradient Richardson number $J_{T}$ and the density ratio $R_{0}$ evolve with time after each accretion event. Qualitatively speaking, we anticipate that the density ratio $R_{0}$ should start from a value close to unity at the onset of the accretion event (where the inverse $\mu$ gradient is the largest), and should then gradually increase towards marginal stability ($R_{0}=1/\tau$) as the inward transport of heavy elements proceeds. Accretion is also expected to generate substantial shear (as the accreting material lands on the star at close to Keplerian speed), which would then gradually decrease as a result of the turbulent momentum transport. Because of these processes, both $J_{T}$ and $R_{0}$ are likely to increase with time after accretion, so whether shear is ever important ($J_{T}<R_{0}$) or not will depend on the initial angular momentum and composition of the accreting material relative to that of the star, as well as the efficiency of both compositional and momentum transport. This will need to be established on a case-by-case basis for each star and each accretion scenario considered. Finding substantial reduction of the fingering efficiency by the shear would be an interesting outcome of this exercise, however, since it could more easily help interpret the observed heavy-element abundances of polluted DA White Dwarfs, that can otherwise only be explained by very large accretion rates (Bauer & Bildsten, 2018).

The simulations were performed using the PADDI code kindly provided by S. Stellmach, on the Hyades supercomputer purchased using an NSF MRI grant. P. G. and J. S. acknowledge funding by NSF-AAG 1517927. A.K. is supported by the Dean’s fellowship and Regents’ fellowship from the Baskin School of Engineering at UC Santa Cruz.