A new stellar mixing process operating below shell convection zones following off-center ignition
DD raft version N ovember
15, 2018
Preprint typeset using L A TEX style emulateapj v. 11 / / A NEW STELLAR MIXING PROCESS OPERATING BELOW SHELL CONVECTION ZONES FOLLOWINGOFF-CENTER IGNITION
M. M oc ´ ak , C asey A. M eakin , E. M¨ uller & L. S iess Draft version November 15, 2018
ABSTRACTDuring most stages of stellar evolution the nuclear burning of lighter to heavier elements results in a ra-dial composition profile which is stabilizing against buoyant acceleration, with light material residing aboveheavier material. However, under some circumstances, such as o ff -center ignition, the composition profileresulting from nuclear burning can be destabilizing, and characterized by an outwardly increasing mean molec-ular weight. The potential for instabilities under these circumstances, and the consequences that they may haveon stellar structural evolution, remain largely unexplored. In this paper we study the development and evolutionof instabilities associated with unstable composition gradients in regions which are initially stable accordingto linear Schwarzschild and Ledoux criteria. In particular, we explore the mixing taking place under variousconditions with multi-dimensional hydrodynamic convection models based on stellar evolutionary calculationsof the core helium flash in a 1.25 M (cid:12) star, the core carbon flash in a 9.3 M (cid:12) star, and of oxygen shell burningin a star with a mass of 23 M (cid:12) . The results of our simulations reveal a mixing process associated with regionshaving outwardly increasing mean molecular weight that reside below convection zones. The mixing is not dueto overshooting from the convection zone, nor is it due directly to thermohaline mixing which operates on atimescale several orders of magnitude larger than the simulated flows. Instead, the mixing appears to be due tothe presence of a wave field induced in the stable layers residing beneath the convection zone which enhancesthe mixing rate by many orders of magnitude and allows a thermohaline type mixing process to operate on adynamical, rather than thermal, timescale. The mixing manifests itself in the form of overdense and cold blob-like structures originating from density fluctuations at the lower boundary of convective shell and ”shooting”down into the core. They are enriched with nuclearly processed material, hence leaving behind traces of highermean molecular weight. In these regions we find that initially smooth composition gradients steepen into stair-step like profiles in which homogeneous, mixed regions are separated by composition jumps. These step likeprofiles are then seen to evolve by a process of interface migration driven by turbulent entrainment. We discussour results in terms of related laboratory phenomena and associated theoretical developments. We also discussthe degree to which the simulated mixing rates depend on the numerical resolution, and what future steps canbe taken to capture the mixing rates accurately. Subject headings:
Stars: evolution – hydrodynamics – convection – mixing INTRODUCTION
The hydrodynamic approach to model certain phases ofstellar evolution, like e.g., a nuclear core flash, by numeri-cally solving the Euler or Navier-Stokes equations is essen-tially built upon first principles in physics, and thus (almost)parameter-free. This approach is advantageous compared to(1D) stellar evolutionary calculations when modeling phe-nomena related to buoyancy or shear driven flow instabilitieswhich are inherently multidimensional in nature. A shortcom-ing of this approach is that it remains impossible to simulate astar over evolutionary timescales. Instead, multi-dimensionalhydrodynamic calculations provide insight into the nature ofmixing processes operating over fluid dynamical timescaleswhich can be used to develop more comprehensive mixingmodels to be used in evolutionary codes. The classic mixinglength theory (MLT) for thermal convection, which remainsthe standard treatment of convection in stellar evolution cal- [email protected] Institut d’Astronomie et d’Astrophysique, Universit´e Libre de Brux-elles, CP 226, 1050 Brussels, Belgium Steward Observatory, University of Arizona,Tucson, AZ, 85721 USA Theoretical Division, Los Alamos National Laboratory, Los Alamos,NM 87545 USA Max-Planck-Institut f¨ur Astrophysik,Postfach 1312, 85741 Garching,Germany culations, is an example of one such mixing model based on,in this case, a very crude picture of turbulent flow. A grow-ing body of work presenting simulations of stellar interiors(e.g. Dearborn et al. 2006; Herwig et al. 2006, 2011; Meakin& Arnett 2007; Arnett et al. 2009; Moc´ak et al. 2008, 2009,2010) o ff ers the promise of moving beyond crude models likeMLT by providing theorists with detailed information aboutthe full non-linear development of fluid instabilities under re-alistic astrophysical conditions.We build on this work in the present paper by providing adetailed account of a new mixing process that we have dis-covered in simulations of shell convection. In particular, wepresent calculations of both helium and carbon shell convec-tion following o ff -center ignition. In both of these evolu-tionary phases we identify a coupling between the convec-tive shell and the layers which reside below the shell. Thiscoupling leads to mixing below the burning shell that has po-tentially significant consequences for the star’s structure andevolution.The paper is organized as follows. In Sect. 2 we describedour hydrodynamics code and input data and present some gen-eral features of the simulated flows. In Sect. 3 we provide adetailed account of the new mixing process that we have iden-tified and discuss its properties in the context of other pro-cesses commonly encountered in stellar interiors, including a r X i v : . [ a s t r o - ph . S R ] A ug Moc´ak et al.
TABLE 1S ome characteristic properties of our initial models : total stellar mass M , stellar population , metal content Z , mass of the mapped model M m , outer ( inner ) radius R m of the mapped model , and nuclear energyproduction rate L m of the mapped model , respectively .Model M Pop.
Z M m R m L m [M (cid:12) ] [M (cid:12) ] [10 cm] [10 L (cid:12) ]M 1 .
25 I 0 .
02 0 .
45 1 . . ∼
1L 9 . .
02 0 .
94 1 . . ∼ .
01O 23 . I ∼ .
02 2 .
67 1 . . ∼ penetration and overshoot, semiconvection, and salt finger-ing. Contact is made with relevant laboratory, geophysical,and fluid mechanics literature. We summarize our findings inSect. 4, and indicate future research directions. INITIAL DATA AND SIMULATION PROPERTIES
Initial data
We have simulated the reactive hydrodynamic evolutionfor three phases of stellar evolution. The initial data usedfor our calculations include: (1) the helium core of a 1.25M (cid:12) star during the peak (maximum core luminosity) of thecore helium flash computed with the GARSTEC code (Weiss& Schlattl 2000, 2007), (2) the carbon-oxygen (C-O) core ofa 9.3 M (cid:12) star at peak of the core carbon flash computed withthe STAREVOL code (Siess 2006), and (3) the core of a 23M (cid:12) star during oxygen shell burning (Meakin & Arnett 2007)computed with the TYCHO code (Arnett et al. 2010). Thethermodynamic structure of all three initial models is qualita-tively similar exhibiting an o ff -center temperature maximumdue to nuclear burning in a partially electron degenerate stel-lar core (see Fig. 1). Some additional properties of the stellarmodels are included in Table 1.The core helium flash occurs after central hydrogen exhaus-tion in the semi-degenerate helium core of low-mass stars(0 . (cid:12) ≤ M (cid:46) . (cid:12) ) due to an o ff -center ignition of he-lium by the triple- α . The core carbon flash ensues after centralhelium exhaustion and is also characterized by o ff -center car-bon ignition in a semi-degenerate carbon-oxygen (CO) coreof a rather massive star (7 M (cid:12) (cid:46) M ≤
11 M (cid:12) ), the domi-nant 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 (cid:38)
12 M (cid:12) ) close to corecollapse, is not ignited o ff -center, but follows an epoch ofcore oxygen burning which leaves behind a silicon-sulfur richsemi-degenerate core. Because of the steep temperature de-pendence of the nuclear energy production rates, convectiondevelops in the burning shell.During the core helium flash the convective shell is en-riched mainly in C which results in a negative mean molec-ular weight gradient below its base (Fig. 1 a) . The situa-tion is similar during the core carbon flash, where the nuclearash from carbon burning ( Ne and O mainly) increases themean molecular weight inside its convective shell relative tothat of the unburned inner CO core (Fig. 1 b). Oxygen shellburning on the other hand, which follows oxygen core burn-ing, results in a positive composition gradient, i.e., the meanmolecular weight increases in the direction of gravity at thebottom of the convective shell (Fig. 1 c).
Properties of hydrodynamic simulations The gradient grows for roughly 10 000 years since the onset of the corehelium flash. TABLE 2S ome properties of the and simulations : number of grid points in r ( N r ), θ ( N θ ), and φ ( N φ ) direction , radial grid resolution ∆ r , angular gridresolution ∆ θ and ∆ φ in theta and phi - direction , and maximumevolutionary time t max of the simulation , respectively .model N r N θ N φ ∆ r ∆ θ ∆ φ t max cm] [ ◦ ] [ ◦ ] [10 s]hefl.2d.1 180 90 - 5.55 1.33 - 30hefl.2d.2 270 180 - 3.70 1 - 30hefl.2d.3 360 240 - 2.77 0.75 - 120hefl.3d 180 90 90 5.55 1.33 1.33 6cafl.2d 360 180 - 1.95 0.5 - 60oxfl.2d 400 320 - 1.75 0.28 - 4.4 The numerical simulations were performed with a modi-fied version of the hydrodynamic code Herakles (Moc´ak et al.2008). The code employs the PPM reconstruction scheme(Colella & Woodward 1984), a Riemann solver for real gasesaccording to Colella & Glaz (1984), and the consistent multi-fluid advection scheme of Plewa & M¨uller (1999). Self-gravity, thermal transport and nuclear burning are included inthe code. Nuclear reaction networks are generated using theREACLIB library developed by Thielemann (private commu-nication). In Table 2 we summarize some characteristic pa-rameters of our 2D and 3D hydrodynamic simulations.The core helium flash simulations were performed on anequidistant spherical polar grid in 2D ( r , θ ) and 3D ( r , θ, φ ).While the angular grid of the axisymmetric 2D simula-tions hefl.2d.2 and hefl.2d.3 covered the full angular range, i.e., ◦ ≤ θ ≤ ◦ , the 2D simulation hefl.2d.1 was restrictedto the angular region 30 ◦ ≤ θ ≤ ◦ . The 3D simula-tion hefl.3d was performed within an angular wedge givenby 30 ◦ ≤ θ ≤ ◦ and − ◦ ≤ φ ≤ + ◦ , respectively.This allowed us to simulate the 3D model with a reasonableamount of computational time using a radial and angular res-olution comparable to that of the 2D model hefl.2d.1 for sev-eral convective turnover timescales. In all core helium flashsimulations the radial grid ranged from r = × cm to r = . × cm. Boundary conditions were reflective in alldirections except for simulations hefl.2d.1 and hefl.3d, wherewe utilized periodic boundary conditions in angular direc-tion(s) to avoid a numerical bias in case of wide angular con-vective structures. Abundance changes due to nuclear burningduring the He-flash are described by a reaction network con-sisting of the four α -nuclei He, C, O, and Ne coupledby seven reactions.The core carbon flash simulation cafl.2d was performed ona 2D equidistant spherical polar grid ( r , θ ) covering a 90 ◦ an-gular wedge centered at the equator ( θ = ◦ ), and a radialgrid ranging from r = × cm to r = × cm. Boundaryconditions were reflective in radial direction and periodic inangular direction. Abundance changes due to nuclear burn-ing were described by a reaction network consisting of H, He, C, N, O, Ne, Ne, Na, and Mg coupled by 17reactions.The oxygen shell burning simulation oxfl.2d was performedon a 2D equidistant spherical polar grid covering a 90 ◦ wedgecentered at the equator, and a radial grid ranging from r = × cm to r = × cm. Boundary conditions werereflective in radial direction and periodic in angular direction.Abundance changes due to nuclear burning were described bya reaction network consisting of neutrons, H, He, C, O, new mixing process operating below shell convection zones 3 F ig . 1.— Temperature T (solid) and mean molecular weight µ (dotted) as a function of radius for the initial core helium flash model (a), the core carbon flashmodel (b), and the oxygen shell burning model (c), respectively . In each of these panels the region of interest below the base of the convection zone is highlightedby the shaded vertical strip.F ig . 2.— Snapshots of the relative angular fluctuations of the mass fraction of He during the core helium flash for model hefl.2d.3 at t ∼ Cduring the core carbon flash for model cafl.2d at t ∼
682 s (b), and of O during oxygen shell burning for model oxfl.2d at t ∼
940 s (c), respectively. Thefluctuations are defined by ∆ X ( A N ) = × [ X ( A N ) − (cid:104) X ( A N ) (cid:105) θ ] / (cid:104) X ( A N ) (cid:105) θ , where A N ∈ { He , C , O } , and (cid:104)(cid:105) θ denotes the angular average at a givenradius. Ne, Na, Mg, Si, P, S, S, Cl, Ar, Ar, K, Ca, Ca, Ti, Ti, Cr, Cr, Fe, Fe, and Ni coupledby 76 reactions.The global character of the flow is illustrated in Figure 2which shows a snapshot of the composition fluctuations in thedeveloped convective flows from the hefl.2d.3 (helium burn-ing), cafld.2d (carbon burning), and oxfl.2d (oxygen burning)calculations. The basic topology of the flow is that of a tur-bulent convective shell surrounded by stably stratified layerswhich host internal waves. In these 2D simulations the shellconvection consists of interacting plumes and convective rollswhich span the depth of the shell. In 3D the flow is betterdescribed by plumes than rolls, since vortices are no longerpinned to meridional planes by the geometry of the calcula-tion (see e.g. Meakin & Arnett 2007). In both the 2D andthe 3D simulations, the internal waves can be seen as nar-row (in radius) bands extending in angle above and below theturbulent convective shell. The lack of a strong signature of internal waves above the convective shell in the helium andoxygen burning models in Figure 2 is a result of the compo-sition profile lacking a gradient there. The signature of wavesis present in other fields which have gradients, such as thepressure, density, and velocity fields.
Helium Flash Convection.—
The hydrodynamic propertiesof shell convection during the core helium flash are describedin detail in Moc´ak et al. (2008, 2009). Convection startsearly ( t < | v | ∼ v MLT (cid:46) × cm s − ), whereas the velocities in our2D models exceed those by up to a factor of four. Turbulententrainment leads to a growth of the width of the convectionzone on a dynamic timescale in both the 2D and 3D models Moc´ak et al. F ig . 3.— Maps of the relative angular fluctuation in carbon mass fraction ∆ X ( C) ≡ × ( X ( C) − (cid:104) X ( C) (cid:105) θ ) / (cid:104) X ( C) (cid:105) θ at the bottom of the He-flashconvection zone, where (cid:104)(cid:105) θ denotes the horizontal average at a given radius. From left to right, the following models are shown: (left) hefl.2d.3 at t = t = t = T max , and the dashed line gives the location from where our mixing begins. due to an exchange between the potential energy contained inthe stratified layers at the boundaries of the convection zoneand the kinetic energy of the turbulent flow inside the convec-tion zone. Carbon Flash Convection. —
An analysis of the hydrody-namic properties of shell convection during the core carbonflash, based on our 2D model cafl.2d, shows that the con-vective flow is dominated by small circular vortices. The an-gular averaged amplitudes of velocities inside the convectionzone | v | ∼ × cm s − exceed those predicted by MLT v MLT ∼ . × cm s − by about an order of magnitude.From authors experience, corresponding 3D convection (notcalculated) would be characterized by smaller velocities by afactor from 2 to 5. The operation of turbulent entrainment atthe convective boundaries has also been found but not inves-tigated in detail in this paper. Oxygen Shell Burning. —
The hydrodynamic properties ofshell convection during oxygen shell burning are discussedin detail in Meakin & Arnett (2007). Our 2D model oxfl.2dreproduces approximatelly the results of these authors. Theangular averaged amplitudes of the convective velocities in-ferred from model oxfl.2d are | v | ∼ × cm s − , whichis 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; formore details see Meakin & Arnett (2007). RESULTS
General Features of the Mixing
In this paper we focus our analysis on a mixing processwhich operates within the stably stratified layers residing be-neath the shell convection zones in both the carbon- andhelium-flash simulations. This mixing process is not observedto operate below the shell convection zone in the oxygen burn-ing case. In this subsection we describe the salient features ofthe mixing for the two shell flash calculations, and contrastthese with the oxygen shell burning case.
Helium Flash Convection.—
In the helium-flash convec-tion simulations we find a mixing process taking place be-low the temperature maximum immediately after convectionstarts. The mixing manifests as the formation of blob-likestructures in the stable layer just below the temperature maxi-mum, which fall inward toward the stellar center. These struc-tures are present in both the 2D and 3D simulations and canbe seen in the snapshots presented in Figure 3 which shows the fluctuations in the helium burning product C for three ofour He-flash simulations.The layers just beneath the base of the convection zoneshow the presence of overdense sinking blob-like structureswhich 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 new compositional step essentiallycoincides with the layer from which the mixing launches.Whereas it almost overlaps with the bottom of convectionzone in early phases of the mixing (Fig. 3; panels b and c ),it stays clearly separated later (Fig. 3; panel a ). The blobs arealways cooler and denser than the ambient matter by roughly0 .
4% and 0 . ∆ ρ/ρ within the region where ourmixing takes place we can derive an analytic estimate for thevelocity 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 ∆ ρ/ρ ∼ × − , inferred from oursimulation (Fig. 4 a), we find v ∼ (cid:115) g Λ ∆ ρρ ∼ × cm s − (1)where g ∼ cm s − is the gravitational acceleration at thebase of the convection zone, and Λ ∼ cm the approxi-mate length of the fingers. The above estimate agrees wellwith the results of an analysis of our simulation, where thegas inside the fingers sinks at a bit smaller velocities rang-ing from ∼ × cm s − up to 3 . × cm s − . It confirmsthat the “buoyancy work” ( ∼ v ) is consistent with the accel-eration of gas seen in the simulation. While the initial back-ground state is dynamically stable (according to Ledoux andSchwarzschild), these flow properties suggests that the mix-ing takes place on a dynamic timescale. This apparent contra-diction is addressed in § t ∼ × s is depicted in Figure 5a, d. Carbon Flash Convection. —
A similar mixing process alsoappears in the carbon-flash model. A snapshot showing var-ious quantities in the lower half of the carbon burning shellconvection zone and the underlying stable layer is presented Note that in Figs. 3 and 4, which show the finger-like structures, themaps are displayed in spherical coordinates ( r and θ in cm and degree) usinga rectangular projection. This makes the finger-like structures appear longerthan they really are. new mixing process operating below shell convection zones 5 F ig . 4.— Snapshots showing the base of the convection zone in the 2D model hefl.2d.3 at t = ff erence betweenthe local and the horizontally averaged value at a given radius of the mean molecular weight µ , the density ρ and the temperature T , respectively. The panels(c,e,f) display the deviation of angular velocity with respect to radius r dv θ / dr , radial velocity v R , and the Ledoux criterion (a positive value implies that the flowis Ledoux or dynamically unstable; see Eq. 2), respectively. The vertical solid line marks the bottom of the convection zone (location of T max ), while the dashedline gives the location from where our mixing begins. Note that only a part of the computational domain is shown. in Figure 5b, e and can be seen to develop in a similar mannerto the He flash model.In this model, however, the initial structure of the stable lay-ers is slightly di ff erent from that in the helium-flash model.When the model was mapped onto our hydrodynamic gridand allowed to adjust to hydrostatic equilibiurm, some narrowbands within the stable layers appeared which were unstableto convection. This was due merely to the di ff erences betweenthe stellar evolution code and hydrodynamics grid as well asthe fact that the stable layer was only marginally stable. Thesenarrow unstable regions quickly mixed, and stabilized, pro-ducing a slight stair stepping and shallow chemical disconti-nuities in the stable layer rather than a smooth profile as in theinitial stellar model. Besides this initial transient, however,the evolution between this model and the helium flash modelproceed in the same manner. Oxygen Shell Burning. —
In contrast to the helium- andcarbon-flash convection cases, the stable layer residing be-neath the oxygen shell burning zone does not undergo notice-able mixing over similar timescales. This is an interestingobservation because the overall structure of the oxygen burn-ing shell is remarkably similar to the flash convection cases,except for the mean molecular weight gradient, which is pos-itive in the oxygen burning case. This fact strongly suggeststhat the molecular weight gradient is responsible, or at least,connected to the mixing seen in the shell convection mod-els. We also point out that the buoyancy frequency N underconvection shell is significantly higher than in the previoustwo cases indicating stronger dynamic stability (Figure 5c, f),which could simply suppress the mixing. Linear Stability of the Mixing Layer
As a starting point to analyze the mixing we consider thelinear stability of the regions in question, the standard ap-proach to determine the mixing properties in 1D stellar evo-lution codes. For adiabatic displacements, the condition forstability is often written in terms of the temperature and com-positions gradients, and reads ∆ ∇ < ( ϕ/δ ) ∇ µ (2)where the stellar structure gradients are ∇ = d ln T / d ln P , ∇ µ = d ln µ/ d ln P , the thermodynamic derivatives are ∇ ad = ( ∂ ln T /∂ ln P ) s ,µ , δ = − ( ∂ ln ρ/∂ ln T ) P ,µ , ϕ = ( ∂ ln ρ/∂ ln µ ) P , s , and ∆ ∇ = ∇ − ∇ ad . This stability condition isknown as the Ledoux criterion for convective stability (Kip-penhahn & Weigert 1990) and it measures the degree to whichthe background is ”superadiabatic”. In a homogeneous com-position medium the stability condition reads ∆ ∇ <
0, whichis known as the Schwarzschild criterion and is equivalent tothe statement that positive entropy gradients ( ds / dr >
0) arestabilizing , negative entropy gradients are destabilizing, andzero entropy gradients result in neutral buoyancy upon dis-placement in the linear regime when pressure equilibrium ismaintained (Landau & Lifshitz 1959). The composition termin Eq. 2 accounts for the buoyancy due to mean molecularweight gradients in the medium.The stability condition (eq.[2]) can be identified with theregion in the ”nabla-plane” below the line described by ∆ ∇ = ( ϕ/δ ) ∇ mu , as illustrated in Figure 6. Above this line the fluidis subject to dynamical instability (e.g., thermal convection,Rayleigh-Taylor instabilities), while below the line the fluid is Moc´ak et al. F ig . 5.— Snapshots showing the relative angular fluctuations of the mean molecular weight ∆ µ/µ = ( µ − (cid:104) µ (cid:105) θ ) / (cid:104) µ (cid:105) θ at the base of the convection zone in the2D model hefl.2d.3 and t = t = t = N , and the solid line is the angular averaged meanmolecular weight µ . The horizontal dotted line corresponds to N = (cid:104)(cid:105) θ denotes the angular average at a given radius.F ig . 6.— Temperature ( ∆ ∇ = ∇ − ∇ ad ) and composition gradient ( φ/δ ) ∇ µ plane. The region above the diagonal is dynamically unstable according tothe Ledoux criterion (Eq. 2). The lower right quadrant is dynamically sta-ble, while the ”semiconvection” and ”thermohaline” regions are dynamicallystable, but unstable on thermal timescales. in a dynamically stable regime. But what about non-adiabatice ff ects? The regions in this plane labeled ”semiconvection”and ”thermohaline” are in fact unstable on thermal timescalesdue to secular instabilities , and their linear growth proper-ties have been analyzed for a sampling of astrophysical con- Secular instabilities are distinguished by growth timescales that are sig-nificantly longer than the dynamical timescale. ditions (Ulrich 1972; Kippenhahn et al. 1980; Grossman &Taam 1996; Charbonnel & Zahn 2007; Siess 2009; Stancli ff eet al. 2009).Kippenhahn et al. (1980) first recognized that o ff -centerhelium and carbon ignition would indeed lead to thermoha-line mixing below the burning shells but the long (thermal)timescale estimated for the growth of this instability com-pared to the rapid evolution of the shell flash led to a dimin-ished interest in this context. In our shell flash simulations,the layers below the shell convection are indeed unstable ac-cording to the thermohaline condition ( ∇ < ∇ ad and ∇ µ < the mixing process operating in our simulations is in factdistinct from direct thermohaline mixing , since it operates ona significantly shorter timescale than the thermal timescale.Instead, we find that the regions which undergo mixing areforced by waves excited within the region by the adjacent tur-bulent convection zone . Therefore, what we have captured inour simulations is a wave induced, turbulent mixing processoperating in a stably stratified layer. Wave Driven Mixing
As described in § § F ig . 7.— A time sequence showing the evolution of relative di ff erence between the local and the horizontally averaged value of mean molecular weight i.e., ∆ µ/µ (a,b,c), horizontally averaged square of the Brunt-V¨ais¨al¨a frequency N (solid), entropy S (dash-dotted) and µ (dashed) in the stable layer residing beneath theconvective shell in the He-flash model hefl.2d.3 at tree times t = × s (a,d), t = s (b,e) and t = s. (c,f). The horizontal dotted line marks zerolevel of N . tion, mixing is instigated which redistributes the entropy andthe composition. This is illustrated in the time sequence pre-sented in Figure 7 for the helium shell flash model hefl.2d.3.There are several notable features in this figure. First,the convection zone remains separated from the underlyingstable layer by a spike in buoyancy frequency N (aroundr ∼ . × cm in panels a and b ) arising from the entropyjump between the convection zone and the underlying layer.This entropy jump is due to the energy deposition by the nu-clear burning in the shell. Note that the stabilizing layer ispresent despite the destabilizing µ − profile. In other words, theentropy profile is stable enough to compensate for the destabi-lizing µ − gradient e ff ect. The spike in buoyancy also preventsovershoot mixing into the underlying stable layers.Another notable feature is the reduction in stability just be-low the buoyancy frequency spike which separates the regionfrom the convection zone as time progresses. This reductionin stability (or decrease in N ) is associated with a flatteningof the entropy profile in the stable layer due to mixing in-duced in the region. The mixing that ensues is also associatedwith a modification of the composition profile, and by t ∼ seconds a new step in the µ -profile has already formed near r ∼ . × cm.Finally, it is worth pointing out that in addition to the mix-ing inside the stable layer there is an overall deepening of theconvection zone due to turbulent entrainment. By t ∼ . × seconds the convective boundary has migrated inward a dis-tance of ∆ r ∼ × cm, from r ∼ . × cm to r ∼ . × cm.In the stable layers, the mixing is strongest in regions wherethe gradients are largest, as expected for a di ff usive mixingprocess. Since compositional mixing at the molecular level is not included in these simulations, nor is heat conduction, thismixing results as a consequence of numerical di ff usion. Westress here that the numerical di ff usion can have a realisticcounterpart in stars, but caution should be paid to the quanti-tative rate of mixing presented in this paper. We will discussthe connection between the simulated results (and numericalmixing) and the situation in real stars below in terms of waveenhanced di ff usion as studied by, e.g., Press (1981). Layering and Interface Migration
A striking phenomena seen in both the He- and C-flash con-vection models is the formation of homogeneous layers in themixing region separated by steep gradient interfaces. Theseare apparent in both Figure 5 and Figure 7. Once these lay-ers have formed, a subsequent phase of evolution takes placewhich involves the migration of the interfaces separating theselayers. This is illustrated in the series of snapshots compiledin Figure 8, which shows the time evolution of the buoyancyfrequency N in the mixing region for the carbon flash model.These N profiles and their peaks are seen to evolve by a pro-cess of corresponding interface migration driven by turbulententrainment. This entrainment is powered by the turbulent ki-netic energy deposited in the region by trapped gravity waveswhich are excited by the overlying convection zone, ratherthan thermally driven turbulence as in the case of a convec-tion zone.The formation of steps and interface migration is a wellknown phenomena in a variety of fluid dynamic systems, in-cluding those studied in geophysical systems, in laboratories,and in simulation. A very close analog to the type of inter-face migration that we see in our simulations is the situationmodeled by Balmforth et al. (1998, see their Figure 6 and Moc´ak et al. F ig . 8.— Profiles of the buoyancy frequency N as a function of radial co-ordinate for di ff erent values of the time. The vertical axis is time and theprofiles are shown spaced at time intervals of ∼
700 seconds. The inset figureat the top is showing the scale for N corresponding to the last profile at time ∼
13, which show striking similarities to our Figure 8). In theirwork, the evolution of buoyancy due to turbulence in a sta-bly stratified layer is modeled in an attempt to understand theorigin of layer formation and the dynamics of interface mi-gration. This situation is a very close analog to the mixingprocess taking place in our simulations, except that the sourceof turbulence induced in the mixing region in our calculationsis due to wave motions excited in the region.Although the quantitative rate of mixing in these regionsremains to be identified through higher resolution calcula-tions and analysis, the formation of steps and their migrationthrough the stable layer provide important clues about howone might treat this mixing in 1D stellar evolution codes (e.g.Balmforth et al. 1998; Zilitinkevich et al. 2007). DISCUSSION AND CONCLUSIONS
We have identified a mixing process operating in our sim-ulations of o ff -center helium- and carbon-flash convection.The mixing results in the redistribution of composition andentropy within the stable layers that reside beneath the shellconvection zones. Associated with the mixing process is theformation of ”fingers” of high mean molecular weight ma-terial which traverse the mixing region. These ”fingers” areessentially tracer of material trasport taking place across thestable layer (e.g., Figure 3).We find that this mixing takes place against a background state which is initially dynamically stable to linear perturba-tions. And while the mixing region is indeed unstable to ther-mohaline mixing on a secular timescale, the mixing processthat we observe operates on a dynamical timescale , thus ex-cluding thermohaline mixing as the root cause. This is re-inforced by the fact that heat conduction was not includedin these calculations, primarily because its e ff ects are tooslow to operate over the time frames simulated. The thermaltimescale for thermohaline mixing (Kippenhahn et al. 1980)estimated for our core helium flash model is of the order of10 years. It is neither expected that the timescale will drasti-cally change (decrease) with improved stellar models nor dur-ing the flash as the thermal structures are in general very sim-ilar (Sweigart & Gross 1978).Our search for a counterpart to the mixing observed inour simulations among terrestrial experiments with qualita-tively similar conditions led us to phenomena like viscous ordensity fingering driven by a Rayleigh-Taylor instability (DeWit 2004, 2008) and buoyancy-driven instabilities resultingfrom di ff erences in mass di ff usion coe ffi cients of chemically-reacting species sitting on top of each other (Almarcha et al.2010). Properties of these phenomena do not fit the mixing inour simulations.As an alternative explanation we describe a scenario inwhich mixing is instigated by the finite amplitude jostlingmotions of the stable layer by the turbulence present in theadjacent convective shell. The perturbations induced in thestable layer by the adjacent convection are primarily in theform of waves (i.e., g-modes). The jostling of the stable layerresults in a di ff usive mixing process therein which leads toa redistribution of the entropy and composition on a dynam-ical timescales. This scenario is qualitatively similar to theenhanced mixing discussed by Press (1981) and Press & Ry-bicki (1981) who found that thermal and mass transport is notonly a function of temperature and composition gradients, butalso a function of the fluid shearing motions setup by the pres-ence of a passing wave field.Indeed, the layer from which the mixing begins coincideswith layers of strongest shear (Fig. 4 c). But the shear andthe finite amplitude perturbations alone do not drive mixingand have to be preceded by the formation of su ffi ciently steepcompositional gradients. In the case of our core helium andcarbon flash models, such gradients are inherent in our initialstellar models at the bottom of convection shells due to o ff -center ignition.The biggest source of uncertainty inherent in this work isthe quantitative rate of mixing. This uncertainty stems fromthe fact that the simulated mixing results from enhanced nu-merical di ff usion, rather than from a resolved velocity fieldand molecular di ff usion processes. Higher resolution simula-tions as well as analytic modeling should be pursued, perhapsin restricted models, in order to quantify the mixing rate dueto this process. This is necessary for including it into a stellarevolution code.One of the pressing goals for stellar evolution is the con-struction of algorithms for 1D codes that capture mixing pro-cesses such as that described in this paper, and which extendmixing due to hydrodynamic processes beyond the bound-aries of convective regions and into stably stratified layers.The need to treat turbulence and mixing in stable layers hasalready been recognized in the geophysical community (e.g.,Fernando 1991; Zilitinkevich et al. 2007), and this work canserve to inform a related e ff ort for stellar evolution. new mixing process operating below shell convection zones 9The simulations were performed at the computer center ofthe Max-Planck-Society in Garching (RZG). Miroslav Moc´akacknowledges financial support from the Communaut´e fran-caise de Belgique - Actions de Recherche Concert´ees, andfrom the Institut d’Astronomie et d’Astrophysique at the Uni- versit´e Libre de Bruxelles (ULB). The authors want to thankto Christophe Almarcha and Anne De Wit for enlightning dis-cussions on chemo-hydrodynamics. We express our gratitudeto Rob Izzard for valuable discussions, and for helpful com-ments on the manuscript. REFERENCESAlmarcha, C., Trevelyan, P. M. J., Riolfo, L. A., Zalts, A., Hasi, C. E.,D’Onofrio, A. D., & De Wit, A. 2010, The Journal of Physical ChemistryLetters, 1, 752Arnett, D., Meakin, C., & Young, P. A. 2009, ApJ, 690, 1715Arnett, D., Meakin, C., & Young, P. A. 2010, ApJ, 710, 1619Balmforth, N. J., Smith, S. G. L., & Young, W. R. 1998, J.Fluid Mech., 335,329Charbonnel, C., & Zahn, J. 2007, A&A, 467, L15Colella, P., & Glaz, H. H. 1984, J.Comp.Phys., 59, 264Colella, P., & Woodward, P. R. 1984, J.Comp.Phys., 54, 174Dearborn, D. S., P., Lattanzio, J. C., & Eggleton, P. P. 2006, ApJ, 639, 405De Wit, A. 2004, These d’Agregation de l’Enseignement Superieur(Habilitation), 1, 0—. 2008, Chimie Nouvelle, 99, 1Fernando, H. 1991, Ann.Rev.Fluid Mech., 23, 455Grossman, S. A., & Taam, R. E. 1996, MNRAS, 283, 1165Herwig, F., Freytag, B., Hueckstaedt, R. M., & Timmes, F. X. 2006, ApJ,642, 1057Herwig, F., Pignatari, M., Woodward, P. R., Porter, D. H., Rockefeller, G.,Fryer, C. L., Bennett, M., Hirschi, R. 2011, ApJ, 727, 88Kippenhahn, R., Ruschenplatt, G., & Thomas, H. 1980, A&A, 91, 175Kippenhahn, R., & Weigert, A. 1990, Stellar Structure and Evolution(Stellar Structure and Evolution, XVI, 468 pp. 192 figs.. Springer-VerlagBerlin Heidelberg New York. Also Astronomy and Astrophysics Library) Landau, L. D., & Lifshitz, E. M. 1959, Fluid mechanics, ed. Landau,L. D. & Lifshitz, E. M.Meakin, C. A., & Arnett, D. 2007, ApJ, 667, 448Moc´ak, M., M¨uller, E., Weiss, A., & Kifonidis, K. 2008, A&A, 490, 265—. 2009, A&A, 501, 659Moc´ak, M., Campbell, S. W., M¨uller, E., & Kifonidis, K. 2010, A&A, 520,114Plewa, T., & M¨uller, E. 1999, A&A, 342, 179Press, W. H. 1981, ApJ, 245, 286Press, W. H., & Rybicki, G. B. 1981, ApJ, 248, 751Siess, L. 2006, A&A, 448, 717—. 2009, A&A, 497, 463Stancli ffff