Chaos and turbulent nucleosynthesis prior to a supernova explosion

With the continuing increase of computing power, it has become possible to numerically simulate — in three dimensions (3D) and with sufficient resolution to allow turbulent flow — aspects of the last stages of evolution of a massive star before core collapse. All such simulations show dynamic behavior beyond that contained in 1D stellar evolutionary models. Some of the dynamic behavior is intrinsic to turbulent convection, and some may be due to interactions between active burning shells. Because the star arranges itself to maintain an approximate balance between cooling by neutrino emission and heating by thermonuclear burning, the stages of burning of carbon, neon, oxygen and silicon are made shorter and more vigorous than they would otherwise be. These stages are crucial for nucleosynthesis, they set the initial conditions for core collapse to form a neutron star or black hole, and they affect any ensuing supernova explosion. It now appears that current ideas about these stages are quantitatively, and perhaps qualitatively, incorrect due to a flaw in the physics used in the 1D stellar evolutionary models: the treatment of convection and mixing.

The detailed numerical modeling of stars begins with work on mechanical calculators by Schwarzschild¹ and Hoyle,² and flourished with the application of digital computers in the work of Henyey et al.,³ Iben,^4,5 Kippenhahn & Wiegert⁶ and others. These efforts contained the assumption that mixing times were short in comparison to nuclear burning times, so that convection zones were completely mixed (homogeneous in composition). A key problem was the choice of method used to join the nearly adiabatic convection of the deep interior with the radiatively dominated, nonadiabatic convection which occurred as the stellar surface was approached. Erika Böhm-Vitense^7,8 provided an algorithm based on the mixing length ideas of Prandtl, which became the standard in the stellar evolution community, and known as “mixing-length theory” (MLT). MLT has one parameter which is actually adjusted, and in this sense is a one parameter theory; this adjustment is done so that a solar model would have the observed radius of the present day sun. Adjustment to get the observed luminosity constrained the solar abundances, at least to the extent that MLT was a physically correct description of solar convection. Böhm-Vitense invented MLT prior to two major discoveries relevant to turbulent flow: the work on the turbulent cascade appeared in the western literature in Kolmogorov,⁹ and work on chaos in a convective roll in Lorenz.¹⁰ A great success of MLT was the description of the “ascent of the giant branch”, in which stars evolve to larger radii and luminosity, with vigorous surface convection regions, using essentially the same value of the mixing length parameter as found from solar calibration.

What about more advanced stages of stellar evolution? Rakavy, Shaviv and Zinamon¹¹ examined the evolution of C and O cores in neutrino cooled stages up to the onset of core collapse instability, assuming complete mixing inside convective zones. Arnett¹² considered He cores (which formed C and O cores by burning He, giving a connection to photon cooled stages), introduced time dependent convection and mixing as an advective process, and evolved the cores into collapse and neutrino trapping.¹³ For shell burning, Eggleton¹⁴ introduced the notion of “convective diffusion” in which the advective action of turbulent convection is approximated by a diffusion operator, and emphasized that this was merely a numerical convenience (see §I G). Weaver, Zimmerman, & Woosley¹⁵ extended this notion to late burning stages of massive stars. Langer, El Eid & Fricke¹⁶ considered semiconvection as well, and offered a modified diffusion coefficient. This became popular and is now regarded as the standard in the field,¹⁷ although we could find no convincing justification for approximating a turbulent process by an heuristic diffusion in the astronomical literature. We note that Daly & Harlow¹⁸ derive the conditions for such diffusion, which include the constraint that the convective scale be small relative to the integral scale; this is not true in the stellar case, for which they are comparable.

Flows in stellar plasma are highly turbulent. The parameter which indicates the onset of turbulence is the Reynolds number

where ℓ is a length characterizing the flow, u a velocity scale characterizing the flow, and ν the kinematic viscosity. The Reynolds number is the ratio of inertial forces to viscous forces. As a representative example we use the Spitzer-Braginskii estimate of kinetic viscosity of a classical plasma (Braginskii;¹⁹ Spitzer²⁰ p. 146),

The coulomb logarithm is often taken to be ln Λ ∼ 10. For stellar conditions, ν ∼ 70, to be compared with ν ∼ 1/7, 1/100, 1/1000 for air, water, and liquid He. The higher viscosity of stellar plasma is more than compensated for by the larger length scales appropriate to stars. In the laboratory, there is a transition from laminar to turbulent flow at about

$Re∼102$ to 2 × 10³, depending upon the geometry of the boundaries. We note ν is essentially the product of a mean-free-path and a microscopic velocity. Microscopic velocities are of the order of the sound speed, so that macroscopic velocities are usually smaller (but not by much, e.g., Mach numbers ∼0.001 to 0.1 are typical). Because the macroscopic length scales characteristic of stars are so much larger than microscopic mean-free-paths (∼10¹⁸, e.g.),

${\cal R}_e \gg 2 \times 10^3$

$Re≫2×103$⁠, the flow in stars generally is expected to be highly turbulent.^21,22 In special cases strong rotation and strong magnetic fields may inhibit turbulence.

In addition to viscosity (the transport of momentum), thermal diffusivity (heat transport) and molecular diffusivity (composition transport) need to be considered. This hierarchy could be summarized by:

where χ, ν and D are the thermal diffusivity, kinematic viscosity and molecular diffusivity, respectively. This order is essentialy determined by the relevant cross-sections for the separate processes; collisional cross sections for ions determine the kinematic viscosity (ions carry most of the momentum in the plasma) and the molecular diffusivity (ions carry the compositional identity). In contrast, radiative cross-sections for electrons are important for thermal diffusivity. In direct numerical simulations (DNS), the zone size is chosen to be smaller than the smallest relevant length scale (Pope,²³ Ch. 9). Given finite computer power, this limits the size of the object that can be simulated. Stars are simply too big to be done with DNS. In this sense stellar simulations are under-resolved. All simulations of stars are large eddy simulations (LES), with the consequent issue of whether the sub-grid physics is properly represented. For turbulence the sub-grid physics implicit in a set of shock-capturing hydrodynamic algorithms (one of which we use) is in fact appropriate (Boris,²⁴ p. 9); such simulations are termed “implicit LES” or ILES. This good fortune is due in part to the nature of the dissipation in a turbulent cascade (Kolmogorov⁹), for which

is the rate that large scale kinetic energy is fed into the cascade. Here v_rms is the average turbulent velocity and ℓ_d is the characteristic linear dimension of the turbulent region. There is no explicit viscosity; the nonlinear flow adjusts itself to “hide” the viscous term in the Navier-Stokes equation, and essentially replaces it by the effective dissipation, −ε. The effective viscosity for the turbulent case has little to do with the viscosity for laminar flow.

What about the other processes? Viallet et al.²⁵ show that the transport properties are dominated by large scale coherent motions which are resolved in the body of the convective region. At the boundaries the current resolutions used may be too coarse, so that the mixing rates seen in the simulations may be affected by the numerical method and the resolution. There may be other physical processes that could remain completely unresolved in 3D simulations (e.g. rotationally induced shear mixing and MHD), for which discovery of a good algorithm (and correct physics) remains a challenge.

Viscosity smooths momentum gradients, so that a finite space-time grid automatically loses detailed information about momentum at the grid scale, and this gives rise to a “numerical viscosity”. Computer turbulence simulations have limits on the values of Reynolds numbers that may be simulated; we may define a “numerical Reynolds number”:

$[Re]num=ru/ΔrΔu∝n4/3$⁠, where n = r/Δr is the number of zones along a linear dimension r and

$\Delta u \sim (\epsilon \Delta r)^{\tfrac{1} {3}}$

$Δu∼(εΔr)13$ is estimated from a Kolmogorov⁹ cascade. For 3D flow, the computational cost is proportional to n⁴ (spacetime is 4D), so that the difference in cost between n ∼ 100 and n ∼ 1000 is a factor of 10,000. For n ∼ 100 numerical viscosity is significant (think molasses) while for n ∼ 1000 the flow is clearly turbulent (think water). The difference in computational cost is roughly equivalent to the power of a laptop compared to a supercomputer (2 cores to 20,000 cores).

John von Neumann²⁶ envisaged using computers to help understand turbulence; this is an example of his idea that numerical experiments, when studied, would inspire mathematical insight, resulting in improved theoretical understanding. Although early 2D simulations (Bazán & Arnett^27–30) indicated that 1D stellar evolution models had problems, 2D hydrodynamics can be dangerously unphysical except in special cases in which there is a physical process to enforce the restricted symmetry. The first 3D simulations with the same microphysics as 1D stellar models, and having adequate resolution to exhibit turbulent flow, were those of Meakin & Arnett,³¹ which modeled a sector of an oxygen burning shell during its early stages. The simulation developed convection quickly, and showed turbulent pulses (about 8 were computed at that time). The convection was not steady-state except in a coarse-grained, average sense. The average rate of buoyancy work was approximately balanced by the rate of kinetic energy dissipation expected from the turbulent cascade (Kolmogorov⁹), but the instantaneous values of dissipation lagged the driving by a significant fraction of a transit time for the convective zone (Arnett, Meakin & Young³²). This dynamic behavior was strikingly similar (see Arnett & Meakin³¹) to the famous, and far simpler model of a convective roll due to Lorenz,¹⁰ which is a classic example of deterministic chaos. A second striking feature of these simulations was the appearance of entrainment at the boundaries of the convective region, which was quantitatively similar to experimental results (see Meakin & Arnett³¹) but less like “overshoot” prescriptions used in stellar modeling. Vigorous wave generation also occurred at the convective zone boundaries.

Figures 1 and 2 give a sense of the turbulent behavior. Both simulations were of oxygen burning in a massive star, prior to silicon burning and core collapse, in the last day of the star's life. The energy available in oxygen burning is still sufficient at this stage to blow the star apart, if released explosively. Figure 1 illustrates turbulent mixing in a stratified burning shell having a depth of four pressure scale heights. The yellow color represents the abundance of

$S32$⁠, which is a product of oxygen burning. It is produced directly by the thermonuclear fusion of

$\rm ^{16}O$

$O16$ in the burning shell, and by erosion of the already burned core (the orange spherical surface) by entrainment. The turbulent region is bounded by stably stratified fluid and exhibits strongly dynamic wave motions. The

$\rm ^{32}S$

$S32$ abundance traces both large scale structure of advective motion, and smaller scale structure resulting from the instability and break-up of large scale motions into the turbulent cascade.

Figure 2 shows an oxygen burning shell with less stratification (two pressure scale heights in depth), illustrating the burning region in more detail. Again, the convective region is enclosed by stably stratified layers above and below. The blue isocontours show

$O16$ mass fraction, with the color map chosen to indicate the difference in value from top to bottom of the convective zone (clear at the bottom and dark at the top). The large distortions in the uppermost opaque contour of

$\rm ^{16}O$

$O16$ are due to surface waves, and illustrate how the unstable motions affect the stably stratified boundary layer. Plumes of material are dredged into the convection zone in localized down drafts. The bottom of the convective region is also active. Energy generation is indicated by a red color map. The region of most active burning (yellow) is thin. Below this another stably stratified region which also is buffeted by turbulent flow. Burning occurs intermittently, in cells which drive plumes of ascending matter, and is steady only in an average sense.

The numerical simulations are a powerful tool with which to study stars, but with limitations: it is not feasible at present to calculate a significant fraction of a stellar lifetime, even for the relatively short-lived massive stars. That takes too many hydrodynamic times, requiring too many time steps. If we can construct a reasonably realistic initial state, we can simulate interesting and rapidly evolving stages of a star's life. We need a better understanding of the convective flow (which we can directly simulate for short times) to develop a theory of stellar convection for long evolutionary times, which goes beyond the currently used MLT with its heuristic diffusion.

The Reynolds averaged Navier Stokes (RANS) equations provide a different and useful representation of the physics of turbulence. Turbulence shows vigorous fluctuations, but a robust underlying regularity. The RANS equations spring from the Navier-Stokes equations, which imply conservation of mass, momentum and energy in a fluid (Landau & Lifshitz³³). In the RANS approach a variable v is represented in terms of its average value ⟨v⟩ and its fluctuating value v′, so v = ⟨v⟩ + v′ and ⟨v′⟩ = 0. The separation of fluctuating and average quantities helps explore the underlying regularities. The averaging operation may involve averaging over time and/or space. While the RANS equations represent the full features of turbulence, they are not closed: each set of moment equations (e.g., for a second order moment ⟨v′v′⟩) generates a higher order moment (e.g., a third order moment ⟨v′v′v′⟩), so that another equation is needed, and so on (Tennekes & Lumley³⁴). The reason behind this lack of closure is the introduction of fluctuations v′ in an insufficiently constrained manner (all such fluctuations are not dynamically self-consistent, so some scheme to eliminate the unphysical ones is required).

An extensive set of mean field equations and data from 3D simulations can be found in Mocák et al.³⁵ 

One theoretical approach to complex problems is to build a simpler, more easily soluble model to capture some of the features of the complete problem. The numerical simulations of turbulence motivate the “turbulent vortex model” (Arnett & Meakin³⁶), both by quantitative analysis and visual inspection. The basic idea of the vortex model is to separate the two distinct aspects of turbulence: the large scale flow as exhibited in the Lorenz convective roll, and a dissipation at the much smaller Kolmogorov scale. The rate at which energy is fed into the turbulent cascade from large to small scales is given by the “four-fifths law” of Kolmogorov (Eq. (4) above), which has been called one of the few exact and fundamental results of turbulence theory.³⁷ This acts in conjuction with buoyant driving at the large (integral) scale, as described by the Lorenz model of a convective roll. A dynamic system may be constructed by equating the time derivative of the mean turbulent kinetic energy to the difference of this driving and damping. This system describes a convective cell (analogous to a thunderstorm cell), not the whole convection region which is comprised of many such cells. Finding the correct ensemble average of an array of such cells is a challenge at present. The turbulent vortex model has the advantage that it may be directly compared to some previous suggestions for time dependent convection (e.g., Gough,³⁸ Arnett,¹² Stellingwerf³⁹), and because it defines an acceleration equation, perhaps it may be generalized to include additional physical effects (composition gradients, rotation, and possibly magnetic fields).

Another, more systematic approach to solving the turbulence problem is to approximate the complete system by enforcing closure on the Reynolds averaged Navier Stokes (RANS) equations at some level. To be soluble an approximation must be made to relate the highest order moments with lower order ones, closing the system. While simple in principle, the practice is difficult. The RANS equations are exact, in principle, prior to closure; the approximation used to close them is the issue. The simulations are well behaved when analysed from the RANS point of view (Viallet et al.²⁵), and promise to further our theoretical understanding, especially if there is an interplay between the more formal and precise RANS approach and the more physically transparent simple modeling.

A similar problem has been an issue in atmospheric science for a long time (e.g., Lumley & Panofsky⁴⁰), and some success has resulted from closure by modifying the equations of third moments (Lumley, Zeman & Siess⁴¹). This is in contrast with MLT used in astrophysics, which only involves modifying the equations of second moments. Canuto⁴² has decryed this lag by the astrophysical community, and presented a formalism for the stellar case which he terms “Reynolds stress models”for these long evolutionary times (RSM), which are based on the geophysical studies and have a third order closure. They have not yet been widely applied, possibly because they are more complicated and there is “no widely accepted stellar evidence that they are superior to MLT for building stellar models”. It now appears, observationally, numerically, and theoretically, that this is not the case (Meakin & Arnett;³¹ Smith & Arnett,⁴³ and below). It is plausible that asteroseismological studies (e.g., Aerts et al.,⁴⁴ Aerts⁴⁵), which now are showing inadequacies of MLT-built stellar models, may be aided by more sophisticated and independently testable treatments of convection. Although our discussion and that of Canuto might appear at first sight to be different because of their very different origins, they seem to point toward the same sort of improved convection model (compare Canuto⁴² and Viallet et al.²⁵). However, stars are plasma and geophysical flows are mostly neutral fluids; these differences may become more apparent when magnetic fields are considered.

The nature of this turbulent velocity field u has a direct impact on stellar nucleosynthesis. Nuclear reactions conserve baryon number, but not mass. The continuity equation

represents the conservation of a quantity ρ in a relativistically invariant way (a four-vector). We define ρ = Σ_iN_iA_im_amu to express baryon number in atomic mass units, where N_i is the volumetric density of species i, A_i is the number of baryons in each i nucleus, and m_amu is the atomic mass unit. This ρ is numerically close to the “mass density”, ρ_m = Σ_iN_im_i, where m_i is the mass of each nucleus of the i species, as used in Landau & Lifshitz³³ for example. The difference is of order of the nuclear binding energy over A_im_amuc², summed over all i species; this is small and may be neglected except for extreme conditions (e.g., core collapse). While ρ is exactly conserved and relativistically invariant, ρ_m is not. Consider the abundance Y_i = N_i/ρN₀ for nucleus i, where N₀ = 1/m_amu is Avogadro's number. Using the continuity equation, the co-moving time derivative in the rest frame may be written as

The “diffusion” term is small for massive stars and may be neglected with regard to its action on the mean fields (however, see Michaud et al.;⁴⁶ Thoul et al.⁴⁷), but is essential for mixing material to the molecular scale in the turbulent convection zones. The “reaction” term is actually comprised of [the sum of all reactions which create nucleus i] minus [the sum of all reactions which destroy nucleus i]; this is a reaction network (Arnett,⁴⁷ Ch. 4). In a given zone in a stellar model, the time derivative of the abundance is

If the rate of mixing is rapid relative to the rate of nuclear reactions, the gradient of Y_i is zero and the second term may be neglected. Using the continuity equation, Eq. (7) and (6) may be written as a conservation equation for Y_i, with source and sink terms for nuclear reactions creating and destroying species i on the right hand side, and the divergence of a composition flux replacing the spatial derivative term, where the composition flux is

Notice that this definition of composition flux does not involve a spatial derivative. An expression for the mean rate of change of composition at a given stellar radius (as opposed to the instantaneous and local expressions of eq. (4) and (5)) can be developed after splitting fields into means and fluctuations, and we find

where the mean Lagrangian derivative is

the turbulent flux of composition of isotope i is

and the tilde and double primes indicate the Favrian mean and fluctuation, respectively (see Viallet et al.²⁵).

An important issue now arises: while we may equate Eq. (7) to local reaction source terms, Eq. (10) requires reaction source terms which account for both mean and fluctuating components of Y_i. Stellar codes almost always use the mean value

$Ỹi$ and ignore contributions due to fluctuations. This is valid when the effects of the fluctuations average to zero, or when the fluctuations are small, as might be expected for fast mixing (turnover time much smaller than reaction time). The late stages of stellar evolution are replete with fast reactions which come into balance with other reactions tending to cancel their effect. Mixing can affect this balance. While convective turnover may be “fast”, there will be some reactions which are faster as well as some which are slower. The fast reactions are local while the slow reactions are a global feature of the convection zone. The net effect of this needs to be explored more deeply, perhaps on a case by case basis. Proton ingestion into He burning regions (see II A) and the Urca process (e.g., Arnett⁴⁸) are interesting examples in which this might be significant.

At present it is more common in stellar evolution codes to replace the turbulent flux term

$F̃i$ by a heuristic operator for the divergence of “convective diffusion” (e.g., Eggleton,¹⁴ Weaver, Zimmerman, & Woosley,¹⁵ Langer, El Eid, & Fricke,¹⁶ Langer¹⁷),

where

$D=13ℓurms$ with ℓ being the mixing length and u_rms being the rms value of the turbulent velocity as estimated from MLT; this smooths composition gradients by mixing. This is equivalent to replacing the radial component of the composition flux by

and setting the other components to zero. This procedure is dangerous from a mathematical standpoint because it increases the order of the spatial derivative by one, and thus introduces spurious additional solutions. It is also dangerous from a physical standpoint because it makes the composition currents large in composition gradients (formally infinite at a contact discontinuity in composition; this must be avoided by setting u_rms to zero appropriately). In contrast, the numerically simulated composition currents are large in the interior of the convection zone, and smaller at the boundaries, a qualitatively different behavior. Because this approach to mixing is heuristic, it gives no help with the issue of the effect of fluctuations in composition in a turbulent burning region, and in addition, could be misleading with regard to estimated rates of mixing and burning. This is especially worrisome for the late stages, e.g., prior to core collapse in massive stars, which are more poorly constrained by observational data.

Why replace MLT with a more complex algorithm? Because there are phenomena which MLT does not describe correctly, and that have observational consequences which may already have been detected. Also, a classical problem in stellar pulsation theory is the role of convection in driving pulsations, which is ignored, as lamented by Unno et al.,⁴⁹ but seems to be related to chaotic convection.³⁶ 

Böhm-Vitense⁸ considered a well-mixed plasma, with uniform composition, so that MLT has no inherent ability to deal with nonuniform composition. In practice stellar evolution codes use the Schwarzschild criterion (independent of composition) or the Ledoux criterion (composition dependent), and sometimes a combination of both; the simulations do not seem to validate this picture. The problem is simply illustrated. Convective flow at low mach number may be regarded as buoyancy driven, where the buoyant acceleration is gρ^′/ρ. For an ideal gas at constant pressure (i.e., low Mach number flow), ρ^′/ρ + T ^′/T + Y ^′/Y = 0, so that fluctuations in composition (which are ignored) come in at the same level as temperature fluctuations (which are considered). Further, in a turbulent medium, composition fluctuations are homogenized, so the convective boundary is a special place, having composition gradients on one side and homogeneity on the other. The simulations show that such regions have complicated dynamic behavior which may not yet be resolved (see however, Woodward et al.⁵⁰), but give interesting suggestions to compare with asteroseismological data (e.g., see the discussion of the “overshoot” parameter in Aerts⁴⁵). This physics is also relevant to proton ingestion and nucleosynthesis in the AGB shell flash (Herwig,⁵¹ Stancliffe et al.,⁵² Woodward et al.⁵⁰), and in the He core flash (Moćak et al.⁵³).

There is a considerable body of observational evidence that many massive stars have eruptive episodes prior to core collapse (Smith & Arnett⁴³). Most obvious are the supernovae of types In and IIn, which require mass ejection just prior to core collapse in order to provide the mass of material which is seen to interact with the supernova ejecta. It has been shown that MLT is a steady state approximation, and that if this restriction is relaxed the convection in late evolution of massive stars is strongly dynamic and chaotic (Smith & Arnett,⁴³ Arnett & Meakin,³⁶ Meakin & Arnett⁵⁶). This may lead to mass loss, either steady but vigorous (Quataert & Shiode,⁵⁴ Shiode & Quataert⁵⁵) or eruptive, as observations indicate.

Gravitational collapse of stellar cores to form neutron stars or black holes is an initial value problem, so that it is important to have realistic initial values, i.e., realistic progenitor star models. The first multidimensional simulation of a massive star for a significant length of time was done by Meakin⁵⁷ (earlier work by Bazàn & Arnett^27–30 was limited to shorter intervals of time by lack of computer power). The results were so startling that extensive checking was required prior to publication (Arnett & Meakin⁵⁸). In these 2D simulations, which included active shells of C, Ne, O and Si burning, the static initial models began to convect, the flow became steadily more vigorous, building radial pulsations, and finally pushing so much matter off the grid that the simulations were terminated. Murphy et al.⁵⁹ had correctly shown that the nuclear burning was relatively stable to the ε-instability. Instead, this instability was driven by the time-dependent turbulent convection (τ-instability, Arnett & Meakin³⁶). All this action occurred near, but prior to, core collapse, and modified the progenitor significantly from the original 1D stellar model. Because 2D simulations may be misleading for 3D problems, these simulations need to be done in 3D, with a larger domain and a full 4π steradian grid. Until such simulations have been done and carefully analyzed, the detailed predictions of 1D stellar evolution codes for core collapse progenitors should be viewed with serious caution.

Couch & Ott⁶⁰ have shown that qualitative changes may result in simulations of neutrino driven explosions; they carried out a controlled set of simulations of the core collapse of a 15 M_⊙ dud of Woosley & Heger.⁶¹ Simply mapping a velocity structure similar to the convective rolls of Meakin & Arnett³¹ and Arnett & Meakin,⁵⁸ and slightly increasing the neutrino heating (5 percent) changed a dud into a successful explosion. Without the changed velocity structure but still with the slightly enhanced neutrino heating, there was no explosion. Couch & Ott⁶⁰ conclude that “initial conditions matter” for the core collapse problem. It is interesting that the same turbulence physics which is important for progenitors is also useful in understanding the behavior of 3D collapse simulations themselves (Murphy & Meakin⁶²). We may ask to what extent are the theoretical difficulties with core collapse supernova models due to incorrect progenitor (“initial”) models?

Meteoritic Anomalies. The meteoritic abundances, especially isotopic ratios, of carbonaceous chondritic meteorites have been a key empirical test for theories of nucleosynthesis. For example, Amari et al.⁶³ use quantitative details of 1D stellar models for estimates of nucleosynthesis yields. Simulations show that rather than one abundance pattern there are several clustered around a mean in a turbulent convection zone (Meakin & Arnett³¹). Beyond this there will be systematic shifts due to the difference between 3D models and 1D models. There will be “extra mixing” because of inadequate treatment of boundaries in 1D. Advection can give large scale transport without microscopic mixing, especially during violently dynamic evolution, so “layering” can be violated. Some of the existing puzzles in interpretation of meteoritic anomalies in abundance may be due to over-simplified and unphysical mixing in 1D models.

Novae. Kelly et al.⁶⁴ have suggested that classical novae may be used as indicators of mixing, and identify key abundance ratios which may be most useful. This analysis is entirely based on a 1D framework of initial and hydrodynamic models. This appears to be a promising approach, but needs to be expanded to include effects discussed above (3D is not 1D); 3D effects in (1) the initial hydrostatic model, (2) the approach to instability and (3) the hydrodynamic response may affect the predictions of nucleosynthesis.

Rotation. A fundamental problem in stellar evolution is the effect of rotation (Maeder & Meynet⁶⁵); stars are observed to have rotation rates ranging from slow to near breakup. Stars are plasma so that rotation should not be discussed without dealing with magnetic fields (the dynamo problem), as the solar surface and pulsars gently remind us. It is thought that gamma-ray bursts (GRB’s) are the consequence of rotating core collapse. The best attempts to model the progenitors of such events (e.g., Groh et al.⁶⁶) use a diffusion model of convection, which seems unlikely to be correct. There is a dimensionless ratio, the Rossby number (ratio of convective velocity to rotational velocity), which separates qualitatively different types of flow. This has not yet been explored in the stellar case for high numerical Reynolds numbers and realistic boundary conditions; see Miesch⁶⁷ for a review of simulations at moderate Reynolds numbers, and Woodward, Herwig & Lin⁵⁰ for the current state of the art.

The Sun. The Sun may provide clues. Goldreich et al.⁶⁸ showed the importance of solar convection in driving “solar-like” oscillations, which are key to helioseismology, and asteroseismology of stars similar to the sun. Helioseismology provides precise measures of the convection zone and the tachocline. Recent work by Balbus⁶⁹ seems to suggest that a simple extension of the turbulent vortex model may be possible, using the “wind equation”, which would describe the differential rotation within the Sun. This seems to be a prudent first step to complete as we progress toward understanding GRB’s.

Moore's law (see, e.g., Gottbrath et al.⁷⁰) has brought us to the point where present day computer power allows ILES simulations with sufficient resolution to give turbulent flow in 3D. We find new phenomena not predicted by standard stellar evolution codes, which use MLT. Theoretical analysis of the simulations promises a deeper understanding of nucleosynthesis, and of the dynamics of the late stages of evolution of massive stars (prior to core collapse and possible supernova explosion), as well as of mixing and burning in stars which are asteroseismological targets.

This work was supported in part by NSF Award 1107445 at the University of Arizona. One of us (DA) wishes to thank the Aspen Center for Physics (ACP) and the Kavli Institute for Theoretical Physics (KITP) for their hospitality and support.