A NEW STELLAR MIXING PROCESS OPERATING BELOW SHELL CONVECTION ZONES FOLLOWING OFF-CENTER IGNITION

We have simulated the reactive hydrodynamic evolution for three phases of stellar evolution. The initial data used for our calculations include: (1) the helium core of a 1.25 M_☉ star during the peak (maximum core luminosity) of the core helium flash computed with the GARSTEC code (Weiss & Schlattl 2000, 2008), (2) the carbon–oxygen (C-O) core of a 9.3 M_☉ star at peak of the core carbon flash computed with the STAREVOL code (Siess 2006), and (3) the core of a 23 M_☉ star during oxygen shell burning (Meakin & Arnett 2007) computed with the TYCHO code (Arnett et al. 2010). The thermodynamic structure of all three initial models is qualitatively similar, exhibiting an off-center temperature maximum due to nuclear burning in a partially electron degenerate stellar core (see Figure 1). Some additional properties of the stellar models are included in Table 1.

Zoom In Zoom Out Reset image size

The core helium flash occurs after central hydrogen exhaustion in the semi-degenerate helium core of low-mass stars (0.7 M_☉ ⩽ M ≲ 2.2 M_☉) due to an off-center ignition of helium by the triple-α. The core carbon flash ensues after central helium exhaustion and is also characterized by off-center carbon ignition in a semi-degenerate C-O core of a rather massive star (7 M_☉ ≲ M ⩽ 11 M_☉), the dominant nuclear reactions being ¹²C (¹²C, α) ²⁰Ne and ¹² C(¹²C, p) ²³Na followed by ¹⁶O (α, γ) ²⁰Ne. Oxygen shell burning, which is typical for massive stars (M ≳ 12 M_☉) close to core collapse, is not ignited off-center, but follows an epoch of core oxygen burning which leaves behind a silicon–sulfur-rich semi-degenerate core. Because of the steep temperature dependence of the nuclear energy production rates, convection develops in the burning shell.

During the core helium flash the convective shell is enriched mainly in ¹²C which results in a negative mean molecular weight gradient below its base (Figure 1(a)).⁵ The situation is similar during the core carbon flash, where the nuclear ash from carbon burning (²⁰Ne and ¹⁶O mainly) increases the mean molecular weight inside its convective shell relative to that of the unburned inner CO core (Figure 1(b)). Oxygen shell burning on the other hand, which follows oxygen core burning, results in a positive composition gradient, i.e., the mean molecular weight increases in the direction of gravity at the bottom of the convective shell (Figure 1(c)).

The numerical simulations were performed with a modified version of the hydrodynamic code Herakles (Mocák et al. 2008). The code employs the piecewise parabolic method reconstruction scheme (Colella & Woodward 1984), a Riemann solver for real gases according to Colella & Glaz (1984), and the consistent multi-fluid advection scheme of Plewa & Müller (1999). Self-gravity, thermal transport, and nuclear burning are included in the code. Nuclear reaction networks are generated using the REACLIB library developed by F. Thielemann (2006, private communication). In Table 2, we summarize some characteristic parameters of our two-dimensional (2D) and three-dimensional (3D) hydrodynamic simulations.

The core helium-flash simulations were performed on an equidistant spherical polar grid in 2D (r, θ) and 3D (r, θ, ϕ). While the angular grid of the axisymmetric 2D simulations hefl.2d.2 and hefl.2d.3 covered the full angular range, i.e., 0° ⩽ θ ⩽ 180°, the 2D simulation hefl.2d.1 was restricted to the angular region 30° ⩽ θ ⩽ 150°. The 3D simulation hefl.3d was performed within an angular wedge given by 30° ⩽ θ ⩽ 150° and −60° ⩽ ϕ ⩽ +60°. This allowed us to simulate the 3D model with a reasonable amount of computational time using a radial and angular resolution comparable to that of the 2D model hefl.2d.1 for several convective turnover timescales. In all core helium-flash simulations the radial grid ranged from r = 2 × 10⁸ cm to r = 1.2 × 10⁹ cm. Boundary conditions were reflective in all directions except for simulations hefl.2d.1 and hefl.3d, where we utilized periodic boundary conditions in angular direction(s) to avoid a numerical bias in case of wide angular convective structures. Abundance changes due to nuclear burning during the He flash are described by a reaction network consisting of the four α-nuclei ⁴He, ¹²C, ¹⁶O, and ²⁰Ne coupled by seven reactions.

The core carbon-flash simulation cafl.2d was performed on a 2D equidistant spherical polar grid (r, θ) covering a 90° angular wedge centered at the equator (θ = 90°) and a radial grid ranging from r = 3 × 10⁸ cm to r = 1 × 10⁹ cm. Boundary conditions were reflective in the radial direction and periodic in the angular direction. Abundance changes due to nuclear burning were described by a reaction network consisting of ¹H, ⁴He, ¹²C, ¹⁴N, ¹⁶O,²⁰Ne, ²²Ne, ²³Na, and ²⁴Mg coupled by 17 reactions.

The oxygen shell burning simulation oxfl.2d was performed on a 2D equidistant spherical polar grid covering a 90° wedge centered at the equator and a radial grid ranging from r = 2 × 10⁸ cm to r = 1 × 10⁹ cm. Boundary conditions were reflective in radial direction and periodic in angular direction. Abundance changes due to nuclear burning were described by a reaction network consisting of neutrons, ¹H, ⁴He, ¹²C, ¹⁶O, ²⁰Ne, ²³Na, ²⁴Mg, ²⁸Si, ³¹P, ³²S, ³⁴S, ³⁵Cl, ³⁶Ar, ³⁸Ar, ³⁹K, ⁴⁰Ca, ⁴²Ca, ⁴⁴Ti, ⁴⁶Ti, ⁴⁸Cr, ⁵⁰Cr, ⁵²Fe, ⁵⁴Fe, and ⁵⁶Ni coupled by 76 reactions.

The global character of the flow is illustrated in Figure 2 which shows a snapshot of the composition fluctuations in the developed convective flows from the hefl.2d.3 (helium burning), cafld.2d (carbon burning), and oxfl.2d (oxygen burning) calculations. The basic topology of the flow is that of a turbulent convective shell surrounded by stably stratified layers which host internal waves. In these 2D simulations the shell convection consists of interacting plumes and convective rolls which span the depth of the shell. In 3D the flow is better described by plumes than rolls, since vortices are no longer pinned to meridional planes by the geometry of the calculation (see, e.g., Meakin & Arnett 2007). In both the 2D and the 3D simulations, the internal waves can be seen as narrow (in radius) bands extending in angle above and below the turbulent convective shell. The lack of a strong signature of internal waves above the convective shell in the helium and oxygen burning models in Figure 2 is a result of the composition profile lacking a gradient there. The signature of waves is present in other fields which have gradients, such as the pressure, density, and velocity fields.

Zoom In Zoom Out Reset image size

Helium-flash Convection. The hydrodynamic properties of shell convection during the core helium flash are described in detail in Mocák et al. (2008, 2009). Convection starts early (t < 1000 s) and quickly extends over the whole convectively unstable region. In axisymmetry (i.e., 2D), the convection is characterized by fast and large circular vortices, while in 3D the convective flow is less ordered showing slower and smaller turbulent features. The convective velocities in our 3D model match those predicted by MLT quite well (|v| ∼ v_MLT ≲ 1 × 10⁶ cm s⁻¹), whereas the velocities in our 2D models exceed those by up to a factor of four. Turbulent entrainment leads to a growth of the width of the convection zone on a dynamic timescale in both the 2D and 3D models due to an exchange between the potential energy contained in the stratified layers at the boundaries of the convection zone and the kinetic energy of the turbulent flow inside the convection zone.

Carbon-flash Convection. An analysis of the hydrodynamic properties of shell convection during the core carbon flash, based on our 2D model cafl.2d, shows that the convective flow is dominated by small circular vortices. The angular averaged amplitudes of velocities inside the convection zone |v| ∼ 4 × 10⁶ cm s⁻¹ exceed those predicted by MLT v_MLT ∼ 1.5 × 10⁵ cm s⁻¹ by about an order of magnitude. From the authors' experience, corresponding 3D convection (not calculated) would be characterized by smaller velocities by a factor from 2 to 5. The operation of turbulent entrainment at the convective boundaries has also been found but not investigated in detail in this paper.

Oxygen Shell Burning. The hydrodynamic properties of shell convection during oxygen shell burning are discussed in detail in Meakin & Arnett (2007). Our 2D model oxfl.2d reproduces approximately the results of these authors. The angular averaged amplitudes of the convective velocities inferred from model oxfl.2d are |v| ∼ 1 × 10⁷ cm s⁻¹, which is roughly equivalent to the velocity v_MLT predicted by MLT. Turbulent entrainment has been detected in model oxfld.2d, too, but we have not further analyzed this phenomenon; for more details see Meakin & Arnett (2007).

In this paper, we focus our analysis on a mixing process which operates within the stably stratified layers residing beneath the shell convection zones in both the carbon- and helium-flash simulations. This mixing process is not observed to operate below the shell convection zone in the oxygen burning case. In this subsection we describe the salient features of the mixing for the two shell flash calculations and contrast these with the oxygen shell burning case.

Helium-flash Convection. In the helium-flash convection simulations we find a mixing process taking place below the temperature maximum immediately after convection starts. The mixing manifests as the formation of blob-like structures in the stable layer just below the temperature maximum, which fall inward toward the stellar center. These structures are present in both the 2D and 3D simulations and can be seen in the snapshots presented in Figure 3 which shows the fluctuations in the helium burning product C¹² for three of our He-flash simulations.

Zoom In Zoom Out Reset image size

The layers just beneath the base of the convection zone show the presence of overdense sinking blob-like structures which are enriched with higher mean molecular weight μ material than the ambient matter into which they penetrate, thereby creating a compositional step as they redistribute the μ.⁶ The beginning of a new compositional step essentially coincides with the layer from which the mixing launches. Whereas it almost overlaps with the bottom of the convection zone in early phases of the mixing (Figure 3, panels (b) and (c)), it stays clearly separated later (Figure 3, panel (a)). The blobs are always cooler and denser than the ambient matter by roughly 0.4% and 0.05%, respectively (Figure 4). Using the value of relative density fluctuations Δρ/ρ within the region where our mixing takes place we can derive an analytic estimate for the velocity of the dense blobs causing the finger-like structures if we assume that the blobs arise from a dynamic instability (Weiss et al. 2004). With Δρ/ρ ∼ 5 × 10⁻⁴, inferred from our simulation (Figure 4(a)), we find

$$\begin{equation} v \sim \sqrt{g \Lambda \frac{\Delta \rho }{\rho }} \sim 7\times 10^5 \mbox{\ cm s${}^{-1}$}, \end{equation} \tag{ 1 }$$

where g ∼ 10⁸ cm s⁻² is the gravitational acceleration at the base of the convection zone and Λ ∼ 10⁷ cm is the approximate length of the fingers. The above estimate agrees well with the results of an analysis of our simulation, where the gas inside the fingers sinks at a bit smaller velocities ranging from ∼1 × 10⁵ cm s⁻¹ up to 3.5 × 10⁵ cm s⁻¹. It confirms that the "buoyancy work" (∼v²) is consistent with the acceleration of gas seen in the simulation. While the initial background state is dynamically stable (according to Ledoux and Schwarzschild), these flow properties suggest that the mixing takes place on a dynamic timescale. This apparent contradiction is addressed in Section 3.2. The development of the flow at t ∼ 9 × 10⁴ s is depicted in Figures 5(a) and (d).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Carbon-flash Convection. A similar mixing process also appears in the carbon-flash model. A snapshot showing various quantities in the lower half of the carbon burning shell convection zone and the underlying stable layer is presented in Figures 5(b) and (e) and can be seen to develop in a similar manner to the He-flash model.

In this model, however, the initial structure of the stable layers is slightly different from that in the helium-flash model. When the model was mapped onto our hydrodynamic grid and allowed to adjust to hydrostatic equilibrium, some narrow bands within the stable layers appeared which were unstable to convection. This was due merely to the differences between the stellar evolution code and hydrodynamics grid as well as the fact that the stable layer was only marginally stable. These narrow unstable regions quickly mixed and stabilized, producing slight stair stepping and shallow chemical discontinuities in the stable layer rather than a smooth profile as in the initial stellar model. Besides this initial transient, however, the evolution between this model and the helium-flash model proceeds in the same manner.

Oxygen Shell Burning. In contrast to the helium- and carbon-flash convection cases, the stable layer residing beneath the oxygen shell burning zone does not undergo noticeable mixing over similar timescales. This is an interesting observation because the overall structure of the oxygen burning shell is remarkably similar to the flash convection cases, except for the mean molecular weight gradient, which is positive in the oxygen burning case. This fact strongly suggests that the molecular weight gradient is responsible or at least connected to the mixing seen in the shell convection models. We also point out that the buoyancy frequency N² under the convection shell is significantly higher than in the previous two cases indicating stronger dynamic stability (Figures 5(c) and (f)), which could simply suppress the mixing.

As a starting point to analyze the mixing, we consider the linear stability of the regions in question, the standard approach to determine the mixing properties in 1D stellar evolution codes. For adiabatic displacements, the condition for stability is often written in terms of the temperature and composition gradients, and reads

$$\begin{equation} \Delta \nabla < (\varphi /\delta)\nabla _\mu, \end{equation} \tag{ 2 }$$

where the stellar structure gradients are ∇ = dln T/dln P, ∇_μ = dln μ/dln P, the thermodynamic derivatives are ∇_ad = (∂ln T/∂ln P)_{s, μ}, δ = −(∂ln ρ/∂ln T)_{P, μ}, φ = (∂ln ρ/∂ln μ)_{P, s}, and Δ∇ = ∇ − ∇_ad. This stability condition is known as the Ledoux criterion for convective stability (Kippenhahn & Weigert 1990) and it measures the degree to which the background is "superadiabatic." In a homogeneous composition medium the stability condition reads Δ∇ < 0, which is known as the Schwarzschild criterion and is equivalent to the statement that positive entropy gradients (ds/dr > 0) are stabilizing, negative entropy gradients are destabilizing, and zero entropy gradients result in neutral buoyancy upon displacement in the linear regime when pressure equilibrium is maintained (Landau & Lifshitz 1959). The composition term in Equation (2) accounts for the buoyancy due to mean molecular weight gradients in the medium.

The stability condition (Equation (2)) can be identified with the region in the "nabla-plane" below the line described by Δ∇ = (φ/δ)∇_mu, as illustrated in Figure 6. Above this line the fluid is subject to dynamical instability (e.g., thermal convection, Rayleigh–Taylor instabilities), while below the line the fluid is in a dynamically stable regime. But what about non-adiabatic effects? The regions in this plane labeled "semiconvection" and "thermohaline" are in fact unstable on thermal timescales due to secular instabilities,⁷ and their linear growth properties have been analyzed for a sampling of astrophysical conditions (Ulrich 1972; Kippenhahn et al. 1980; Grossman & Taam 1996; Charbonnel & Zahn 2007; Siess 2009; Stancliffe et al. 2009).

Zoom In Zoom Out Reset image size

Kippenhahn et al. (1980) first recognized that off-center helium and carbon ignition would indeed lead to thermohaline mixing below the burning shells but the long (thermal) timescale estimated for the growth of this instability compared to the rapid evolution of the shell flash led to a diminished interest in this context. In our shell flash simulations, the layers below the shell convection are indeed unstable according to the thermohaline condition (∇ < ∇_ad and ∇_μ < 0) but the mixing process operating in our simulations is in fact distinct from direct thermohaline mixing, since it operates on a significantly shorter timescale than the thermal timescale. Instead, we find that the regions which undergo mixing are forced by waves excited within the region by the adjacent turbulent convection zone. Therefore, what we have captured in our simulations is a wave-induced, turbulent mixing process operating in a stably stratified layer.

As described in Sections 3.1 and 3.2, the regions immediately below the convective shells are initially (linearly) dynamically stable. Furthermore, the absence of heat conduction in these simulations precludes these regions from undergoing thermally driven instabilities, such as semiconvection and thermohaline mixing. Therefore, the instabilities which develop in these regions are most naturally attributable to the jostling of the layers by turbulence in the adjacent convection zone. The nature of the jostling is predominantly in the form of internal wave motions (i.e., gravity waves).

As the stable layer is jostled by the neighboring convection, mixing is instigated, which redistributes the entropy and the composition. This is illustrated in the time sequence presented in Figure 7 for the helium shell flash model hefl.2d.3.

Zoom In Zoom Out Reset image size

There are several notable features in this figure. First, the convection zone remains separated from the underlying stable layer by a spike in buoyancy frequency N² (around r ∼ 4.7 × 10⁸ cm in panels (a) and (b)) arising from the entropy jump between the convection zone and the underlying layer. This entropy jump is due to the energy deposition by the nuclear burning in the shell. Note that the stabilizing layer is present despite the destabilizing μ-profile. In other words, the entropy profile is stable enough to compensate for the destabilizingμ-gradient effect. The spike in buoyancy also prevents overshoot mixing into the underlying stable layers.

Another notable feature is the reduction in stability just below the buoyancy frequency spike which separates the region from the convection zone as time progresses. This reduction in stability (or decrease in N²) is associated with a flattening of the entropy profile in the stable layer due to mixing induced in the region. The mixing that ensues is also associated with a modification of the composition profile, and by t ∼ 10⁴ s a new step in the μ-profile has already formed near r ∼ 4.5 × 10⁸ cm.

Finally, it is worth pointing out that in addition to the mixing inside the stable layer there is an overall deepening of the convection zone due to turbulent entrainment. By t ∼ 9.5 × 10⁴ s, the convective boundary has migrated inward a distance of Δr ∼ 2 × 10⁷cm, from r ∼ 4.8 × 10⁸ cm to r ∼ 4.6 × 10⁸ cm.

In the stable layers, the mixing is strongest in regions where the gradients are largest, as expected for a diffusive mixing process. Since neither compositional mixing at the molecular level nor heat conduction is included in these simulations, this mixing results as a consequence of numerical diffusion. We stress here that the numerical diffusion can have a realistic counterpart in stars, but caution should be paid to the quantitative rate of mixing presented in this paper. We will discuss the connection between the simulated results (and numerical mixing) and the situation in real stars below in terms of wave-enhanced diffusion as studied by, e.g., Press (1981).

A striking phenomenon seen in both the He- and C-flash convection models is the formation of homogeneous layers in the mixing region separated by steep gradient interfaces. These are apparent in both Figures 5 and 7. Once these layers have formed, a subsequent phase of evolution takes place which involves the migration of the interfaces separating these layers. This is illustrated in the series of snapshots compiled in Figure 8, which shows the time evolution of the buoyancy frequency N² in the mixing region for the carbon-flash model. These N² profiles and their peaks are seen to evolve by a process of corresponding interface migration driven by turbulent entrainment. This entrainment is powered by the turbulent kinetic energy deposited in the region by trapped gravity waves which are excited by the overlying convection zone, rather than thermally driven turbulence as in the case of a convection zone.

Zoom In Zoom Out Reset image size

The formation of steps and interface migration is a well-known phenomena in a variety of fluid dynamic systems, including those studied in geophysical systems, in laboratories, and in simulation. A very close analog to the type of interface migration that we see in our simulations is the situation modeled by Balmforth et al. (1998; see their Figures 6 and 13, which show striking similarities to our Figure 8). In their work, the evolution of buoyancy due to turbulence in a stably stratified layer is modeled in an attempt to understand the origin of layer formation and the dynamics of interface migration. This situation is a very close analog to the mixing process taking place in our simulations, except that the source of turbulence induced in the mixing region in our calculations is due to wave motions excited in the region.

Although the quantitative rate of mixing in these regions remains to be identified through higher resolution calculations and analysis, the formation of steps and their migration through the stable layer provide important clues about how one might treat this mixing in 1D stellar evolution codes (e.g., Balmforth et al. 1998; Zilitinkevich et al. 2007).