and spectropolarimetry

Proceedings IAU Symposium No. 307, 2014

G. Meynet, C. Georgy, J.H. Groh &amp; Ph. Stee, eds.

c © 2014 International Astronomical Union DOI: 00.0000/X000000000000000X

## 3D and Some Other Things Missing from the Theory of Massive Star Evolution

## W. David Arnett 1

1

Steward Observatory University of Arizona Tucson AZ 85721, USA

email: darnett@as.arizona.edu

Abstract. This is a sketch of a 321D approximation which is nonlocal, and thus has nonzero fluxes of KE (to be published in more detail elsewhere). We plan to add this as an option to MESA. Inclusion of KE fluxes seems to help resolve the solar abundance problem (Asplund et al. 2009). Smaller cores may ease the explosion problems with core collapse supernova simulations.

Keywords. convection, turbulence, Sun: abundances, stars: evolution, supernovae

## 1. Introduction

Stars are three dimensional (3D), turbulent plasma, and much more complex than the simplified one dimensional (1D) models we use for stellar evolution. Computer power is inadequate † at present for adequately resolved (i.e., turbulent) 3D simulations of whole stars. In his review talk, Meakin (2015) has illustrated this complexity and shown how it may be tamed by use of 3D simulations and Reynolds-Averaged Navier-Stokes (RANS) equations.

A minimalist, more approximate step may be easier to implement in stellar evolutionary codes, and instructive for the closure problem. Formally, the RANS equations are incomplete unless taken to infinite order; they must be closed by truncation at low order to be useful. This may be due to the nature of the Reynolds averaging, which allows all fluctuations rather than restriction only to dynamically consistent ones. Closure requires additional information to remove these new, extraneous solutions. As a complement to the RANS approach, approximations which focus on dynamics are emphasized here. In the cascade, turbulent kinetic energy and momentum are concentrated in the largest eddies. We introduce a simple dynamically self-consistent model which contains the largest eddies and the Kolmogorov cascade, and examine the consequences.

Here we will concentrate on the following issues:

- ( a ) the turbulent cascade (Kolmogorov 1941, 1962),
- ( b ) chaos (Lorenz 1963) and time dependence (Arnett &amp; Meakin 2011b),
- ( c ) boundary layers (Prandtl &amp; Tietjens 1934) and surfaces of separation (Landau &amp; Lifshitz 1959),
- ( d ) combined mixing and burning, and
- ( e ) sensitivity of core collapse to progenitor structure

The important and related issues of coherent treatment of pulsations, eruptions and explosions, rotation and magnetic fields, and variable mass and angular momentum must be deferred.

† But perhaps close; see Herwig et al. (2013).

W. David Arnett

Table 1. A Few Examples of Three-dimensional Simulations of Convection in Stars.

| Attribute                                                                                                                           | 3D Atmospheres                                                          | Solar Convection Zone                                | Stellar Interiors                                                          |
|-------------------------------------------------------------------------------------------------------------------------------------|-------------------------------------------------------------------------|------------------------------------------------------|----------------------------------------------------------------------------|
| representative photosphere composition gradient nuclear burning magnetic field rotation driving geometry boundary inside grid hydro | Stein & Nordlund 1 Y(yes) N N Y N top box in star top only compressible | J. Toomre 2 N(no) N N Y Y top CZ in box No anelastic | Meakin & Arnett 3 N(no) Y Y N N top or bottom box in star Yes compressible |

Notes:

- 1 See Stein &amp; Nordlund (1998), and Magic et al. (2013, 2014) which also refers to other recent work.
- 2 See Brun et al. (2011, 2004), and many papers in the IAUS271 proceesings (Brun 2011).
- 3 See Meakin &amp; Arnett (2007); Viallet et al. (2013), and Arnett et al. (2014) for an overview.

Erika B¨hm-Vitense developed mixing-length theory in the 1950's (Vitense 1953; B¨hmo o Vitense 1958), prior to the publication in the west of Andrey Kolmogorov's theory of the turbulent cascade (Kolmogorov 1941, 1962). MLT might have been different had she been aware of the original work Kolmogorov had done in 1941. Edward Lorenz showed that a simple convective roll had chaotic behavior (a strange attractor, Lorenz (1963)). Ludwig Prandtl developed the theory of boundary layers (Prandtl &amp; Tietjens 1934). All these ideas will be relevant to our discussion, which is based upon theory and 3D simulations.

## 1.1. 3D Simulations

Numerical approaches for fluid dynamics simulations include implicit large eddy simulations (ILES) and direct numerical simulation (DNS). DNS resolves the smallest relevant scales, and procedes to larger scales until limited by computer power available. At present this limit is far smaller than stellar scales, so DNS is important in stellar problems for points of principle (e.g., Wood et al. (2013)). ILES includes the largest scales in the system, and procedes to smaller scales until limited by computer power available. Implict is the assumption that sub-grid scale phenomena are correctly treated. This seems to be valid for 3D turbulence using state-of-the-art methods (Grinstein, F. F. et al. 2007).

In order to be useful in a 1D stellar evolution code, 3D information must be projected onto a 1D coordinate system, hence '321D' . Table 1 gives a brief comparison of fea-† tures of some 3D simulations which are relevant to designing 1D stellar algorithms. The simulations have different strengths and weaknesses, and tend to complement each other.

3D atmospheres. Simulations of stellar atmospheres in 3D were pioneered by Stein &amp; Nordlund (1998), see also Magic et al. (2013) and references therein. These simulations are one of the great successes of radiation hydrodynamics, removing the necessity for the micro-turbulence and macro-turbulent fudge factors previously used to calculate stellar spectra. It is tempting to use MLT to fit such 3D simulations as a 'stellar engineering' exercise to connect atmospheres to interiors, because stellar evolution codes are almost always formulated in the language of MLT. Magic et al. (2014) give a clear discussion of this process, and show that MLT must be modified in at least one respect to make the identification: a ram pressure term must be added to MLT. It appears that such fits necessarily ignore a quantity important for stellar interiors which comes from the same

† The acronym '321D' for 'projection of three dimensions to one dimension' is due to John Lattanzio.

term but in the energy equation. This is the flux of turbulent kinetic energy, which is incorrectly defined as zero in MLT.

These 3D stellar atmospheres generally do not contain the whole convective zone, but use a lower boundary condition which has been shown to have little effect on the predicted spectra. By the same token, this means that these simulations are not a sensitive probe of the convection at the bottom of the convection zone.

Convection, rotation and MHD. Juri Toomre and his students and collaborators have pioneered anelastic simulations of the solar convection zone, focusing on global MHD and differential rotation. For 2D see (Hurlburt et al. 1984), and 3D (Brun et al. 2004). Unlike the 'box in a star' grids used in 3D atmospheres and stellar interiors, such global solutions make more serious demands on computer resources because of their larger extent, so that, other things being equal, their zoning tends to be coarser. Recent results suggest that the numerical viscosity is sufficiently low so that the flows are becoming realistically turbulent (Brun et al. 2011). Merging of the insights from these simulations with stellar evolution with rotation is a present challenge.

Stellar Interiors. Another option for use of limited computational resourses is to examine deep interiors in which composition gradients and nuclear burning occur. By using a 'box in a star' approach it is possible to include convective zone boundaries within the computational grid. The treatment of radiation flow is simpler (radiative diffusion). Due to neutrino cooling the thermal time scales become shorter, and make thermal relaxation easier to deal with; also higher luminosity as found in deep layers of red giants reduces the thermal time scale. Deep zones (large stratification) and pressure dilatation (Viallet et al. 2013) have been examined. The 'box in a star' approach truncates the lowest order modes, so that a natural complement is the pioneering work of Paul Woodward (e.g., Herwig et al. (2013)) to put a 'whole star in a box'.

## 2. A 321D Algorithm

This procedes in several steps:

- ( a ) add turbulent cascade,
- ( b ) add dynamic (acceleration) equation for integral scale,
- ( c ) balance between driving, damping, and the role of turbulent KE flux,
- ( d ) make quantitative connections to MLT and the Lorenz model,
- ( e ) use the steady-state Lorenz model to approximate average behavior,
- ( ) f use acceleration equation to define boundary behavior, and
- ( g ) add composition effects.

Because an acceleration equation is used, it is straight-forward (in principle at least) to add inertial forces (centrifugal and coriolis), Lorenz forces and differential rotation.

## 2.1. The Turbulent Cascade

Arnett et al. (2014) estimate the Reynolds number to be Re ∼ 10 18 at the base of the solar convection zone. Numerical simulations and laboratory experiments become turbulent for Re ∼ 10 , so fluid flows in stars are strongly turbulent if, as we assume 3 for the moment, rotational and magnetic field effects may be neglected. This special, simplier case is thought to be widely but not universally appropriate to stellar interiors.

For homogeneous, isotropic, and steady-state turbulence, the Kolmogorov relation between the dissipation rate of turbulent kinetic energy, velocity, and length scale is,

$$\epsilon & \sim v ^ { 3 } / \ell.$$

This is a global constraint, and applies to each length scale λ in the turbulent cascade,

for all scales λ , or,

$$\Delta v _ { \lambda } & \sim ( \epsilon \lambda ) ^ { \frac { 1 } { 3 } }.$$

so that the velocity variation across a scale λ is ∆ v λ , and increases as λ 1 3 . A description of the cascade needs both large and small scales; Eq. 2.3 implies that the largest (integral) scales have most of the KE and momentum while the smallest have the fastest relaxation times. Simulations confirm this (Arnett et al. 2009).

## 2.2. Dynamics

In MLT the buoyant force is approximately integrated over a mixing length to obtain an average velocity u (e.g., Smith &amp; Arnett (2014)),

$$u ^ { 2 } = g \beta _ { T } \Delta \nabla \left ( \frac { \ell _ { M L T } ^ { 2 } } { 8 H _ { P } } \right ).$$

Working backward, this may be expressed as

$$d u / d t = g \beta _ { T } \Delta \nabla - u / \tau,$$

where glyph[lscript] d ≡ glyph[lscript] 2 MLT / 8 H P and τ = glyph[lscript] d / u | | , and for ∆ ∇ &gt; 0. Multiplying by u gives a kinetic energy equation,

$$d ( u ^ { 2 } / 2 ) / d t = u \cdot g \beta _ { T } \Delta \nabla - u ^ { 2 } / \tau,$$

for which the steady-state solution † is Eq. 2.4 (only positive ∆ ∇ are allowed).

glyph[negationslash]

Had it been available, B¨hm-Vitense might have identified the damping term with the o Kolmogorov value (Eq. 2.1). However, Kolmogorov found the damping length glyph[lscript] d to be the depth of the turbulent region, so that it is not a free parameter, unlike MLT. There is a further issue; glyph[epsilon1] is the average dissipation rate, not the instantaneous value ( u /glyph[lscript] 3 d ) which fluctuates over time; that is, u = v except on average over τ (see Fig. 4 in Meakin &amp; Arnett (2007)); we return to this below.

## 2.3. Kinetic Energy Flux

The flow is relative to the grid of the background stellar evolution code, so

$$d ( u ^ { 2 } / 2 ) / d t = \partial _ { t } ( \mathbf u \cdot \mathbf u ) / 2 + \nabla \cdot \mathbf F _ { \text{KE} },$$

where F KE = ρ u u ( · u ) / 2 is a flux of kinetic energy. The generation of the divergence of a kinetic energy flux in this way is robust for dynamic models; it occurs in the more precise RANS approach as well (Meakin &amp; Arnett 2007). It does not occur in MLT, which assumes symmetry in velocity between upflows and downflows, a condition that is violated at even weak levels of stratification (Meakin &amp; Arnett 2010; Viallet et al. 2013). Alternatively, MLT is equivalent to Eq. 2.5 with du/dt = 0, the local approximation.

Viallet et al. (2013); Moc´k et al. (2014) find a a global balance between buoyant driving (at the largest scales) and turbulent damping (at the smallest). Because of the large separation of these length scales, they are weakly coupled, giving rise to fluctuating behavior typical of turbulence. This balance (on average only) is a new condition beyond MLT, and allows the elimination of the free parameter in MLT, the mixing length, or α = glyph[lscript]/H P . The actual flow patterns that correspond to this global balance depend upon the details of the turbulent energy input; nuclear heating and photospheric cooling give different flows (Meakin &amp; Arnett 2010), i.e., fluxes of kinetic energy.

- † Care must be taken with the sign of the transit time τ and the deceleration for negative u .

$$\epsilon \sim ( \Delta v _ { \lambda } ) ^ { 3 } / \lambda,$$

We choose a vector equation for the integral scale velocity u ,

$$\partial \mathbf u / \partial t + ( \mathbf u \cdot \nabla ) \mathbf u = \mathbf g \beta _ { T } \Delta \nabla - \mathbf u / \tau.$$

The damping term is chosen to give the Kolmogorov expression for glyph[epsilon1] , the turbulent dissipation, if Eq. 2.8 is dotted by u and averaged over a turnover time τ .

## 2.4. Connections to MLT and Lorenz

In the local, steady-state, limiting case, the left-hand side of Eq. 2.8 vanishes, and an equation similar to Eq. 2.4 results, but with the mixing length replaced by the turbulent damping length (essentially the lesser of the depth of the convective zone or 4 pressure scale heights, Arnett &amp; Meakin (2011b)). With this change, the cubic equation of B¨hmo Vitense may be derived (Smith &amp; Arnett 2014), and we recover a form of MLT.

If it is assumed that the integral scale motion is that of a convective roll, Eq. 2.8 may be reduced to the form of the classic Lorenz equations , but with a nonlinear damping term provided by the Kolmogorov cascade (Arnett &amp; Meakin 2011b). Because of the time lag, the modified equations become even more unstable than the original ones.

For 321D the stellar evolution code must be supplied a smooth time-averaged value for the convective variables; we approximate that by the steady state. The weak coupling between large scale driving and dissipation at the small scale causes time dependent fluctuations of significant amplitude in luminosity and turbulent velocity; see Meakin &amp; Arnett (2007), Fig. 4, for fluctuations in KE in an oxygen burning shell. The term ∂ u /∂t in Eq. 2.8 is needed for chaotic fluctuations, wave generation and large scale dynamic behavior. A stellar evolution code must step over the shorter turnover time scales (weather) to solve for the evolutionary times (climate); this requires an average over active cells, so that there is a cellular structure in both space and time. The steady-state limit of the Lorenz equation gives a reasonable approximation to its average behavior, filtering out the fluctuations (Arnett &amp; Meakin 2011b); we apply the same approximation to Eq. 2.8 for slow stages of stellar evolution.

## 2.5. Boundary Layers

Fluid motion in a star may be separated into two fundamentally different flows: solenoidal flow (divergence free) and potential flow (curl free) (Landau &amp; Lifshitz 1959). Fig. 1 shows the striking separation in flow velocities at such boundary layers (see Prandtl &amp; Tietjens (1934). We do not claim that we have resolved these boundary layers yet , but their † structure and nature are important to the rate at which turbulent flow moves into or from non-turbulent regions-the entrainment rate.

Peter Eggleton took an early step in dealing with steep gradients in composition with the introduction of an admittedly ad-hoc diffusion operator (Eggleton 1973); this numerically advantageous procedure has been widely adopted for stellar evolution.

Amore physical picture results from consideration of the dynamics of the motion. At its most elemental level, the velocity vector must turn at boundaries to maintain solenoidal flow. Most of the momentum is at the largest scales, where dissipation is least and the flow almost adiabatic, so this velocity must be the integral scale velocity. The magnitude of the acceleration required is just the centrifugal value v /w 2 where w is the width of the turning region and v the relevant velocity. Using Eq. 2.8 in the steady state limit, and taking w ∼ ∆ r &lt;&lt; glyph[lscript] , the radial component of the acceleration equation becomes

$$u _ { r } \partial u _ { r } / \partial r \sim \Delta ( \frac { 1 } { 2 } u _ { r } ^ { 2 } ) / \Delta r \sim g \beta _ { T } \Delta \nabla.$$

† We have fewer than 10 zones across the lower boundary layer; numerical viscosity may affect the computed entrainment rate.

Sorting out the signs we see that the buoyancy must be negative for the radial kinetic energy to decrease, so that negative buoyancy decelerates the motion, and gives the required turning. This defines a layer at which the velocity goes to zero on average. The 'overshoot' region has a width w ; this material is mixed because it moves back into the convective region after it turns.

There are fluctuations, so that the layer undulates, generating waves in the neighboring stable region. In a coordinate system moving with the undulating layer, there is a boundary between the mixed (overshoot) region and the overlying layer which is stably stratified. The overshoot region has solenoidal flow which flow in the stably stratified region is potential (wave) flow. By Kelvin's theorem (conservation of circulation), entropy change (e.g., dissipation, radiative diffusion) is required to move matter from one region to the other. Most of the turbulent momentum is in the largest scale convective motions, the integral scale. In simulations we see that these couple well with gravity waves † , whose speed increases with wavelength (Landau &amp; Lifshitz (1959), § 12), so that the longer wavelengths carry most of the energy. The shorter wavelengths are more dissipative; they are generated by nonlinear interaction of the long wavelengths (wave breaking; think water waves). This is an example of a mechanism for changing entropy at the boundary, allowing the turbulent region to grow. If this example is representative, such entrainment rates are not universal but depend upon local conditions at the boundary as well as turbulent velocities. See also Viallet et al. (2013) discussion of radiative heating at bottom of a deep (strongly stratified) convection zone.

- † Sound waves couple well only for higher Mach number flow (Landau &amp; Lifshitz (1959), § 64).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0326v1-figure-5510e2d3-9fb2-49f8-b210-a4cfcd89c29c.jpg)

A snapshot of turbulent convection Viallet et al. 20/3

Figure 1. Boundary Layers enclosing Convection (after Viallet et al. (2013).

Grolsch (2015) approaches this issue from a quite different point of view, but comes to some similar conclusions.

## 2.6. Mixing and Burning

Change in composition due to nuclear burning is a fundamental feature of stellar evolution. In general stars do not have uniform composition. In stellar evolution, mixing is determined according to the criterion of Schwarzschild, or of Ledoux. The Schwarzschild criterion is defined as S = ∇-∇ e , and the Ledoux criterion is L = S -∇ Y , so ∆ ∇ = L which reverts to S if there is no composition gradient ( ∇ Y = 0).

Mixing motions (solenoidal flow) result from buoyancy, pressure perturbations, or differential rotation. We focus on the buoyant acceleration, -g ρ /ρ. ′ Composition gradients enter the buoyancy on an equal basis with entropy gradients; for the simplest case of an ideal gas, ρ /ρ ′ = T /T ′ + Y /Y ′ -P /P, ′ where Y = 1 /µ is the number of free particles per baryon, or inverse mean molecular weight. For a more general equation of state, there are mutiplicative factors of order unity on the r.h.s. Traditionally P /P ′ is taken to be zero even though this is incorrect for strongly stratified cconvecton (Viallet et al. 2013).

Brunt frequency, which is important for for asteroseismology, is N 2 = -L ( gβ T /H P ), indicating a fundamental connection with boundary conditions (Aerts, C. et al. 2010).

## 3. Solar Abundances

In the steady-state but non-local case, the ( u ·∇ ) u term gives a coupling between driving regions and damping regions in the form of a flux of turbulent kinetic energy (Meakin &amp; Arnett 2007, 2010); this is assumed to be zero by symmetry in MLT (a major flaw). Even moderate stratification breaks the symmetry, and gives a finite flux of turbulent kinetic energy.

Any model of the Sun which uses MLT will neglect the effect of turbulent kinetic energy fluxes. To estimate the sign and size of the necessary changes, we scale from the simulations in Viallet et al. (2013) for a convection zone which like the Sun is highly stratified. The luminosity due to kinetic energy is L KE ∼ -0 35 . L near the base of the convection zone, which is significant. This negative luminosity must be balanced by increased positive radiative luminosity, to maintain a solar luminosity, so L new /L glyph[circledot] = 1 35. This may . obtained by a reduction in opacity. The opacity is well known, and depends upon the composition, primarily the metal abundance. If we simply assume that the opacity scales with metalicity, κ new /κ old ∼ 3 / 4, and the actual metalicity should shift from the stellar evolution value (from the 'standard solar model', SSM) to z new ∼ 0 02(3 . / 4) ∼ 0 015, . which is the value determined from 3D atmospheres. Asplund et al. (2009) conclude that the discrepancy between abundances determined from 3D atmospheres, and from stellar evolution is unidentified. The argument just given suggests that there is a significant flaw in the physics of SSM, correction for which tends to resolve the discrepancy. Until stellar models including kinetic energy flux are constructed, it seems reasonable to take the 3D atmospheric abundances as our best estimate, and attribute the disagreement to use of a local convection theory.

## 4. Sensitivity to Progenitor Structure

Simulations of core collapse may be sensitive to convection algorithms used to construct precollapse models. Couch &amp; Ott (2013) get increased tendency toward explosion simply by adding a nonradial velocity component to the progenitor model (Arnett &amp; Meakin

2011a), as turbulence requires. This also suggests changes to the initial models used by Kochanek (2014).

Fig. 2 shows two TYCHO models at the end of oxygen burning. The density structure is sensitive to the assumed entrainment physics. The two curves represent stars of 15 M glyph[circledot] , with only the entrainment physics changed. Not only are the sizes of the cores changed, but the carbon/oxygen ratio is affected by entrainment as well.

Core-collapse simulations often fail because of the high rate of infall of mantle matter onto the newly formed, nuclear density core. This matter must be photo-dissociated by the explosion shock if an explosion is to occur. Lower rates of infall aid the explosion. The rate of infall is

$$\dot { M } _ { i n } = u _ { i n } 4 \pi r _ { i n } ^ { 2 } \rho _ { i n }.$$

The infall velocity is essentially the local sound speed in the progenitor mantle; matter falls in as a rarefaction wave moves out. To make the rate of mass infall ˙ M in small, the initial density ρ in should be small. The critical time occurs during the infall of the mantle (labeled 'Infall'); the envelope falls in too slowly to affect the explosion shock. Clearly the curve in Fig. 2 with the larger core will have larger ρ in , and be harder to explode. This large core resulted from taking the maximum entrainment rate that was energetically allowed, throughout the evolution from the main sequence (for simplicity no mass loss or binary stripping were considered).

The small core case resulted from an entrainment rate of 0.01 of the maximum, so that the boundary dynamics was almost elastic, and consistent with analytic estimates. Smaller rates and smaller cores are possible. The highest resolution 3D simulations we

Figure 2. Changes in Density structure of 15 M glyph[circledot] with Entrainment.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0326v1-figure-99b4a4bf-df0e-4d92-ae6f-61872240ef8e.jpg)

have done were the most nearly elastic in the boundary layers, so we regard this as the more reasonable case. It has smaller cores that conventional progenitor models; ad-hoc diffusion smoothes gradients, leading to larger cores.

## 5. Conclusions

This is a sketch of a 321D approximation which is nonlocal, and thus has nonzero fluxes of KE (to be published in more detail elsewhere). We plan to add this as an option to MESA. Inclusion of KE fluxes seems to help resolve the solar abundance problem (Asplund et al. 2009). Smaller cores may ease the explosion problems with core collapse supernova simulations.

## References

Aerts, C., Christensen-Dalsgaard, J., &amp; Kurtz, D. W. 2010, Asteroseismology , Springer

Arnett, D., Meakin, C., &amp; Young, P. A. 2009, ApJ 690, 1715

Arnett, W. D. &amp; Meakin, C. 2011a,

ApJ

733, 78

Arnett, W. D. &amp; Meakin, C. 2011b, ApJ 741, 33

Arnett, W. D., Meakin, C., &amp; Viallet, M. 2014, AIP Advances 4(4), 041010

Asplund, M., Grevesse, N., Sauval, A. J., &amp; Scott, P. 2009, ARAA 47, 481

B¨ ohm-Vitense, E. 1958, ZfA 46, 108

Brun, A. S. (ed.) 2011, Astrophysical Dynamics: From Stars to Galaxies , Vol. S 271, International Astronomical Union, CUP

Brun, A. S., Miesch, M. S., &amp; Toomre, J. 2004, ApJ 614, 1073

Brun, A. S., Miesch, M. S., &amp; Toomre, J. 2011, ApJ 742, 79

Couch, S. M. &amp; Ott, C. D. 2013, ApJ 778, L7

Eggleton, P. P. 1973, MNRAS 163, 279

Grinstein, F. F., Margolin, L. G., &amp; Rider, W. J. (eds.) 2007, Implicit Large Eddy Simulations , Cambridge University Press

Grolsch, A. 2015, in Proceedings of this conference

Herwig, F., Woodward, P. R., Lin, P.-H., Knox, M., &amp; Fryer, C. 2013, ArXiv e-prints

Hurlburt, N. E., Toomre, J., &amp; Massaguer, J. M. 1984, ApJ 282, 557

Kochanek, C. S. 2014, MNRAS

Kolmogorov, A. 1941, Akademiia Nauk SSSR Doklady 30, 301

Kolmogorov, A. N. 1962, Journal of Fluid Mechanics 13, 82

Landau, L. D. &amp; Lifshitz, E. M. 1959, Fluid mechanics , Pergamon

Lorenz, E. N. 1963, Journal of Atmospheric Sciences 20, 130

Maeder, Andr´ e 1999, Physics, Formation and Evolution of Rotating Stars , Springer

Magic, Z., Collet, R., Asplund, M., et al. 2013, A&amp;A 557, A26

Magic, Z., Weiss, A., &amp; Asplund, M. 2014, ArXiv e-prints

Meakin, C. 2015, in Proceedings of this conference

Meakin, C. A. &amp; Arnett, D. 2007, ApJ 667, 448

Meakin, C. A. &amp; Arnett, W. D. 2010, Ap&amp;SS 328, 221

Moc´ ak, M., Meakin, C., Viallet, M., &amp; Arnett, D. 2014, ArXiv e-prints

Prandtl, L. &amp; Tietjens, O. G. 1934, Applied Hydro-&amp; Aeromechanics , Dover Publications, Inc.

Smith, N. &amp; Arnett, W. D. 2014, ApJ 785, 82

Stein, R. F. &amp; Nordlund, ˚ . 1998, A ApJ 499, 914

Viallet, M., Meakin, C., Arnett, D., &amp; Moc´k, M. 2013, a ApJ 769, 1

Vitense, E. 1953, ZfA 32, 135

Wood, T. S., Garaud, P., &amp; Stellmach, S. 2013, ApJ 768, 157

## Discussion

Andr´ e Maeder: Analytical studies suggest that transport processes by meridional circulation and by shear diffusion are strongly affected by the horizontal turbulence in differentially rotating stars. What can the 3D simulations tell us about the possible value of the horizontal turbulence?

David Arnett: Great question! We do not seem to see Jean-Paul Zahn's strongly asymmetric diffusive mixing velocities. However the overturn is effective at mixing the whole convective region, except for Si burning, which has 'cellular' burning regions. The shellular approximation may be fine for earlier stages, up to and including oxygen burning, if we can get the boundary motion-the entrainment rates-right. The shellular approximation will break down during explosive or eruptive events due to 3D instabilities, of course.

Now suppose we spin up a non-rotating star. Turbulence is turbulence, whether driven by buoyancy or by differential rotation. At first we get slight distortions of spherical equipotentials, leading to meridional circulation, and the solenoidal flow also reacts to the inertial accelerations. Eq. 2.8 is a 3D vector equation, and contains these effects; the velocity includes the velocity of meridional circulation and the horizontal turbulence, at least in principle. An interesting issue is the relative importance of dissipation in the boundary layers versus in the bulk of the turbulent flow.

As the spin increases (Rossby number decreases) it will have an increasingly important influence on flow, convection, and angular momentum transport, and the importance of MHD becomes an issue (our code does not have MHD active at present). One can think of a continuous sequence from non-rotating star to accretion disk, parameterized by angular momentum, and the literature contains examples of 3D simulations all along this sequence. That said, the results are complex, and it seems that we still have much to learn.

Victor Khalak: How may stellar rotation be introduced into the equation for turbulence?

David Arnett: It is already in Eq. 2.8 if we do the usual transformation to a rotating frame. This is a deceptively simple vector equation (the Navier-Stokes equation in the turbulent limit), which has deep connections to a lot of theoretical work on rotating stars. Andr´ Maeder has written a beautiful development of the physics of rotating stars e (Maeder, Andr´ e 1999) which starts from the Navier-Stokes equation (his Eq. 1.2).

Ehsan Moravveji: In the very vicinity of the fully mixed core and the radiative layers (braking and/or entrainment layer?), what is the behavior of ∆ ∇ term in your energy balance equation?

David Arnett: Thanks, this is an important point that I went through too quickly. The answer may be implicit in § 2.3 and § 2.4. The ∆ ∇ is a factor in the buoyant acceleration, and can change sign. It depends upon both the temperature gradient and the composition gradient. However mixing makes the composition gradient tend to zero, which changes ∆ ∇ itself; this problem must be solved implicitly. In the simulations there is a transition layer in which the composition does change, but it may not yet be resolved numerically (the 'boundary layer' in Fig. 1). This is the layer where the velocity field changes from solenoidal to potential flow, i.e., from convection to waves. See also the contribution by Arlette Grolsch (Grolsch 2015).

The simulations show additional complexity: the boundary is dynamic and has vigorous and fluctuating wave motion, features not in MLT. If we assume that the boundary dynamics has only a slow secular variation on average (an 'entrainment' model, Meakin &amp; Arnett (2007)), we might use Eq. 2.9 and asteroseismology to try to make progress.
W. David Arnett
3D and Some Other Things Missing from the Theory of Massive Star Evolution
This is a sketch of a 321D approximation for stellar convection which is
nonlocal, and thus has nonzero fluxes of KE (to be published in more detail
elsewhere). Boundary conditions are discussed in a fluid dynamics context
(i.e., predictions for overshoot, semiconvection and entrainment are analyzed).
We plan to add this as an option to MESA. Inclusion of KE fluxes seems to help
resolve the solar abundance problem (Asplund 2009). Smaller cores may ease the
explosion problems with core collapse supernova simulations.