Double-diffusive mixing in stellar interiors in the presence of horizontal gradients

# Double-diffusive mixing in stellar interiors in the presence of horizontal gradients

## Abstract

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.

hydrodynamics – instabilities – stars : interiors – stars : evolution

## 1 Introduction

Double-diffusive instabilities have long been known to drive mixing in stars. Schwarzschild & Härm (1958) and Kato (1966) first considered the combined effects of an unstable entropy gradient and a stable compositional gradient (gradient, hereafter) to study what is commonly known as “semi-convection”. Meanwhile, Ulrich (1972) and Kippenhahn et al. (1980) investigated the opposite scenario of a stable entropy gradient with an unstable gradient, known as “thermohaline” or “fingering” convection. In both cases, the overall density gradient is stable to dynamical convection (so the system is Ledoux-stable) but the rapid diffusion of heat compared with the slower molecular diffusion of composition can destabilize the fluid, and trigger some level of turbulent mixing (see the reviews by Garaud, 2013; Radko, 2013, for instance).

For the next few decades, however, the lack of experimental evidence prevented astrophysicists from estimating how much mixing these instabilities could generate. Only recently has it become possible to follow their nonlinear development numerically, and quantify the induced turbulent fluxes. Much progress has been made in this respect thanks to three-dimensional (3D) direct numerical simulations (DNS), both in the case of semi-convection (Mirouh et al., 2012; Wood et al., 2013) and of fingering convection (Traxler et al., 2011; Brown et al., 2013). In this work, we focus on the latter.

Brown et al. (2013) proposed a parameter-free, experimentally calibrated model for mixing by fingering convection in non-rotating stellar interiors, that has been successfully tested against 3D DNS. The only limitation of this model is that it assumes a purely one-dimensional background, with both entropy- and -gradients aligned with gravity. This assumption is consistent with the assumptions of most 1D stellar models, but does not necessarily represent actual conditions in stellar interiors. Indeed, a number of situations may arise in which weak horizontal gradients are present. As we demonstrate here, these can have significant impacts on mixing by double-diffusive instabilities.

Horizontal gradients are expected whenever a star is rotating. For instance, rotation affects the convective efficiency differently at different latitudes, presumably causing pole-to-equator entropy gradients in adjacent radiative regions. It also generates large-scale meridional flows, which burrow into the nearby radiative region and drive further inhomogeneities (e.g. Spiegel & Zahn, 1992; Wood & Brummell, 2012). In radiation zones, the centrifugal force causes a misalignment between isotherms and isobars, which also break the spherical symmetry – an effect at the source of Eddington-Sweet flows. Even if the star is rotating too slowly for this misalignment to be significant today, its spin-down history will have generated radial differential rotation, which is always accompanied by large-scale meridional flows, and therefore by lateral entropy gradients (Zahn, 1992; Oglethorpe & Garaud, 2013). Finally, interactions with binary companions cause tides, which are also sources of asymmetry. In these examples, horizontal entropy gradients are presumably also accompanied by horizontal compositional gradients, which prompts the question of how fingering instabilities are modified in their presence.

This question was in fact already addressed and mostly solved (at least in the linear regime) by Holyer (1983), who was interested in its application to mixing in the Earth’s oceans – the original thermohaline case. In this letter, we apply in Section 2 the linear stability analysis of Holyer (1983) to determine what effect horizontal gradients may have on double-diffusive mixing in stars. In Section 3 we carry out 3D DNS experiments of the process to study the nonlinear evolution of the instability, and the resulting turbulent mixing.

## 2 Model setup and linear stability analysis

We build on the local model for fingering convection used by Traxler et al. (2011), and Brown et al. (2013), and include the effects of horizontal gradients as in the work of Holyer (1983). We model a small region within a star using a local Cartesian coordinate system, with gravity defining the vertical direction: . We use the Boussinesq approximation (Spiegel & Veronis, 1960), which is valid as long as the domain height is much smaller than the local scaleheight.

The background temperature and composition profiles are given by and , where the local means , , and the local gradients , , and are all constant. In addition, we assume that the local adiabatic temperature gradient, , is constant. The total fields and are then defined as and .

For simplicity and consistency with the oceanographic case, we define the potential temperature as the temperature in excess of a theoretical adiabatic stratification: . Defining

and as before, we then have . Potential temperature plays a role analogous to that of entropy in the Boussinesq approximation.

Density perturbations are related to and through the linearized equation of state:

 ~ρρ0=−α~θ+β~μ , (2)

where and are the coefficients of thermal expansion and compositional contraction. Following Holyer (1983), we assume that there is no net horizontal density gradient, which implies

 αθ0x=βμ0x . (3)

So-called “density-compensating” horizontal gradients are common in the ocean (e.g. Stommel & Fedorov, 1967). Presumably, non-compensating gradients would drive fast horizontal flows that rapidly equalize density along isobars, turning the system into a density-compensated one. We assume here that similar arguments apply to stellar interiors.

When non-dimensionalized in as Traxler et al. (2011) based on: (1) the expected finger width, so , where is the thermal diffusivity and is the viscosity, (2) the thermal diffusion timescale across a finger, so , and (3) with and , the governing equations are:

 1Pr(∂u∂t+u⋅∇u)=−∇p+(θ−μ)ez+∇2u , ∂θ∂t+u⋅∇θ+su+w=∇2θ , ∂μ∂t+u⋅∇μ+su+wR0=τ∇2μ , ∇⋅u=0 , (4)

where is the non-dimensional velocity field, and where for simplicity from here on, , and represent the non-dimensional pressure, (potential) temperature and compositional perturbations away from the background.

Four non-dimensional numbers appear: the Prandtl number , the diffusivity ratio (where is the compositional diffusivity), and finally

 R0=αθ0zβμ0z and s=θ0xθ0z . (5)

is the density ratio, which measures the strength of the stabilizing vertical potential temperature (entropy) gradient compared with the destabilizing -gradient. Meanwhile, is the slope of constant background surfaces. The slope of planes of constant (background) is then equal to .

In stellar interiors, and can vary from down to , and is typically slightly greater than . The slope could reasonably be assumed to be between and , depending on the mechanism that generates the horizontal gradients. Finally, can take any value between one (near standard convective regions) and infinity, provided the region supports an inverse -gradient. It is important to note that, as a result, the slope of the constant- planes could in theory either be much smaller or much larger than one. However, we do not believe that the case is physically realistic in stars, and therefore always assume in what follows that .

Baines & Gill (1969) showed that in the absence of horizontal gradients, a system is only unstable to fingering convection if . Cases where are unstable to the Ledoux criterion, and therefore subject to dynamical convection. Cases where are completely stable to fingering. Holyer (1983) showed that horizontal gradients always act to destabilize the system. In stellar interiors, where , this statement applies regardless of the value of the density ratio. To see this, we seek normal form solutions of (4) of the kind . Looking for the most unstable modes of this system can be done without loss of generality assuming that (Holyer, 1983). This yields a cubic equation for the growth rate , of the form with

 a2=|k|2(1+Pr+τ) , (6) a1=|k|4(τPr+τ+Pr)+Prl2|k|2(1−1R0) , a0=|k|6τPr+Pr[l2(τ−1R0)+kls(1−τ)] ,

where and .

Figure 1 illustrates the behavior of the growth rate of the fastest growing mode, that is, the solution of the aforementioned cubic maximized over all values of and , as a function of the reduced density ratio defined as

 r=R0−1τ−1−1 . (7)

The case corresponds to the standard fingering instability, and illustrates how the unstable range is limited to , which corresponds to . However, we also find that unstable modes exist for as soon as . Secondly, we see that the growth rate of the fastest-growing mode for reaches a finite asymptotic value that scales with . In fact, it can be shown using asymptotic analysis that for , for and small and of the same order, and for , the fastest-growing modes have the following properties: is real, and have opposite signs, and

 Missing or unrecognized delimiter for \right (8)

where the functions , and are of order unity (for of order unity), and are related to one another. The actual expressions for , and are not particularly illuminating for our present purposes, however.

Since , we see from Equation (8) that while remains constant as decreases, tends to 0. Given that the unstable modes can simply be viewed as fingers whose orientation is inclined at an angle with respect to the horizontal, we see that they become more and more inclined as . Their width, which is reasonably well-approximated by , increases as decreases.

Figure 2 explains the basic instability mechanism. Note that all angles are here greatly exaggerated, to illustrate the process more clearly. The background potential temperature and compositional fields are stratified both horizontally and vertically with , leading to high-/high- material in the top right of the figure, and low-/low- material in the bottom left. Since density is constant on horizontal surfaces, any parcel of fluid can be displaced to the left without loss or gain of potential energy (1). Once displaced, the parcel rapidly loses heat to the surrounding fluid, but retains its higher content (2), as long as . It becomes denser than the surroundings, and sinks (3). The entire process takes place in reality in a continuous fashion along inclined planes, as shown below in Section 3 (see Figure 3), and causes a net transport of high- material from the top right to the bottom left of the domain (4). Note that similar steps (not shown), starting from the bottom left and moving a parcel first to the right, lead to an upward and rightward moving parcel.

## 3 Numerical simulations

In order to study the nonlinear development of the system, and estimate the rate of turbulent mixing induced by the instability, we turn to 3D DNS of Equations (4). We use the same code as in Traxler et al. (2011b) and Brown et al. (2013), modified by adding the terms containing in the temperature and compositional equations. This code assumes that all perturbations are triply-periodic in space. We present the results of a single simulation, with parameters: , , and . This corresponds to a reduced density ratio , which would be stable to fingering convection if was 0.

At these parameter values, the fastest-growing mode has , and . We therefore need a computational domain of horizontal length that is at least equal to in order to contain it. In what follows, we choose . This domain is fairly thin in the direction to save on computational time, but nevertheless sufficiently thick to capture the 3D dynamics of the problem (see Garaud and Brummell, in preparation). The resolution used has effective grid-points.

In the oceanographic case, secondary instabilities form that tend to have even larger aspect ratios (smaller ) than the linear modes. Containing them in the computational domain becomes prohibitive. Simeonov & Stern (2007) proposed that one can still model them numerically by using an inclined computational domain, with the tilt angle chosen so that the secondary modes are perfectly horizontal in the new reference frame. Preliminary simulations suggest that these modes at low have the same angle as the primary modes. We thus tilt the domain by the angle (see above) in the numerical simulations presented below. For clarity, any quantity measured in the tilted domain from here on is referred to with a prime (e.g. , , ).

Figure 3 shows successive snapshots of our simulation. At early times the behavior of the instability is as expected from linear theory. We find that the modes are horizontal in the inclined frame, and therefore inclined at an angle from the true horizontal. Alternating inclined fingers, whose width is roughly 30, carry high- material downward, and low- material upward. Very rapidly, shear develops between the fingers and destabilizes them, locally mixing both entropy and composition. This weakens the stratification so that the local density ratio drops below , at which point regions that are properly fingering-unstable, separated by thin diffusive interfaces, appear. In the oceanographic context, these structures are called “lateral intrusions” (see the review by Radko, 2013, Chapter 7). The over-dense part of each intrusion propagates in this plot downward and to the left, while the under dense part propagates to the right and upward. Later, we see that the fingering region from one intrusion impinges on the one below, and the weakest of the two disappears. This results in a progressive coarsening of the system, which exhibits fewer and fewer yet gradually stronger intrusions with time.

To better understand the structure of the intrusions, we show in Figure 4 various profiles at time (see final snapshot in Figure 3). We clearly see extended regions where the total temperature and compositional fields increase with , which are associated with weak along-intrusion flows, and a local density ratio that remains smaller than . These are the fingering layers observed in Figure 3. Separating them are thin diffusive interfaces where both gradients reverse sign. These interfaces could in principle be semi-convectively unstable, in the sense that they support both an unstable potential temperature gradient, and a stable compositional gradient, and that their density ratio lies above the critical value (Baines & Gill, 1969). However, the strong interfacial shear, combined with the fact that the region is quite thin, presumably stabilizes the interface against semi-convection (at least, in this simulation).

The fingering regions cause substantial vertical transport, bringing high-/high- material downward, and low-/low- material upward. This net flux must be accompanied by an equivalent interfacial flux, since we do not observe any pile-up just above the interfaces. The interfacial flux is primarily caused by shear instabilities, as pure diffusion across the interface actually acts in the opposite direction (because of the reversal of the sign of both gradients). Our findings are generally consistent with those in the oceanographic context (Simeonov & Stern, 2007).

Since is presumably very small in real stars, we can get an approximation of the amount of vertical transport in the system (i.e. along ) simply by considering the turbulent fluxes perpendicular to the intrusions (i.e. along ). The resulting (dimensional) potential temperature and compositional fluxes, in the tilted domain, are then simply given by

 FT≃−κTθ0z(1−⟨w′θ′⟩)≡−DTθ0z and Fμ≃−κμμ0z(1−R0τ⟨w′μ′⟩)≡−Dμμ0z , (9)

where the angular brackets denote a box-average. This defines the effective diffusivities and . Figure 5 shows that the turbulent fluxes gradually increase as the intrusions coarsen. While remains close to one, is of the order of 30 at the end of the simulation, showing that compositional mixing is significantly enhanced by the intrusive instability, but that heat transport is not. Since we have not been able to follow the coarsening beyond the 3-intrusion phase, it remains to be seen whether the latter has stalled, or simply proceeds on a much longer timescale than we were able to integrate the simulation for.

## 4 Conclusions and prospects

We have shown that horizontal entropy and compositional gradients in a stably-stratified radiation zone can trigger new double-diffusive instabilities, that evolve to form stacks of fingering layers separated by sheared diffusive interfaces, called intrusions. Turbulent compositional mixing in this regime is significantly more efficient than pure diffusion, and could explain some of the long-standing discrepancies between theory and observations of stellar abundances, as in the case of low-mass RGB stars for instance (e.g. Charbonnel & Zahn, 2007).

Much work remains to be done, however. First and foremost, we need to extend our numerical analysis towards more realistic parameters (lower , and ), and carry out a comprehensive study of parameter space as in Brown et al. (2013). This is needed to determine how the cross-intrusion transport rates scale with governing parameters, and with the intrusion heights. This task is computationally expensive, since the already-large domain size and resolution need to be increased further when , or decrease. We also need to study the conditions under which the small-scale primary instability (see Section 2) transitions to the large-scale intrusive instability. While many models have been proposed to explain the formation of oceanic intrusions (see, e.g. Walsh & Ruddick, 2000; Simeonov & Stern, 2007), it remains unclear whether they apply in astrophysics. The unified mean-field theory developed by Traxler et al. (2011b) is a promising approach in this respect. We also need to determine whether the intrusion coarsening process observed in our simulations proceeds indefinitely, or eventually stalls for a given intrusion height. Finally, we must determine whether any of the effects that were invoked as the possible origin of horizontal gradients (such as rotation, or slow large-scale flows) could have a sufficiently large influence on the development of both primary and secondary intrusive instabilities to invalidate our results. Nevertheless, these findings open an interesting set of new possibilities for mixing in stellar interiors, that had not been explored so far.

### References

1. Baines, P., & Gill, A. 1969, J. Fluid Mech., 37
2. Brown, J. M., Garaud, P., & Stellmach, S. 2013, ApJ, 768, 34
3. Charbonnel, C., & Zahn, J. 2007, Astron. Astrophys., 467
4. Garaud, P. 2013, in EAS Publications Series, Vol. 63, EAS Publications Series, 285–295
5. Holyer, J. Y. 1983, Journal of Fluid Mechanics, 137, 347
6. Kato, S. 1966, PASJ, 18, 374
7. Kippenhahn, R., Ruschenplatt, G., & Thomas, H. 1980, A&A, 91, 175
8. Mirouh, G. M., Garaud, P., Stellmach, S., Traxler, A. L., & Wood, T. S. 2012, ApJ, 750, 61
9. Oglethorpe, R. L. F., & Garaud, P. 2013, ApJ, 778, 166
10. Radko, T. 2013, Double-diffusive convection (Cambridge University Press)
11. Schwarzschild, M., & Härm, R. 1958, ApJ, 128, 348
12. Simeonov, J., & Stern, M. E. 2007, Journal of Physical Oceanography, 37, 625
13. Spiegel, E. A., & Veronis, G. 1960, ApJ, 131, 442
14. Spiegel, E. A., & Zahn, J.-P. 1992, A&A, 265, 106
15. Stommel, H., & Fedorov, K. N. 1967, Tellus, 19, 306
16. Traxler, A., Garaud, P., & Stellmach, S. 2011, ApJ, 728, L29
17. Traxler, A., Stellmach, S., Garaud, P., Radko, T., & Brummell, N. 2011b, Journal of Fluid Mechanics, 677, 530
18. Ulrich, R. K. 1972, ApJ, 172, 165
19. Walsh, D., & Ruddick, B. 2000, Journal of Physical Oceanography, 30, 2231
20. Wood, T. S., & Brummell, N. H. 2012, ApJ, 755, 99
21. Wood, T. S., Garaud, P., & Stellmach, S. 2013, ApJ, 768, 157
22. Zahn, J.-P. 1992, A&A, 265, 115
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters