Turbulent mixing in the radiative regions of stars is usually either ignored or crudely accounted for in most stellar evolution models. However, there is growing theoretical and observational evidence that such mixing is present and can affect various aspects of a star's life. In this work, we present a first attempt at quantifying mixing by horizontal shear instabilities in stars using Direct Numerical Simulations. The shear is driven by a body force, and rapidly becomes unstable. At saturation, we find that several distinct dynamical regimes exist, depending on the relative importance of stratification and thermal diffusion (viscosity can in principle also matter, but is usually negligible in most stellar applications). In each of the regimes identified, we put forward a certain number of theoretically motivated scaling laws for the turbulent vertical eddy scale, the typical turbulent diffusion coefficient, and the typical amplitude of temperature fluctuations (among other quantities). Based on our findings, we predict that the majority of stars should fall into one of two categories: high P\'eclet number stratified turbulence, and low P\'eclet number stratified turbulence. The latter is presented in detail in a related paper by Cope et al. (2020), while the former is discussed here. Applying our results to the best-known stellar shear layer, namely the solar tachocline, we find that it should lie in the high P\'eclet number stratified turbulence regime, and predict a substantial amount of vertical mixing for temperature, momentum and composition. Taken as is, the new turbulence model predictions are incompatible with the Spiegel & Zahn (1992) model of the solar tachocline. However, we also show that rotation and magnetic fields are likely to affect the turbulence, and need to be taken into account in future studies.

The solar tachocline is a shear layer located at the base of the solar convection zone. The horizontal shear in the tachocline is likely turbulent, and it is often assumed that this turbulence would be strongly anisotropic as a result of the local stratification. What role this turbulence plays in the tachocline dynamics, however, remains to be determined. In particular, it is not clear whether it would result in a turbulent eddy diffusivity, or anti-diffusivity, or something else entirely. In this paper, we present the first direct numerical simulations of turbulence in horizontal shear flows at low Prandtl number, in an idealized model that ignores rotation and magnetic fields. We find that several regimes exist, depending on the relative importance of the stratification, viscosity and thermal diffusivity. Our results suggest that the tachocline is in the stratified turbulence regime, which has very specific properties controlled by a balance between buoyancy, inertia, and thermal diffusion.

We consider the dynamics of a vertically stratified, horizontally-forced Kolmogorov flow. Motivated by astrophysical systems where the Prandtl number is often asymptotically small, our focus is the little-studied limit of high Reynolds number but low P\'eclet number (which is defined to be the product of the Reynolds number and the Prandtl number). Through a linear stability analysis, we demonstrate that the stability of two-dimensional modes to infinitesimal perturbations is independent of the stratification, whilst three-dimensional modes are always unstable in the limit of strong stratification and strong thermal diffusion. The subsequent nonlinear evolution and transition to turbulence is studied numerically using direct numerical simulations. For sufficiently large Reynolds numbers, four distinct dynamical regimes naturally emerge, depending upon the strength of the background stratification. By considering dominant balances in the governing equations, we derive scaling laws for each regime which explain the numerical data.

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.

Convective overshooting in super asymptotic giant branch stars has been suggested to lead to the formation of hybrid white dwarfs with carbon-oxygen cores and oxygen-neon mantles. As the white dwarf cools, this core-mantle configuration becomes convectively unstable and should mix. This mixing has been previously studied using stellar evolution calculations, but these made the approximation that convection did not affect the temperature profile of the mixed region. In this work, we perform direct numerical simulations of an idealized problem representing the core-mantle interface of the hybrid white dwarf. We demonstrate that, while the resulting structure within the convection zone is somewhat different than what is assumed in the stellar evolution calculations, the two approaches yield similar results for the size and growth of the mixed region. These hybrid white dwarfs have been invoked as progenitors of various peculiar thermonuclear supernovae. This lends further support to the idea that if these hybrid white dwarfs form then they should be fully mixed by the time of explosion. These effects should be included in the progenitor evolution in order to more accurately characterize the signatures of these events.

Double-diffusive convection has been well studied in geophysical contexts, but detailed investigations of the regimes characteristic of stellar or planetary interiors have only recently become feasible. Since most astrophysical fluids are electrically conducting, it is possible that magnetic fields play a role in either enhancing or suppressing double-diffusive convection, but to date there have been no numerical investigations of such possibilities. Here we study the effects of a vertical background magnetic field (aligned with the gravitational axis) on the linear stability and nonlinear saturation of fingering (thermohaline) convection, through a combination of theoretical work and direct numerical simulations (DNSs). We find that a vertical magnetic field rigidifies the fingers along the vertical direction which has the remarkable effect of enhancing vertical mixing. We propose a simple analytical model for mixing by magnetized fingering convection, and argue that magnetic effects may help explain discrepancies between theoretical and observed mixing rates in low-mass red giant branch (RGB) stars. Other implications of our findings are also discussed.

We study the dynamics associated with the extension of turbulent convective motions from a convection zone (CZ) into a stable region (RZ) that lies below the latter. For that purpose, we have run a series of three-dimensional direct numerical simulations solving the Navier-Stokes equations under the Boussinesq approximation in a spherical shell geometry. We observe that the overshooting of the turbulent motions into the stably stratified region depends on three different parameters: the relative stability of the RZ, the transition width between the two, and the intensity of the turbulence. In the cases studied, these motions manage to partially alter the thermal stratification and induce thermal mixing, but not so efficiently as to extend the nominal CZ further down into the stable region. We find that the kinetic energy below the convection zone can be modeled by a half-Gaussian profile whose amplitude and width can be predicted a priori for all of our simulations. We examine different dynamical lengthscales related to the depth of the extension of the motions into the RZ, and we find that they all scale remarkably well with a lengthscale that stems from a simple energetic argument. We discuss the implications of our findings for 1D stellar evolution calculations.

This paper continues the systematic investigation of diffusive shear instabilities initiated in Part I of this series. In this work, we primarily focus on quantifying the impact of non-local mixing, which is not taken into account in Zahn's mixing model \citep{Zahn92}. We present the results of direct numerical simulations in a new model setup designed to contain coexisting laminar and turbulent shear layers. As in Part I, we use the Low P\'eclet Number approximation of \citet{Lign1999} to model the evolution of the perturbations. Our main findings are twofold. First, turbulence is not necessarily generated whenever Zahn's nonlinear criterion \citep{Zahn1974} $J{\rm Pr} < (J{\rm Pr})_c$ is satisfied, where $J=N^2/S^2$ is the local gradient Richardson number, ${\rm Pr} = \nu/ \kappa_T$ is the Prandtl number, and $(J{\rm Pr})_c \simeq 0.007$. We have demonstrated that the presence or absence of turbulent mixing in this limit hysteretically depends on the history of the shear layer. Second, Zahn's nonlinear instability criterion only approximately locates the edge of the turbulent layer, and mixing beyond the region where $J{\rm Pr} < (J{\rm Pr})_c$ can also take place in a manner analogous to convective overshoot. We found that the turbulent kinetic energy decays roughly exponentially beyond the edge of the shear-unstable region, on a lengthscale $\delta$ that is directly proportional to the scale of the turbulent eddies, which are themselves of the order of the Zahn scale (see Part I). Our results suggest that mixing by diffusive shear instabilities should be modeled with more care than is currently standard in stellar evolution codes.

We study the effects of rotation on the growth and saturation of the double-diffusive fingering (thermohaline) instability at low Prandtl number. Using direct numerical simulations, we estimate the compositional transport rates as a function of the relevant non-dimensional parameters - the Rossby number, inversely proportional to the rotation rate, and the density ratio which measures the relative thermal and compositional stratifications. Within our explored range of parameters, we generally find rotation to have little effect on vertical transport. However, we also present one exceptional case where a cyclonic large scale vortex (LSV) is observed at low density ratio and fairly low Rossby number. The LSV leads to significant enhancement in the fingering transport rates by concentrating compositionally dense downflows at its core. We argue that the formation of such LSVs could be relevant to solving the missing mixing problem in RGB stars.

Being able to account for the missing mixing in stellar radiative zones is a key step toward a better understanding of stellar evolution. Zahn (1974) argued that thermally diffusive shear-induced turbulence might be responsible for some of this mixing. In Part I and Part II of this series of papers we showed that Zahn's (1974, 1992) mixing model applies when the properties of the turbulence are local. But we also discovered limitations of the model when this locality condition fails, in particular near the edge of a turbulent region. In this paper, we propose a second-order closure model for the transport of momentum and chemical species by shear-induced turbulence in strongly stratified, thermally diffusive environments (the so-called low P\'eclet number limit), which builds upon the work of Garaud \& Ogilvie (2005). Comparison against direct numerical simulations (DNSs) shows that the model is able to predict the vertical profiles of the mean flow and of the stress tensor (including the momentum transport) in diffusive shear flows, often with a reasonably good precision, and at least within a factor of order unity in the worst case scenario. The model is sufficiently simple to be implemented in stellar evolution codes, and all the model constants have been calibrated against DNSs. While significant limitations to its use remain (e.g. it can only be used in the low P\'eclet number, slowly rotating limit), we argue that it is more reliable than most of the astrophysical prescriptions that are used in stellar evolution models today.

We present Direct Numerical Simulations of the transport of heat and heavy elements across a double-diffusive interface or a double-diffusive staircase, in conditions that are close to those one may expect to find near the boundary between the heavy-element rich core and the hydrogen-helium envelope of giant planets such as Jupiter. We find that the non-dimensional ratio of the buoyancy flux associated with heavy element transport to the buoyancy flux associated with heat transport lies roughly between 0.5 and 1, which is much larger than previous estimates derived by analogy with geophysical double-diffusive convection. Using these results in combination with a core-erosion model proposed by Guillot et al. (2004), we find that the entire core of Jupiter would be eroded within less than 1Myr assuming that the core-envelope boundary is composed of a single interface. We also propose an alternative model that is more appropriate in the presence of a well-established double-diffusive staircase, and find that in this limit a large fraction of the core could be preserved. These findings are interesting in the context of Juno's recent results, but call for further modeling efforts to better understand the process of core erosion from first principles.

We examine the dynamics associated with weakly compressible convection in a spherical shell by running 3D direct numerical simulations using the Boussinesq formalism [1]. Motivated by problems in astrophysics, we assume the existence of a finite adiabatic temperature gradient $\nabla T_{\rm{ad}}$ and use mixed boundary conditions for the temperature with fixed flux at the inner boundary and fixed temperature at the outer boundary. This setup is intrinsically more asymmetric than the more standard case of Rayleigh-B\'{e}nard convection in liquids between parallel plates with fixed temperature boundary conditions. Conditions where there is substantial asymmetry can cause a dramatic change in the nature of convection and we demonstrate that this is the case here. The flows can become pressure- rather than buoyancy- dominated leading to anomalous heat transport by upflows. Counter-intuitively, the background temperature gradient $\nabla\bar{T}$ can develop a subadiabatic layer (where $\boldsymbol{g}\cdot\nabla\bar{T}<\boldsymbol{g}\cdot\nabla T_{\rm{ad}}$, where $\boldsymbol{g}$ is gravity) although convection remains vigorous at every point across the shell. This indicates a high degree of non-locality.

Shear-induced turbulence could play a significant role in mixing momentum and chemical species in stellar radiation zones, as discussed by Zahn (1974). In this paper we analyze the results of direct numerical simulations of stratified plane Couette flows, in the limit of rapid thermal diffusion, to measure the turbulent diffusivity and turbulent viscosity as a function of the local shear and the local stratification. We find that the stability criterion proposed by Zahn (1974), namely that the product of the gradient Richardson number and the Prandtl number must be smaller than a critical values $(J\Pr)_c$ for instability, adequately accounts for the transition to turbulence in the flow, with $(J\Pr)_c \simeq 0.007$. This result recovers and confirms the prior findings of Prat et al. (2016). Zahn's model for the turbulent diffusivity and viscosity (Zahn 1992), namely that the mixing coefficient should be proportional to the ratio of the thermal diffusivity to the gradient Richardson number, does not satisfactorily match our numerical data when applied as is. It fails (as expected) in the limit of large stratification where the Richardson number exceeds the aforementioned threshold for instability, but it also fails in the limit of low stratification where the turbulent eddy scale becomes limited by the computational domain size. We propose a revised model for turbulent mixing by diffusive stratified shear instabilities, that now properly accounts for both limits, fits our data satisfactorily, and recovers Zahn's 1992 model in the limit of large Reynolds numbers.

Oscillatory double-diffusive convection (ODDC, more traditionally called semiconvection) is a form of linear double-diffusive instability that occurs in fluids that are unstably stratified in temperature (Schwarzschild unstable), but stably stratified in chemical composition (Ledoux stable). This scenario is thought to be quite common in the interiors of stars and giant planets, and understanding the transport of heat and chemical species by ODDC is of great importance to stellar and planetary evolution models. Fluids unstable to ODDC have a tendency to form convective thermo-compositional layers which significantly enhance the fluxes of temperature and chemical composition compared with microscopic diffusion. Although a number of recent studies have focused on studying properties of both layered and non-layered ODDC, few have addressed how additional physical processes such as global rotation affect its dynamics. In this work we study first how rotation affects the linear stability properties of rotating ODDC. Using direct numerical simulations we then analyze the effect of rotation on properties of layered and non-layered ODDC, and study how the angle of the rotation axis with respect to the direction of gravity affects layering. We find that rotating systems can be broadly grouped into two categories, based on the strength of rotation. Qualitative behavior in the more weakly rotating group is similar to non-rotating ODDC, but strongly rotating systems become dominated by vortices that are invariant in the direction of the rotation vector and strongly influence transport. We find that whenever layers form, rotation always acts to reduce thermal and compositional transport.

This work presents numerical results on the transport of heat and chemical species by shear-induced turbulence in strongly stratified but thermally diffusive environments. The shear instabilities driven in this regime are sometimes called "secular" shear instabilities, and can take place even when the gradient Richardson number of the flow (the square of the ratio of the buoyancy frequency to the shearing rate) is large, provided the P\'eclet number (the ratio of the thermal diffusion timescale to the turnover timescale of the turbulent eddies) is small. We have identified a set of simple criteria to determine whether these instabilities can take place or not. Generally speaking, we find that they may be relevant whenever the thermal diffusivity of the fluid is very large (typically larger than $10^{14}$cm$^2$/s), which is the case in the outer layers of high-mass stars ($M\ge 10 M_\odot$) for instance. Using a simple model setup in which the shear is forced by a spatially sinusoidal, constant-amplitude body-force, we have identified several regimes ranging from effectively unstratified to very strongly stratified, each with its own set of dynamical properties. Unless the system is in one of the two extreme regimes (effectively unstratified or completely stable), however, we find that (1) only about 10% of the input power is used towards heat transport, while the remaining 90% is viscously dissipated; (2) that the effective compositional mixing coefficient is well-approximated by the model of Zahn (1992), with $D \simeq 0.02 \kappa_T /J$ where $\kappa_T$ is the thermal diffusivity and $J$ is the gradient Richardson number. These results need to be confirmed, however, with simulations in different model setups and at higher effective Reynolds number.

Fingering convection (otherwise known as thermohaline convection) is an instability that occurs in stellar radiative interiors in the presence of unstable compositional gradients. Numerical simulations have been used in order to estimate the efficiency of mixing induced by this instability. However, fully three-dimensional (3D) computations in the parameter regime appropriate for stellar astrophysics (i.e. low Prandtl number) are prohibitively expensive. This raises the question of whether two-dimensional (2D) simulations could be used instead to achieve the same goals. In this work, we address this issue by comparing the outcome of 2D and 3D simulations of fingering convection at low Prandtl number. We find that 2D simulations are never appropriate. However, we also find that the required 3D computational domain does not have to be very wide: the third dimension need only contain a minimum of two wavelengths of the fastest-growing linearly unstable mode to capture the essentially 3D dynamics of small-scale fingering. Narrow domains, however, should still be used with caution since they could limit the subsequent development of any large-scale dynamics typically associated with fingering convection.

This work addresses the question of the stability of stratified, spatially periodic shear flows at low P\'eclet number but high Reynolds number. This little-studied limit is motivated by astrophysical systems, where the Prandtl number is often very small. Furthermore, it can be studied using a reduced set of "low-P\'eclet-number equations" proposed by Lignieres [Astronomy & Astrophysics, 348, 933-939, 1999]. Through a linear stability analysis, we first determine the conditions for instability to infinitesimal perturbations. We formally extend Squire's theorem to the low-P\'eclet-number equations, which shows that the first unstable mode is always two-dimensional. We then perform an energy stability analysis of the low-P\'eclet-number equations and prove that for a given value of the Reynolds number, above a critical strength of the stratification, any smooth periodic shear flow is stable to perturbations of arbitrary amplitude. In that parameter regime, the flow can only be laminar and turbulent mixing does not take place. Finding that the conditions for linear and energy stability are different, we thus identify a region in parameter space where finite-amplitude instabilities could exist. Using direct numerical simulations, we indeed find that the system is subject to such finite-amplitude instabilities. We determine numerically how far into the linearly stable region of parameter space turbulence can be sustained.

Oscillatory double-diffusive convection (ODDC) (also known as semi- convection) refers to a type of double diffusive instability that occurs in regions of planetary and stellar interiors which have a destabilizing thermal stratification and a stabilizing mean molecular weight stratification. In this series of papers, we use an extensive suite of three-dimensional (3D) numerical simulations to quantify the transport of heat and chemical species by ODDC. Rosenblum et al. (2011) first showed that ODDC can either spontaneously form layers, which significantly enhance the transport of heat and chemical species compared to mi- croscopic transport, or remain in a state dominated by large scale gravity waves, in which there is a more modest enhancement of the turbulent transport rates. Subsequent studies in this series have focused on identifying under what condi- tions layers form (Mirouh et al. 2012), and quantifying transport through layered systems (Wood et al. 2013). Here we proceed to characterize transport through systems that are unstable to the ODDC instability, but do not undergo spon- taneous layer formation. We measure the thermal and compositional fluxes in non-layered ODDC from both 2D and 3D numerical simulations and show that 3D simulations are well approximated by similar simulations in a 2D domain. We find that the turbulent mixing rate in this regime is weak and can, to a first level approximation, be neglected. We conclude by summarizing the findings of papers I through III into a single prescription for transport by ODDC.

Semiconvection - mixing that occurs in regions that are stable when considering compositional gradients, but unstable when ignoring them - is shown to have the greatest potential impact on main sequence stars with masses in the range 1.2 - 1.7 solar masses. We present the first stellar evolution calculations using a prescription for semiconvection derived from extrapolation of direct numerical simulations of double-diffusive mixing down to stellar parameters. The dominant mode of semiconvection in stars is layered semiconvection, where the layer height is an adjustable parameter analogous to the mixing length in convection. The rate of mixing across the semiconvective region is sensitively dependent on the layer height. We find that there is a critical layer height that separates weak semiconvective mixing (where evolution is well-approximated by using the Ledoux criterion) from strong semiconvective mixing (where evolution is well-approximated by using the Schwarzschild criterion). This critical layer height is much smaller than the minimum layer height expected from simulations so we predict that for realistic layer heights, the evolution is nearly the same as a model ran with the Schwarzschild criterion. We also investigate the effects of compositional gradient smoothing, finding that it causes convective cores to artificially shrink in the absence of additional mixing beyond the convective boundary. Layered semiconvection with realistic layer heights provides enough such mixing that stars will still evolve as if the Schwarzschild criterion is employed. Finally, we discuss the potential of detecting such semiconvection and its implication on convective core sizes in solar-like oscillators.

Fingering convection (or thermohaline convection) is a weak yet important kind of mixing that occurs in stably-stratified stellar radiation zones in the presence of an inverse mean-molecular-weight gradient. Brown et al. (2013) recently proposed a new model for mixing by fingering convection, which contains no free parameter, and was found to fit the results of direct numerical simulations in almost all cases. Notably, however, they found that mixing was substantially enhanced above their predicted values in the few cases where large-scale gravity waves, followed by thermo-compositional layering, grew spontaneously from the fingering convection. This effect is well-known in the oceanographic context, and is attributed to the excitation of the so-called "collective instability". In this work, we build on the results of Brown et al. (2013) and of Traxler et al. (2011b) to determine the conditions under which the collective instability may be expected. We find that it is only relevant in stellar regions which have a relatively large Prandtl number (the ratio of the kinematic viscosity to the thermal diffusivity), $O(10^{-3})$ or larger. This implies that the collective instability cannot occur in main sequence stars, where the Prandtl number is always much smaller than this (except in the outer layers of surface convection zones where fingering is irrelevant anyway). It could in principle be excited in regions of high electron degeneracy, during He core flash, or in the interiors of white dwarfs. We discuss the implications of our findings for these objects, both from a theoretical and from an observational point of view.

Regions of stellar and planetary interiors that are unstable according to the Schwarzschild criterion, but stable according to the Ledoux criterion, are subject to a form of oscillatory double-diffusive (ODD) convection often called "semi-convection". In this series of papers, we use an extensive suite of three-dimensional (3D) numerical simulations to quantify the transport of heat and composition by ODD convection, and ultimately propose a new 1D prescription that can be used in stellar and planetary structure and evolution models. The first paper in this series demonstrated that under certain conditions ODD convection spontaneously transitions from an initially homogeneously turbulent state into a staircase of convective layers, which results in a substantial increase in the transport of heat and composition. Here, we present simulations of ODD convection in this layered regime, we describe the dynamical behavior of the layers, and we derive empirical scaling laws for the transport through layered convection.

We have identified an important source of mixing in stellar radiation zones, that would arise whenever two conditions are satisfied: (1) the presence of an inverse vertical compositional gradient, and (2) the presence of density-compensating horizontal gradients of temperature (alternatively, entropy) and composition. The former can be caused naturally by any off-center burning process, by atomic diffusion, or by surface accretion. The latter could be caused by rotation, tides, meridional flows, etc. The linear instability and its nonlinear development have been well-studied in the oceanographic context. It is known to drive the formation of stacks of fingering layers separated by diffusive interfaces, called intrusions. Using 3D numerical simulations of the process in the astrophysically-relevant region of parameter space, we find similar results, and demonstrate that the material transport in the intrusive regime can be highly enhanced compared with pure diffusion, even in systems which would otherwise be stable to fingering (thermohaline) convection.

Iron-rich layers are known to form in the stellar subsurface through a combination of gravitational settling and radiative levitation. Their presence, nature and detailed structure can affect the excitation process of various stellar pulsation modes, and must therefore be modeled carefully in order to better interpret Kepler asteroseismic data. In this paper, we study the interplay between atomic diffusion and fingering convection in A-type stars, and its role in the establishment and evolution of iron accumulation layers. To do so, we use a combination of three-dimensional idealized numerical simulations of fingering convection, and one-dimensional realistic stellar models. Using the three-dimensional simulations, we first validate the mixing prescription for fingering convection recently proposed by Brown et al. (2013), and identify what system parameters (total mass of iron, iron diffusivity, thermal diffusivity, etc.) play a role in the overall evolution of the layer. We then implement the Brown et al. (2013) prescription in the Toulouse-Geneva Evolution code to study the evolution of the iron abundance profile beneath the stellar surface. We find, as first discussed by Th\'eado et al. (2009), that when the concurrent settling of helium is ignored, this accumulation rapidly causes an inversion in the mean molecular weight profile, which then drives fingering convection. The latter mixes iron with the surrounding material very efficiently, and the resulting iron layer is very weak. However, taking helium settling into account partially stabilizes the iron profile against fingering convection, and a large iron overabundance can accumulate. The opacity also increases significantly as a result, and in some cases ultimately triggers dynamical convection.

It has long been known that solar-type stars undergo significant spin-down, via magnetic braking, during their Main-Sequence lifetimes. However, magnetic braking only operates on the surface layers; it is not yet completely understood how angular momentum is transported within the star, and how rapidly the spin-down information is communicated to the deep interior. In this work, we use insight from recent progress in understanding internal solar dynamics to model the interior of other solar-type stars. We assume, following Gough and McIntyre (1998), that the bulk of the radiation zone of these stars is held in uniform rotation by the presence of an embedded large-scale primordial field, confined below a stably-stratified, magnetic-free tachocline by large-scale meridional flows downwelling from the convection zone. We derive simple equations to describe the response of this model interior to spin-down of the surface layers, that are identical to the two-zone model of MacGregor and Brenner (1991), with a coupling timescale proportional to the local Eddington-Sweet timescale across the tachocline. This timescale depends both on the rotation rate of the star and on the thickness of the tachocline, and can vary from a few hundred thousand years to a few Gyr, depending on stellar properties. Qualitative predictions of the model appear to be consistent with observations, although depend sensitively on the assumed functional dependence of the tachocline thickness on the stellar rotation rate.

Much progress has recently been made in understanding and quantifying vertical mixing induced by double-diffusive instabilities such as fingering convection (usually called thermohaline convection) and oscillatory double-diffusive convection (a process closely related to semiconvection). This was prompted in parts by advances in supercomputing, which allow us to run Direct Numerical Simulations of these processes at parameter values approaching those relevant in stellar interiors, and in parts by recent theoretical developments in oceanography where such instabilities also occur. In this paper I summarize these recent findings, and propose new mixing parametrizations for both processes that can easily be implemented in stellar evolution codes.

We perform coagulation & fragmentation simulations to understand grain growth in T Tauri & brown dwarf discs. We present a physically-motivated approach using a probability distribution function for the collision velocities and separating the deterministic & stochastic velocities. We find growth to larger sizes compared to other models. Furthermore, if brown dwarf discs are scaled-down versions of T Tauri discs (in terms of stellar & disc mass, and disc radius), growth at the same location with respect to the outer edge occurs to similar sizes in both discs.

A region of a star that is stable to convection according to the Ledoux criterion may nevertheless undergo additional mixing if the mean molecular weight increases with radius. This process is called fingering (thermohaline) convection and may account for some of the unexplained mixing in stars such as those that have been polluted by planetary infall and those burning helium-3. We propose a new model for mixing by fingering convection in the parameter regime relevant for stellar (and planetary) interiors. Our theory is based on physical principles and supported by three-dimensional direct numerical simulations. We also discuss the possibility of formation of thermocompositional staircases in fingering regions, and their role in enhancing mixing. Finally, we provide a simple algorithm to implement this theory in one-dimensional stellar codes, such as KEPLER and MESA.

We present the first numerical simulations of the solar interior to exhibit a tachocline consistent with the Gough and McIntyre (1998) model. We find nonlinear, axisymmetric, steady-state numerical solutions in which: (1) a large-scale primordial field is confined within the radiation zone by downwelling meridional flows that are gyroscopically pumped in the convection zone (2) the radiation zone is in almost-uniform rotation, with a rotation rate consistent with observations (3) the bulk of the tachocline is magnetic free, in thermal-wind balance and in thermal equilibrium and (4) the interaction between the field and the flows takes place within a very thin magnetic boundary layer, the tachopause, located at the bottom of the tachocline. We show that the thickness of the tachocline scales with the amplitude of the meridional flows exactly as predicted by Gough and McIntyre. We also determine the parameter conditions under which such solutions can be obtained, and provide a simple explanation for the failure of previous numerical attempts at reproducing the Gough and McIntyre model. Finally, we discuss the implications of our findings for future numerical models of the solar interior, and for future observations of the Sun and other stars.

Planet formation occurs within the gas and dust rich environments of protoplanetary disks. Observations of these objects show that the growth of primordial sub micron sized particles into larger aggregates occurs at the earliest stages of the disks. However, theoretical models of particle growth that use the Smoluchowski equation to describe collisional coagulation and fragmentation have so far failed to produce large particles while maintaining a significant populations of small grains. This has been generally attributed to the existence of two barriers impeding growth due to bouncing and fragmentation of colliding particles. In this paper, we demonstrate that the importance of these barriers has been artificially inflated through the use of simplified models that do not take into account the stochastic nature of the particle motions within the gas disk. We present a new approach in which the relative velocities between two particles is described by a probability distribution function that models both deterministic motion and stochastic motion. Taking both into account can give quite different results to what has been considered recently in other studies. We demonstrate the vital effect of two "ingredients" for particle growth: the proper implementation of a velocity distribution function that overcomes the bouncing barrier and, in combination with mass transfer in high-mass-ratio collisions, boosts the growth of larger particles beyond the fragmentation barrier. A robust result of our simulations is the emergence of two particle populations (small and large), potentially explaining simultaneously a number of long-standing problems in protoplanetary disks, including planetesimal formation close to the central star, the presence of mm to cm size particles far out in the disk, and the persistence of micron-size grains for millions of years.

The process referred to as "semi-convection" in astrophysics and "double-diffusive convection in the diffusive regime" in Earth and planetary sciences, occurs in stellar and planetary interiors in regions which are stable according to the Ledoux criterion but unstable according to the Schwarzschild criterion. In this series of papers, we analyze the results of an extensive suite of 3D numerical simulations of the process, and ultimately propose a new 1D prescription for heat and compositional transport in this regime which can be used in stellar or planetary structure and evolution models. In a preliminary study of the phenomenon, Rosenblum et al. (2011) showed that, after saturation of the primary instability, a system can evolve in one of two possible ways: the induced turbulence either remains homogeneous, with very weak transport properties, or transitions into a thermo-compositional staircase where the transport rate is much larger (albeit still smaller than in standard convection). In this paper, we show that this dichotomous behavior is a robust property of semi-convection across a wide region of parameter space. We propose a simple semi-analytical criterion to determine whether layer formation is expected or not, and at what rate it proceeds, as a function of the background stratification and of the diffusion parameters (viscosity, thermal diffusivity and compositional diffusivity) only. The theoretical criterion matches the outcome of our numerical simulations very adequately in the numerically accessible "planetary" parameter regime, and can easily be extrapolated to the stellar parameter regime. Subsequent papers will address more specifically the question of quantifying transport in the layered case and in the non-layered case.

To date, no self-consistent numerical simulation of the solar interior has succeeded in reproducing the observed thinness of the solar tachocline, and the persistence of uniform rotation beneath it. Although it is known that the uniform rotation can be explained by the presence of a global-scale confined magnetic field, numerical simulations have thus far failed to produce any solution where such a field remains confined against outward diffusion. We argue that the problem lies in the choice of parameters for which these numerical simulations have been performed. We construct a simple analytical magneto-hydrodynamic model of the solar interior and identify several distinct parameter regimes. For realistic solar parameter values, our results are in broad agreement with the tachocline model of Gough & McIntyre. In this regime, meridional flows driven at the base of the convection zone are of sufficient amplitude to hold back the interior magnetic field against diffusion. For the parameter values used in existing numerical simulations, on the other hand, we find that meridional flows are significantly weaker and, we argue, unable to confine the interior field. We propose a method for selecting parameter values in future numerical models.

Double-diffusive convection, often referred to as semi-convection in astrophysics, occurs in thermally and compositionally stratified systems which are stable according to the Ledoux-criterion but unstable according to the Schwarzchild criterion. This process has been given relatively little attention so far, and its properties remain poorly constrained. In this paper, we present and analyze a set of three-dimensional simulations of this phenomenon in a Cartesian domain under the Boussinesq approximation. We find that in some cases the double-diffusive convection saturates into a state of homogeneous turbulence, but with turbulent fluxes several orders of magnitude smaller than those expected from direct overturning convection. In other cases the system rapidly and spontaneously develops closely-packed thermo-compositional layers, which later successively merge until a single layer is left. We compare the output of our simulations with an existing theory of layer formation in the oceanographic context, and find very good agreement between the model and our results. The thermal and compositional mixing rates increase significantly during layer formation, and increase even further with each merger. We find that the heat flux through the staircase is a simple function of the layer height. We conclude by proposing a new approach to studying transport by double-diffusive convection in astrophysics.

We present the first three-dimensional simulations of fingering convection performed in a parameter regime close to the one relevant for astrophysics, and reveal the existence of simple asymptotic scaling laws for turbulent heat and compositional transport. These laws can straightforwardly be extrapolated to the true astrophysical regime. Our investigation also indicates that thermocompositional "staircases," a key consequence of fingering convection in the ocean, cannot form spontaneously in the fingering regime in stellar interiors. Our proposed empirically-determined transport laws thus provide simple prescriptions for mixing by fingering convection in a variety of astrophysical situations, and should, from here on, be used preferentially over older and less accurate parameterizations. They also establish that fingering convection does not provide sufficient extra mixing to explain observed chemical abundances in RGB stars.

When a planet falls onto the surface of its host star, the added high-metallicity material does not remain in the surface layers, as often assumed, but is diluted into the interior through fingering (thermohaline) convection. Until now, however, the timescale over which this process happens remained very poorly constrained. Using recently-measured turbulent mixing rates for fingering convection, I provide reliable numerical and semi-analytical estimates for the rate at which the added heavy elements drain into the interior. I find that the relative metallicity enhancement post-infall drops by a factor of ten over a timescale which depends only on the structure of the host star, and decreases very rapidly with increasing stellar mass (from about 1Gyr for a 1.3M_sun star to 10Myr for a 1.5M_sun star). This result offers an elegant explanation to the lack of observed trend between metallicity and convection zone mass in planet-bearing stars. More crucially, it strongly suggests that the statistically-higher metallicity of planet-bearing stars must be of primordial origin. Finally, the fingering region is found to extend deeply into the star, a result which would provide a simple theoretical explanation to the measurements of higher lithium depletion rates in planet-bearing stars.

The maintenance of large-scale differential rotation in stellar convective regions by rotationally influenced convective stresses also drives large-scale meridional flows by angular--momentum conservation. This process is an example of ``gyroscopic pumping'', and has recently been studied in detail in the solar context. An important question concerns the extent to which these gyroscopically pumped meridional flows penetrate into nearby stably stratified (radiative) regions, since they could potentially be an important source of non-local mixing. Here we present an extensive study of the gyroscopic pumping mechanism, using a combination of analytical calculations and numerical simulations both in Cartesian geometry and in spherical geometry. The various methods, when compared with one another, provide physical insight into the process itself, as well as increasingly sophisticated means of estimating the gyroscopic pumping rate. As an example of application, we investigate the effects of this large-scale mixing process on the surface abundances of the light elements Li and Be for stars in the mass range 1.3-1.5 solar masses (so-called ``Li-dip stars''). We find that gyroscopic pumping is a very efficient mechanism for circulating material between the surface and the deep interior, so much in fact that it over-estimates Li and Be depletion by orders of magnitude for stars on the hot side of the dip.However, when the diffusion of chemical species back into the surface convection zone is taken into account, a good fit with observed surface abundances of Li and Be as a function of stellar mass in the Hyades cluster can be found for reasonable choices of model parameters.

We propose a closure model for the transport of entropy and momentum in astrophysical turbulence, intended for application to rotating stellar convective regions. Our closure model is first presented in the Boussinesq formalism, and compared with laboratory and numerical experimental results on Rayleigh-Benard convection and Homogeneous Rayleigh-Benard convection. The predicted angular momentum transport properties of the turbulence in the slowly rotating case recover the well-known Lambda-effect, with an amplitude uniquely related to the convective heat flux. The model is then extended to the anelastic case as well as the fully compressible case. In the special case of spherical symmetry, the predicted radial heat flux is equivalent to that of mixing-length theory. For rotating stars, our model describes the coupled transport of heat and angular momentum, and provides a unified formalism in which to study both differential rotation and thermal inhomogeneities in stellar convection zones.

The solar convection zone exhibits a strong level of differential rotation, whereby the rotation period of the polar regions is about 25-30% longer than the equatorial regions. The Coriolis force associated with these zonal flows perpetually "pumps" the convection zone fluid, and maintains a quasi-steady circulation, poleward near the surface. What is the influence of this meridional circulation on the underlying radiative zone, and in particular, does it provide a significant source of mixing between the two regions? In Paper I, we began to study this question by assuming a fixed meridional flow pattern in the convection zone and calculating its penetration depth into the radiative zone. We found that the amount of mixing caused depends very sensitively on the assumed flow structure near the radiative--convective interface. We continue this study here by including a simple model for the convection zone "pump", and calculating in a self-consistent manner the meridional flows generated in the whole Sun. We find that the global circulation timescale depends in a crucial way on two factors: the overall stratification of the radiative zone as measured by the Rossby number times the square root of the Prandtl number, and, for weakly stratified systems, the presence or absence of stresses within the radiative zone capable of breaking the Taylor-Proudman constraint. We conclude by discussing the consequences of our findings for the solar interior and argue that a potentially important mechanism for mixing in Main Sequence stars has so far been neglected.

The rotation rate of the solar radiative zone is an important diagnostic for angular-momentum transport in the tachocline and below. In this paper we study the contribution of viscous and magnetic stresses to the global angular-momentum balance. By considering a simple linearized toy model, we discuss the effects of field geometry and applied boundary conditions on the predicted rotation profile and rotation rate of the radiative interior. We compare these analytical predictions with fully nonlinear simulations of the dynamics of the radiative interior, as well as with observations. We discuss the implications of these results as constraints on models of the solar interior.

We present a detailed numerical study of the Gough & McIntyre model for the solar tachocline. This model explains the uniformity of the rotation profile observed in the bulk of the radiative zone by the presence of a large-scale primordial magnetic field confined below the tachocline by flows originating from within the convection zone. We attribute the failure of previous numerical attempts at reproducing even qualitatively Gough & McIntyre's idea to the use of boundary conditions which inappropriately model the radiative--convective interface. We emphasize the key role of flows downwelling from the convection zone in confining the assumed internal field. We carefully select the range of parameters used in the simulations to guarantee a faithful representation of the hierarchy of expected lengthscales. We then present, for the first time, a fully nonlinear and self-consistent numerical solution of the Gough & McIntyre model which qualitatively satisfies the following set of observational constraints: (i) the quenching of the large-scale differential rotation below the tachocline - including in the polar regions - as seen by helioseismology (ii) the confinement of the large-scale meridional flows to the uppermost layers of the radiative zone as required by observed light element abundances and suggested by helioseismic sound-speed data.

We examine a physical process that leads to the efficient formation of gas giant planets around intermediate mass stars. In the gaseous protoplanetary disks surrounding rapidly-accreting intermediate-mass stars we show that the midplane temperature (heated primarily by turbulent dissipation) can reach > 1000 K out to 1 AU. Thermal ionization of this hot gas couples the disk to the magnetic field, allowing the magneto-rotational instability (MRI) to generate turbulence and transport angular momentum. Further from the central star the ionization fraction decreases, decoupling the disk from the magnetic field and reducing the efficiency of angular momentum transport. As the disk evolves towards a quasi-steady state, a local maximum in the surface density and in the midplane pressure both develop at the inner edge of the MRI-dead zone, trapping inwardly migrating solid bodies. Small particles accumulate and coagulate into planetesimals which grow rapidly until they reach isolation mass. In contrast to the situation around solar type stars, we show that the isolation mass for cores at this critical radius around the more massive stars is large enough to promote the accretion of significant amounts of gas prior to disk depletion. Through this process, we anticipate a prolific production of gas giants at ~1 AU around intermediate-mass stars.

Turbulent motions in the interior of a star play an important role in its evolution, since they transport chemical species, thermal energy and angular momentum. Our overall goal is to construct a practical turbulent closure model for convective transport that can be used in a multi-dimensional stellar evolution calculation including the effects of rotation, shear and magnetic fields. Here, we focus on the first step of this task: capturing the well-known transition from radiative heat transport to turbulent convection with and without rotation, as well as the asymptotic relationship between turbulent and radiative transport in the limit of large Rayleigh number. We extend the closure model developed by Ogilvie (2003) and Garaud and Ogilvie (2005) to include heat transport and compare it with experimental results of Rayleigh-Benard convection.

Meridional flows with velocities of a few meters per second are observed in the uppermost regions of the solar convection zone. The amplitude and pattern of the flows deeper in the solar interior, in particular near the top of the radiative region, are of crucial importance to a wide range of solar magnetohydrodynamical processes. In this paper, we provide a systematic study of the penetration of large-scale meridional flows from the convection zone into the radiative zone. In particular, we study the effects of the assumed boundary conditions applied at the convective-radiative interface on the deeper flows. Using simplified analytical models in conjunction with more complete numerical methods, we show that penetration of the convectively-driven meridional flows into the deeper interior is not necessarily limited to a shallow Ekman depth but can penetrate much deeper, depending on how the convective-radiative interface flows are modeled.

This series of papers investigates the early stages of planet formation by modeling the evolution of the gas and solid content of protostellar disks from the early T Tauri phase until complete dispersal of the gas. In this first paper, I present a new set of simplified equations modeling the growth and migration of various species of grains in a gaseous protostellar disk evolving as a result of the combined effects of viscous accretion and photo-evaporation from the central star. Using the assumption that the grain size distribution function always maintains a power-law structure approximating the average outcome of the exact coagulation/shattering equation, the model focuses on the calculation of the growth rate of the largest grains only. The coupled evolution equations for the maximum grain size, the surface density of the gas and the surface density of solids are then presented and solved self-consistently using a standard 1+1 dimensional formalism. I show that the global evolution of solids is controlled by a leaky reservoir of small grains at large radii, and propose an empirically derived evolution equation for the total mass of solids, which can be used to estimate the total heavy element retention efficiency in the planet formation paradigm. Consistency with observation of the total mass of solids in the Minimum Solar Nebula augmented with the mass of the Oort cloud sets strong upper limit on the initial grain size distribution, as well as on the turbulent parameter $\alphat$. Detailed comparisons with SED observations are presented in a following paper.

We construct a new set of self-consistent analytical disk models by taking into account both viscous and radiative sources of thermal energy. We analyze the non-isothermal structure of the disk across the mid-plane for optically thick disks, and use the standard two-temperature model in the case of optically thin disks. We deduce a set of general formula for the relationship between the mass accretion rate and the surface density profile. Our results recover those of Chiang & Goldreich in the optically thin regions, but extend their work for the opaque regions of the disk. For the purpose of illustration, we apply our theory in this paper to determine the structure of protostellar disks around T Tauri stars under a state of steady accretion and derive the corresponding radial distribution function of various disk properties such as surface density and temperature near the mid-plane. We calculate the position of the snow line around a sun-like T Tauri star, and deduce that it can evolve from well outside 10 AU during FU Orionis outbursts, to about 4 AU during passive accretion phase, to the present-day orbital radius of Venus and finally re-expand to over 2.2 AU during the protostellar- to-debris disk transition. This non-monotonous evolution of the snow line may provide some novel and deterministic explanation for the total water content and its isotopic composition of both Venus and the Earth. In the optically thin, outermost regions of the disk we find that the surface density profile of the dust varies roughly as 1/r, which is consistent with mm observations of spatially resolved disk of Mundy et al. (2000).

We investigate the nonlinear dynamics of turbulent shear flows, with and without rotation, in the context of a simple but physically motivated closure of the equation governing the evolution of the Reynolds stress tensor. We show that the model naturally accounts for some familiar phenomena in parallel shear flows, such as the subcritical transition to turbulence at a finite Reynolds number and the occurrence of a universal velocity profile close to a wall at large Reynolds number. For rotating shear flows we find that, depending on the Rayleigh discriminant of the system, the model predicts either linear instability or nonlinear instability or complete stability as the Reynolds number is increased to large values. We investigate the properties of Couette-Taylor flows for varying inner and outer cylinder rotation rates and identify the region of linear instability (similar to Taylor's), as well as regions of finite-amplitude instability qualitatively compatible with recent experiments. We also discuss quantitative predictions of the model in comparison with a range of experimental torque measurements. Finally, we consider the relevance of this work to the question of the hydrodynamic stability of astrophysical accretion discs.

We study the interaction between gas and dust particles in a protoplanetary disk, comparing analytical and numerical results. We first calculate analytically the trajectories of individual particles undergoing gas drag in the disk, in the asymptotic cases of very small particles (Epstein regime) and very large particles (Stokes regime). Using a Boltzmann averaging method, we then infer their collective behavior. We compare the results of this analytical formulation against numerical computations of a large number of particles. Using successive moments of the Boltzmann equation, we derive the equivalent fluid equations for the average motion of the particles; these are intrinsically different in the Epstein and Stokes regimes. We are also able to study analytically the temporal evolution of a collection of particles with a given initial size-distribution provided collisions are ignored.

We have developed a numerical code with which we study the effects of 2D perturbations on stellar structure. We present new numerical and analytical results on the heating of a main-sequence star in a binary system by its companion.

A quasi-steady state model of the consequences of rotation on the hydrodynamical structure of a stellar radiative zone is derived, by studying in particular the role of centrifugal and baroclinic driving of meridional motions in angular-momentum transport. This nonlinear problem is solved numerically assuming axisymmetry of the system, and within some limits, it is shown that there exist simple analytical solutions. The limit of slow rotation recovers Eddington-Sweet theory, whereas it is shown that in the limit of rapid rotation, the system settles into a geostrophic equilibrium. The behaviour of the system is found to be controlled by one parameter only, linked to the Prantl number, the stratification and the rotation rate of the star.

Gough & McIntyre have suggested that the dynamics of the solar tachocline are dominated by the advection-diffusion balance between the differential rotation, a large-scale primordial field and baroclinicly driven meridional motions. This paper presents the first part of a study of the tachocline, in which a model of the rotation profile below the convection zone is constructed along the lines suggested by Gough & McIntyre and solved numerically. In this first part, a reduced model of the tachocline is derived in which the effects of compressibility and energy transport on the system are neglected; the meridional motions are driven instead by Ekman-Hartmann pumping. It is shown that there exists only a narrow range of magnetic field strengths for which the system can achieve a nearly uniform rotation. The results are discussed with respect to observations and to the limitations of this initial approach. A following paper combines the effects of realistic baroclinic driving and stratification with a model that follows closely the lines of work of Gough & McIntyre.