3D Simulations and MLT: I. Renzini's Critique

Renzini (1987) wrote an influential critique of mixing-length theory (MLT) as used in stellar evolution codes, and concluded that three-dimensional (3D) fluid dynamical simulations were needed to clarify several important issues. We have critically explored the limitations of the numerical methods and conclude that they are approaching the required accuracy. Implicit large eddy simulations (ILES) automatically connect large scale turbulence to a Kolmogorov cascade below the grid scale, allowing turbulent boundary layers to remove singularities that appear in the theory. Interactions between coherent structures give multi-modal behavior, driving intermittency and fluctuations. Reynolds averaging (RA) allows us to abstract the essential features of this dynamical behavior of boundaries which are appropriate to stellar evolution, and consider how they relate static boundary conditions (Richardson, Schwarzschild or Ledoux). We clarify several questions concerning when and why MLT works, and does not work, using both analytical theory and 3D high resolution numerical simulations. The composition gradients and boundary layer structure which are produced by our simulations suggest a self-consistent approach to boundary layers, removing the need for ad hoc procedures for 'convective overshooting' and `semi-convection'. In a companion paper we quantify the adequacy of our numerical resolution, determine of the length scale of dissipation (the `mixing length') without astronomical calibration, quantify agreement with the four-fifths law of Kolmogorov for weak stratification, and extend MLT to deal with strong stratification.

Thermonuclear burning in stars gives composition change, which drives stellar evolution. Convection is a fundamentally important process because it (1) moves energy and (2) gives mixing of fuels and ashes. The complete problem is daunting (Arnett & Meakin 2016;Mocák, et al. 2018); at present it is not possible to directly simulate turbulent convection in stars over evolutionary times. This paper is an early step and simplifications are needed. We first consider convection separately from the effects of rotation, binary companions, and magnetic fields, and consider only a sector of a convective shell. This is essentially the classic Böhm-Vitense problem, which is standard for stellar evolution. We also consider turbulence, wave generation, and boundary physics, each of which were not part of the Böhm-Vitense problem, but should have been.
In §2 we summarize Renzini's critique, and refer back to these specific problems in subsequent sections. We use this critique as a foundation on which we illustrate the direct relevance of 3D simulations to stellar evolutionary practice. In §3 we present our methods of simulation and analysis, which allow us to study flows of relatively strong turbulence (Reynolds numbers up to ∼ 10 4 ). In §4 we reject the "blob" picture of MLT, examine the length scale for turbulent damping, the balance between driving and damping, and introduce the time-dependent turbulent kinetic energy (TKE) equation. In §5 we present simulation results for the evolution of the TKE, including multi-modal behavior and wave generation. In §6 we discuss composition gradi-ents and boundary layer structure. 1 In §7 we summarize this paper.

RENZINI'S CRITIQUE
In his classic critique of MLT, Renzini (1987) focused on convective "overshooting", which denotes attempts to deal with boundaries in a merely local theory. In discussing some proposed overshooting algorithms, he identified several fundamental problems with MLT that he labeled as "embarrassing": 1. The ends problem: infinite accelerations and decelerations are required at the beginning and end of the MLT trajectory. MLT does not define a boundary, so that additional physics must be assumed, usually involving the Schwarzschild or the Ledoux linear stability condition.
2. The two lengths problem: are the path length and the size of the "blob" the same?
3. The resolution problem: are lengths resolved which are smaller than the mixing length (or a pressure scale height H P )? Is a convective cell of the order of the zone size?
4. The fluctuations problem: turbulent fluctuations are ignored even if not small, so that MLT does not deal with a "storm of the century" event, nor the accumulated effects of fluctuations.
5. The origin problem: what is the flow pattern near the center of the star? Conservation of baryons requires that any flow into the central regions must be balanced by a flow out; if the radial velocity is non-zero at the origin, then its gradient must be zero.
6. The braking problem: what causes the flow to turn, and be contained in the convection zone? This requires buoyancy braking, which is not in MLT; and is often patched by "overshoot" prescriptions.
7. The dynamics problem: Renzini's "wind and water line" problem, or "boundary layer" dynamics. How are waves generated? How do convective boundaries grow and recede?
8. The non-locality problem: what are the turbulent trajectories, and do distant regions affect local motion?
9. The flux of turbulent kinetic energy problem: this is ignored in MLT but has been demonstrated to be non-negligible for strongly-stratified flows (e.g., surface convection zones in stars).
10. The composition problem: composition is assumed to be homogeneous in MLT, which is violated in regions having nuclear burning, for example.
We will take these problems one by one, illustrating the strengths and weaknesses of MLT, and the improvements suggested by 3D simulations and their theoretical analysis. It has become traditional to complain about the flaws of MLT, but as Renzini (1987) emphasized, MLT works surprisingly well in some respects. We focus on the questions: why does MLT work at all? and what is still missing?
It is appropriate here to warn the reader that two assumptions of MLT are violated by 3D simulations: (1) the turbulent velocity field is nonlocal, and (2) the net radial turbulent kinetic energy flux is not always negligible. Renzini's "embarrassments" are affected by these differences, which modify physics at boundaries.

SIMULATION METHODS
Our set of 3D numerical simulations 2 are "box-in-star" computations which range from "very low resolution" (128 3 zones) to "very high resolution" (now 1536 × 1028 2 and 1536 3 zones). These are implicit large eddy simulations (Woodward 2007;Grinstein, Magolin & Rider 2007;Apsden, et al. 2008, ILES). They include two stable layers sandwiching a turbulent convective region; see Fig 1, which shows a cross-sectional slice through a representative case. Velocity magnitudes are shown (red is high; white is moderate; blue is low). This figure illustrates the complexity of the turbulent flow, and strikingly shows the boundary layers which form at top and bottom of the turbulence. The boundary layers are thin, and not step functions but sigmoids (Cristini, et al. 2017). Coherent structuresrolls and plumes-are strongly dynamic: forming, breaking apart, and re-forming elsewhere. Both boundary layers are also dynamic. They bend and stretch, and radiate g-mode waves (which are most clearly visible at the top). There is intermittency (Tennekes & Lumley 1972) in both space and time, as the "patchiness" in Fig. 1 indicates; intermittency is related to nonlinear interactions of coherent structures (Warhaft 2002). This is confirmed by movies of (1) the evolution in time ("Very high resolution movie of the C-shell"), and of (2) a fly-through of the computations at a given instant in time Figure 1. Vertical slice of a 3D, 1024 3 simulation of a carbon burning shell (Cristini, et al. 2017). Velocity magnitudes are shown (red is high; white is medium; blue is low). This figure illustrates the complexity of the turbulent flow, and strikingly shows the boundary layers which form at top and bottom. The structures-rolls and plumes-are intermittent: forming, breaking apart, and re-forming elsewhere. Both boundary layers are also dynamic, and radiate g-mode waves (most clearly visible at the top).
("Carbon shell (1024 3 ) simulation: fly-through movie"). 3 To the extent that the simulations are in a statistical steady state in time, and statistically homogeneous in space, the movies will have a similar visual appearance (as they do). Turbulence has a space-time isotropy. This allows averaging procedures to be robust.
Unless stated to the contrary, our simulations include a full reaction network of appropriate size for the problem at hand (Arnett 1996), with reaction rates coupled to the thermodynamic and compositional fluctuations, so that nuclear heating and neutrino cooling are dynamic variables. Simplifications are only made after explicit testing of their validity.

Methods: ILES and RA
The simulations extend in resolution from above the integral scale of turbulence, down to well inside the inertial range of the cascade (Cristini, et al. 2017). The full cascade is represented because the method used (PPM, Colella & Woodward (1984)) solves the Riemann problem for nonlinear flow at the individual zone level. This method and others of its class (finite-volume monotonic solvers, Grinstein, Magolin & Rider (2007); Leveque (2002)) automatically result in a Kolmogorov description of the turbulent cascade down to a dissipation scale at a sub-grid level. These ILES methods use Riemann solvers to capture shocks at the scale of a zone ∆r. They relate the change in turbulent kinetic energy (∆v) 2 and traversal time ∆r/∆v across the shock, to the rate of specific entropy production across the shock − ∂S ∂t ∼ (∆v) 3 /∆r. This is the Kolmogorov expression for a turbulent cascade, so these methods automatically (implicitly) match motion to a turbulent cascade The time evolution and the flythrough movies may be found at http: //www.astro.keele.ac.uk/shyne/321D/ convection-and-convective-boundary-mixing/ visualisations. at the grid scale. 4 Conservation of mass, momentum and energy are enforced to machine accuracy, so that numerical error is concentrated in the calculated shape of the flow field, not the conservation laws.
These simulations solve the Euler equations, not the Navier-Stokes equations, so the question of convective instability is determined by the Reynolds number. The formal Rayleigh number may be infinity (no explicit radiative diffusion 5 ). The 1D model used for guidance (Cristini, et al. 2017) had a very high Péclet number (Pe ∼ 10 8 ), but the turbulent cascade implies that the results are insensitive to Pe (see also Orvendahl, et al. (2018)). The convective Mach numbers are small (< 0.02).
For these simulations the effective Reynolds numbers are roughly Re ∼ n 4/3 , where n is the number of zones across the turbulent layer, so Re ∼ 600 to more than 2 × 10 4 . This allows us to compute from the integral scale of fully turbulent flow, down into the inertial range of the cascade 6 , with no additional assumptions concerning flow geometry.
To tame the turbulent fluctuations we integrate over angle (or the horizontal dimension of the convective region), and over several turnover times, to obtain the average behavior over a spherical shell.
A novel and key feature of our procedure is that it avoids the classical closure problem of the Reynolds-averaged Navier-Stokes (RANS) equations (Tritton 1988). The numerical simulations do not have the un-physical nondynamic fluctuations found in unconstrained statistical methods, and therefore provide exact averages, limited only by the granularity of the space-time grid. To be precise we call our method "Reynolds-Averaged ILES", or RA-ILES (Mocák, et al. 2018). The RA-ILES method allows an accurate and separate assessment of dissipation due to (1) turbulence, and (2) resolution error, as we show in the companion paper. Comparison of simulations with different resolution shows the effects of finite zoning. For more detail see Arnett, et al. (2015); Viallet, et al. (2011);Cristini, et al. (2017Cristini, et al. ( , 2018; Mocák, et al. (2018).
The ILES approach is a natural complement to the more traditional method, Direct Numerical Simulation (DNS), which resolves the dissipation scale (Pope 2000) but strains to encompass the larger scales, and generally cannot deal with turbulence at those scales which generate motion in stars. In contrast, ILES resolves these large scales, but approximates by a Kolmogorov cascade the behavior downward to a much smaller dissipation scale. Kolmogorov theory is one of the success stories of turbulence (Frisch 1995), but it too is not perfect as it does not contain the direct influence of the large scales on small scales (Warhaft 2002). In ILES this is accomplished at least in part by the nonlinear interaction of large scales, giving intermittency (Tennekes & Lumley (1972) ;Holmes, Lumley & Berkooz (1996), §3). DNS and ILES have different strengths and weaknesses, so that taken together they provide a fuller picture.
RA-ILES adds to ILES the power to detect errors due to finite granularity in the computational spacetime, and allows us to quantify the relative importance of different physical effects in our problem (Mocák, et al. 2018). This aids in the construction of mathematical models of much lower complexity, whose solutions nevertheless approximate the full numerical simulations. We shall freely use such approximations to illustrate the underlying physics. Thus, our procedure has three components: numerical (ILES), analysis of numerical simulation (RA), and analytic approximation (321D).

Limitations: rotation and MHD
Even a perfect solution to the Böhm-Vitense problem would not solve the more general issue of convection in stars. Featherstone & Hindman (2016) suggest that supergranulation is a rotationallyconstrained flow; to add rotation goes beyond the Böhm-Vitense formulation. Rotation requires a star-in-box approach to capture the largest scale; see pioneering simulations of . Rotation forces non-locality ) and symmetry breaking (Viallet, et al. 2013).
To the extent that magnetohydrodynamics (MHD) is important, the symmetry is broken upon which Kolmogorov (and we) rely; that is an issue for future research 7 . Turbulence makes and stretches vortices (e.g., Pope (2000)); a seed magnetic field will be compressed (doing work against magnetic pressure) and stretched (doing work against magnetic tension); (see Parker 1979;Davidson 2001). Fluid kinetic energy is therefore converted into magnetic field energy, and fluid flow is retarded, giving a dynamo. Magnetic fields are buoyant and will tend to rise in a gravitational field. Unless all the field escapes, it strengthens. Stronger magnetic fields tend to stop the flow, which then no longer generates magnetic field, allowing convection to reassert itself. There is a tendency for a magnetic cycle, reminesent of the solar cycle.
Stars are made of high energy-density (HED) plasma 8 , so that magnetic fields will be ubiquitous in stars, but what geometry, strength and dynamic behavior will they have? Geometry is evidently important: 2D turbulent flow develops a reverse cascade in which strong vortices form, merge, and grow in size, while in 3D such vortices are shredded.
We have chosen the simplest version of this complex problem, that of non-rotating, nonmagnetic Böhm-Vitense convection. We have added turbulence, stratification, non-uniform composition, time dependence, non-locality, and boundary physics, but incorporating rotation and magnetic fields remains a challenge for the future. We regard our RA-ILES method as an interesting approximation which, unlike DNS, can resolve the integral scale for stellar turbulence.

SOME INSIGHTS FROM THE TURBULENT KINETIC ENERGY EQUATION
Simulations of 3D turbulence have a very large number of degrees of freedom; to understand them we will frequently use mathematical models which have solutions similar to the 3D simulations, but a much reduced set of degrees of freedom. In spirit this reduction parallels the method of Holmes, Lumley & Berkooz (1996), but in contrast uses RA-ILES analysis to identify dominant terms ( §3.1). To better understand these results we construct even simpler analytic solutions, more at the level of (but better than) MLT.

Blobs?
One major flaw in MLT is the idea of a convective "blob", which is created, accelerated by buoyancy along a ballistic trajectory of length , and then dissolved into the background. As Renzini (1987) discusses, this causes problems at the beginning and the end of the trajectory (the ends problem), and introduces a free parameter , which is taken to be both the size of the blob and the distance it travels.
No such "blob" behavior is seen in 3D simulations 9 of stellar atmospheres (Stein & Nordlund 1989), or interiors Meakin & Arnett 2007b); instead the behavior is that of a complex and dynamic combination of rolls, plumes and other shapes (Fig. 1). Renzini's problems 1, 2, 3 and 4 in §2 are all related to flaws in this "blob" concept.

Length scale for convective velocity
Baryon conservation and the assumption of a quasi-static, quasi-spherical background place strong general constraints on the nature of convective flow in stars: mass flux balance implies, for a spherically symmetric star, where H ρ is the density scale height 10 defined as −(∂ ln ρ/∂r) −1 . Eq. 2 connects the convective velocity structure to a length scale without any MLT assumption. For a medium of uniform density, H ρ → ∞, and the scale becomes the size of the turbulent medium, which is finite. This does give the characteristic length scale for a quasi-steady flow in a stratified medium (such as stellar convection). However, this equation is linear in velocity, so that the velocity scale is not constrained; another equation is required to determine it (involving convective enthalpy flux, for example). Any form of "rotational" flow (Eq. 1) can remove the ends problem because such flow turns back to remain in a finite volume. Lorenz (1963) showed that the simplest version of such a flow (a 2D convective roll) is an example of deterministic chaos; the flow chaotically flips from clockwise to counter-clockwise. In 3D the angular momentum vector is not restricted to only two orientations (as in 2D), but may wander through 4π steradians 11 . Such instabilities associated with one or several strange attractors seem to act as seeds for turbulence (see §5).

Length scale for turbulent damping
The viscosity ν of stellar plasma is roughly comparable to that of common fluids, but stars are so large that stellar Reynolds numbers 12 are far larger than we encounter terrestrially, making their flows highly turbulent (Arnett & Meakin 2016). Kolmogorov (1962Kolmogorov ( , 1941 showed that K , the rate of flow of turbulent specific kinetic energy through a 3D turbulent cascade, is determined by the large scale flows (Landau & Lifshitz (1959), §31), where v is the speed, d = αH P the linear size of such flows, and the transit time is τ = d /|v|. The flows at such scales are affected by boundary conditions, so that α is not a universal constant, and Eq. 3 may apply only on average (Arnett, Meakin & Young 2009). This corresponds to a drag (negative acceleration) of . The mixing-length assumption provides a representation of the effect of a drag term, but no estimate of its value.

Balance between driving and damping
The chaotic driving of turbulence causes large fluctuations, and requires that we average instantaneous properties to obtain useful 11 See Arnett & Meakin (2011) for a roll model, and also Gabriel & Belkacem (2018) for a plume model. 12 The Reynolds number is the product of the length and velocity scales, divided by the viscosity (Landau & Lifshitz 1959, Eq. 19.1).
variables for stellar evolution ( §2.6 in Arnett, et al. (2015), and Meakin & Arnett (2007b)). Fluctuations are not part of MLT; see the fluctuations problem in §2. We do a double average, over angles (spherical shells), and over several turnover times, which are short compared to evolutionary times 13 . Here ∆r cz is the depth of the convection zone. When performed on even modestly resolved numerical simulations of convection, such averaging shows a balance over the turbulent region, between (1) large scale driving and (2) dissipation at the small-scale end of the turbulent cascade. Convection in strong stratification (∆r cz 2H P ) is also driven by "pressure dilatation" as well as buoyancy (Viallet, et al. 2013).
For weakly-stratified convection zones (and MLT) the buoyancy driving (the work done by buoyant acceleration B acting on the turbulent velocity v) is where the gravitational acceleration vector is g, the super-adiabatic excess is ∆∇ = ∇ − ∇ e , and β T is the thermal expansion coefficient (Kippenhahn & Weigert 1990). The entropy excess ∆∇ may contain contributions from composition differences, which are ignored in MLT. To evaluate these, the issue of mixing must be solved consistently with that of convection, complicating the problem (Arnett, et al. 2015;Woodward et al. 2015;Mocák, et al. 2018). The rate of dissipation for turbulent kinetic energy due to the turbulent cascade is, on av- which is essentially the Kolmogorov value for homogeneous isotropic turbulence 14 (Frisch 1995). The Brunt-Väisälä frequency squared (Eq. 6.17, Kippenhahn & Weigert 1990) is so that we may write an equation whose solutions can be related to the simulation results. For −H P N 2 > 0 we have convective flow which is turbulent. This is basically a statement of Newtonian mechanics, with driving by buoyancy and damping by drag 15 . Gough (1977) gives a historical context going back to Prandtl and to Biermann. The early attempts (and most of the recent ones) have used a kinetic theory model, in which the mixing length was a sort of mean free path. To connect with numerical results, we prefer a model of the momentum equation for fluid dynamics, involving structures such as waves, convective rolls, or plumes.
Taking the dot product of Eq. 9 with v gives a turbulent kinetic energy equation, for which the steady-state solution 16 is a bal-14 For our purposes this may be adequately accurate, but the idea of homogeneous isotropic turbulence may be over-simplified; intermittency and anisotropy observed at small scales may require refinement of these ideas (Warhaft 2002). See §3.1. 15 If N 2 0 we are in the wave regime, so that v/τ must be negligible, and we have the wave equation. For N 2 ≈ 0 but not N 2 0, buoyancy braking may be operative, and Kolmogorov damping happens. This is a surface of separation, where Prandtl theory has a singularity. 16 Care must be taken (for negative v) with the sign of the transit time τ and the deceleration.
ance between driving and damping, with d ≡ 2 M LT /8H P for MLT, and ∆∇ > 0. In Eq. 10, negative values of ∆∇ are allowed; this permits buoyant deceleration. The singularities in MLT at the convective zone boundaries ( §9 in Gough 1977), and in boundary layers ( §40 in Landau & Lifshitz 1959) are removed by the introduction of the total time derivative of the specific turbulent kinetic energy. The singularities in this case occur in Prandtl's equations for a boundary layer as the velocity perpendicular to the surface goes to zero. In a star the motion does not go to zero but becomes wave-like 17 rather than turbulent (e.g., §5.3; Fig. 4).
The flow is relative to the grid of the background stellar evolution model, so the comoving time derivative of turbulent kinetic energy may also be written as where F K = ρv(v · v)/2 is a flux of turbulent kinetic energy. The generation of a divergence of the kinetic energy flux in this way is robust for dynamic models; it occurred in the precise RANS approach as well (Meakin & Arnett 2007b). This flux acts to spread locallydriven turbulence more evenly over the turbulent region, to be dissipated as Kolmogorov suggested. Eq. 10 and 11 together comprise a spartan form of the turbulent kinetic energy equation (Meakin & Arnett 2007b). Our RA-ILES numerical simulations satisfy Eq. 7 and 9 (Arnett, Meakin & Young 2009;Arnett & Meakin 2011) and as well as a local balance on average, of driving by acceleration and deceleration, over a weakly-stratified convection zone (Meakin & Arnett 2007b;Arnett, et al. 2015); 17 Because the Mach numbers are low, gravity (g-)modes dominate over compressional (p-)modes.
phase lags between driving and damping cause the pulses seen in Fig. 2 (below).
Thus, on average, which is something like MLT: see Arnett, et al. (2015) for a discussion which allows buoyancy braking. This equation is nonlinear in velocity and therefore can set a velocity scale for Eq. 2, but with the Kolomogorov damping length d replacing the mixing length.
Use of the turbulent cascade removes the two lengths problem. The damping length is not a size, but a measure of the rate of damping due to turbulence. Further, the cascade resolves all scales down to the dissipation scale of Kolmogorov, so that the resolution problem also disappears.
These results are reminicent of some previous work, e.g., Canuto & Mazzitelli (1991); Canuto, Goldman & Mazzitelli (1996) who developed a theory of full-spectrum turbulence, and Canuto (2011a,b,c,d,e) who made further progress with a Reynolds-stress approach. These earlier suggestions, like the 3D simulations, shift the focus from the original MLT picture of blobs to one involving a turbulent cascade.

EVOLUTION OF TURBULENT KINETIC ENERGY
Kolmogorov dissipation is derived for a homogeneous, isotropic turbulent medium (Frisch 1995), where boundaries are not important. We accept this hint, and first examine the behavior of specific turbulent kinetic energy in bulk, reserving a discussion of boundaries until later ( §6).
How much kinetic energy is involved in the convective motion? This is Renzini's flux of turbulent kinetic energy problem. Fig. 2, top pane, shows the evolution in time of specific turbulent kinetic energy (TKE) in the carbon-burning simulations (Cristini, et al. 2017(Cristini, et al. , 2018, for resolutions of 128 3 , 256 3 , and 512 3 . Multi-mode behavior, like that seen in Meakin & Arnett (2007b) (Fig 4) and Arnett, Meakin & Young (2009) (Fig. 5), is evident. Simulations of 768 3 , 1024 3 , and 1536 3 were not continued for such long times, but they tracked the 256 3 and 512 3 results; the short 1024 3 case is shown in cyan. All simulations are consistent with an approach to a quasisteady state, but with significant and continuing fluctuations around that average value. Such behavior, while typical of turbulent flows, is not in MLT (Renzini's fluctuation problem), but is a feature of high resolution simulations; see also Fig. 5 in Woodward et al. (2015).

Initial transient
Each simulation is initiated from a hydrostatic state, recalculated on the grid to machine accuracy. It is then overlaid with very small random perturbations in density, with amplitudes which were 10 −3 of the initial static value. Convection grows gently in the unstable region, gradually forming a non-linear chaotic flow (the turbulent cascade plus coherent structures). The convective Mach numbers rise to ∼ 10 −2 ; these kinetic energies are 10 4 times larger than implied by the initial seeds. Such small perturbations quickly disappear in the stable regions 18 .
Initially there is no turbulent dissipation, only driving by buoyancy. A large first pulse develops due to this delay in the turbulent cascade 19 . After this pulse, driving balances dissipation on average, but not exactly: there is  (Cristini, et al. 2017(Cristini, et al. , 2018, versus radius for 18,000 seconds. This shows "lrez" (128 3 , red dots) and "mrez" (256 3 , green dashes), and "hrez" (512 3 , solid blue line) simulations. Fluctuations in TKE are seen, as found by Meakin & Arnett (2007b). Multi-mode behavior is evident. The 128 3 simulation shows effects of lower resolution (higher numerical viscosity), but 256 3 and 512 3 seem unaffected. Wave kinetic energy (same colors but thin lines) is much smaller, and is plotted (multiplied by a factor of 25) along the bottom. Higher resolution simulations (768 3 , 1024 3 , and 1536 3 have not been done for such long times, but agree fairly well with the 256 3 and 512 3 cases shown. A small segment in cyan of the 1024 3 simulation is shown for illustration. (bottom two panes) Mode analysis of 256 3 simulation, for frequencies 1.6, 0.23, 1.8, 0.81, and 2.1 mHz. A trend of (5.26 × 10 3 t 2 − 4.41 × 10 7 t + 1.56 × 10 12 ), t in s, was divided out; it corresponds to an evolutionary change (see text). While the pulses are not sine waves, they are periodic, so we separate out a few (5) frequencies that are shorter than the turnover time, as expected for a turbulent cascade. Even five modes begin to capture the "pulses" moderately well (bottom panel). These low order, multi-mode fluctuations are robust features of the TKE. a phase lag (Arnett & Meakin 2011), so that the pulses do not disappear because damping always lags driving. 20 There is an evolutionary growth in TKE due to a slight mismatch of thermal balance between initial conditions based on a 1D MLT model, and the energetically scaled 3D model (Cristini, et al. 2017). The heating was 10 3 of its realistic value; although speeding the evolution correspondingly, the heating per turnover is still very small in comparison to the internal energy. This trend is shown as a line in the middle pane, and is removed in order to extract frequencies of the dominant modes (bottom panes). A few (five) modes are sufficient to capture most of the multi-mode behavior, as suggested by Holmes, Lumley & Berkooz (1996).

Multi-modal behavior
Movies of the simulations support the view that the fluctuations in Fig. 2 are caused by turbulent break-up of multiple 3D rolls. Such multi-modal behavior also causes motion of the convective boundaries, driving waves into neighboring regions, and causing variations in TKE. These pulses are clearly seen in Fig. 2, and do not attenuate noticeably over ∼ 30 turnover times.
Oxygen burning requires no scaling of the heating rate because the consumption of fuel may be explicitly followed. However, the interaction of turbulence, burning, and mixing of Ne and O is more complex than generally realized (Mocák, et al. 2018), so that we simplify at this point by focusing on C burning, which is similar but has no Ne ingestion. On these time scales (2 × 10 4 s), little carbon is consumed, even with the enhanced rate. These carbon burning simulations need no explicit nuclear evolution over this time scale.
The middle pane in Fig. 2 shows the original TKE curve and a flatter one after "detrending" to remove the effect of a slow thermal evolution (see caption). The bottom pane shows the reconstruction for 3, 4, and 5 frequencies, which resembles the "detrended" curve increasingly well; at ten frequencies (not shown) the fit is excellent, but even a few modes are sufficient to capture the basic behavior. A proper orthogonal decomposition should allow a better representation (Holmes, Lumley & Berkooz 1996), but for now we simply want to emphasize the robust nature of the pulses in time. The pulses, which are comparable to a transit time in duration, imply that statistical estimates of turbulent properties will be modulated by multi-mode behavior, unless averaged over many transit times.

Wave generation
The kinetic energy in waves is shown, in the top pane in Fig. 2, as thin lines, scaled up by a factor of 25 for better visibility and using the same color code as for resolutions. Arnett, et al. (2015) showed that the boundary of convection has a particular, dynamically required structure, which implies a particular rate of wave generation. In order for matter to turn and remain in the turbulent zone, outgoing flows must be decelerated. Buoyancy braking requires that the square of the Brunt-Väisälä frequency change sign (Eq. 9). This means waves are supported, with their amplitude depending upon the stellar structure and the vigor of convection. The buoyancy braking provides a direct connection between convective and wave motion. Thus, the physics of a convective boundary requires the generation of waves. In this case the energy in the waves is small relative to that in convective motion. See §6, Fig. 4, where at r < 0.43 × 10 9 cm, the blue curve gives evidence for a background of g-mode waves.

Linear Stability Theory
Fig. 2 contains another implication for stellar physics. In linear stability theory (Aerts, et al. 2010;Unno, et al., 1989) it is not possible to include a realistic treatment of convective driving and damping because these terms are inherently non-local and non-linear (see discussion of the τ -mechanism in Arnett & Meakin (2011)). This issue is sometimes called the "time-dependent convection" problem. Actually, all stellar convection is "time-dependent", to the extent that it based on turbulent flow. Because of intermittency, convection is only "steady-state" in an average sense, if at all. The multi-modal behavior of the simulations is most simply described as an interaction of a few chaotic 3D convective rolls (Lorenz (1963); Holmes, Lumley & Berkooz (1996); Arnett & Meakin (2011), and Fig. 2). Turbulent convection drives waves which, once launched, may then be described by linear theory. A sequence of increasingly stronger waves (increasing convective Mach number) transforms pulsations into explosions.

Resolution
With each new initial resolution, synchronization is lost, so that each simulation is an independent member of a ensemble; all of which are attracted to the cascade in the long term. The RA-ILES integrated values are reproduced with surprising accuracy even at crude resolution (128 3 ), and are well resolved at (256 3 ) and above (Meakin & Arnett 2007b;Viallet, et al. 2013;Arnett, et al. 2015;Cristini, et al. 2017Cristini, et al. , 2018. The lowest resolution case (128 3 ) has the highest numerical dissipation; it settles toward a quasi-steady state with the lowest TKE in Fig. 2. This higher dissipation may also be responsible for the reduced amplitude of the TKE peaks in 128 3 relative to, e.g., the 256 3 simulation. While high resolution and long evolution are both desirable, they are in conflict for a finite computer budget, so the highest resolution runs are relatively short.

COMPOSITION AND BOUNDARIES
The study of turbulent kinetic energy has led us back to the issue of convective boundaries, which is not a part of MLT, but which we now address.

Errors in dynamics at stellar boundaries
Boundaries of stellar convective regions are assumed to be adjacent to regions of convective stability ("radiative" regions assumed to have zero flow velocity). Thus convective flow must be joined to zero flow, which leads to the infinite accelerations and decelerations in MLT (the ends problem of §2).
The treatment of boundaries requires a consideration of dynamics not included in MLT. Because of fluctuations, convective flow at a boundary can only be zero on average. This implies that the boundary moves, and is a source of waves. Pure radial motion gives compression, and p-mode waves. However, non-radial motion is also nonzero, giving shear and gravity-modes. At low Mach numbers the g-modes dominate; only these are visible in Fig. 1.
Terrestrial experiments and numerical simulations show that convective and nonconvective regions are separated by a boundary layer (Fig. 1), which is required to join rotational (∇ × v = 0) and potential (∇ × v = 0) flow patterns (Landau & Lifshitz (1959), §44). Rotational flow is associated with mixing and potential flow is not. Such layers require large gradients to perform the joining of these very different flows, and thus the layers are thin. These layers provide buoyancy braking to turn the flow, which implies negative buoyancy and thus a Brunt-Väisälä frequency which is real (Lighthill (1978), §4.1), so that these layers support waves (see Eq. 9). They are well mixed by the turbulent flow. In such layers convective fluctuations must be coupled to wave production. Beyond these layers, in the radiative zone which may support a composition gradient, the Brunt-Väisälä frequency is also real, so internal waves are supported, and directly coupled to the boundary layer fluctuations. Because they can carry vorticity, g-mode waves can induce some energeticallylimited mixing beyond the boundary layer.
"One of the properties of the region of rotational turbulent flow is that the exchange of fluid between this region and the surrounding space can occur in only one direction. The fluid can enter this region from the region of potential flow, but can never leave it." (Landau & Lifshitz 1959, §34). Mixing, which increases entropy, gives the uni-directional nature of the exchange. Such entrainment of material from a radiative zone is a necessary aspect of convective boundaries, and one not included in MLT. This does not preclude receding convective regions, which could shed regions of decaying turbulence; also see Holmes, Lumley & Berkooz (1996).

Distinct boundary regions
We construct a cartoon to clarify discussion of the average properties; see Fig. 3 which shows the boundary of a convection region (convective to the left, radiative to the right). There are five distinct regions, only two of which (1 and 5) are recognized in MLT: 1. Convection. Inside r 1 , the superadiabatic excess is positive, as are the buoyant acceleration and the enthalpy flux. In this region there is a balance between buoyant driving and turbulent dissipation ( §4.4). This is the "convective" region of MLT (Böhm-Vitense 1958).
2. Braking. Between r 1 and r 2 the superadiabatic excess becomes negative, and the flow can turn, and reverse direction. This is the "braking" region which is required dynamically. It has a negative enthalpy flux, is well mixed, and generates waves. This region does not exist in MLT, and is discussed in Arnett, et al. (2015).
3. Shear. Between r 2 and r 3 is a boundary layer such as seen in Fig. 1. As horizontal components of the velocity grow at the expense of radial components (baryon conservation), the boundary layer develops high shear, and is subject to KH instability (Drazin 2002).
4. Composition gradient. Between r 3 and r 4 is the region that can develop a composition gradient, perhaps maintained at the margin of KH instability (Drazin 2002). This is not part of the original MLT (Böhm-Vitense 1958); it is the composition problem of §2. A "fix" must be added (e.g., "semi-convection" and/or "overshoot" regions). This has led to a search for algorithms to extend MLT, whose justification is that they are at least empirically desirable if not always self-consistent.
5. Radiative. Beyond r 4 is the "radiative" region of MLT. Unlike MLT, waves are automatically implied by flexing of the boundary layers (r 1 to r 4 ).
This crude cartoon is to be understood as a snapshot, which flexes and bends; its average properties are to be identified with stellar properties. The layers between r 1 and r 4 may be relatively thin in radial extent (see Fig. 4). Regions 2 and 3 are not defined within MLT; they are part of the ends problem of §2. In part region 4 is not defined within MLT because MLT assumes homogeneous composition (Böhm-Vitense 1958) braking marginal shear waves Figure 3. Two-dimensional schematic of the average structure of a convective boundary. The length b corresponds to the radius of curvature needed to reverse (contain) the flow (vr → −vr). The centrifugal acceleration is provided by negative buoyancy. The radial direction is denoted by r and the transverse by h. The boundary layer lies between r2 and r3. Imagine the whole system undulating due to turbulent fluctuations.
physics must be assumed concerning composition gradients to bridge this gap. This region is often treated as "semi-convective", or partially mixed by "overshoot". Regions 2, 3, 4, and 5 all support waves and relate to the wind and waterline problem ( §2). Any improvement over MLT will certainly affect regions 2, 3 and 4.
6.3. Composition profile Fig. 4 shows a composition profile of 16 O at the lower boundary of the oxygen burning shell. To the extent that they are resolved the 3D ILES simulations automatically produce dynamically self-consistent boundary physics for turbulent flow.
The rms velocity profiles are shown for radial (red) and horizontal (blue) motion. The dominance of horizontal velocity (blue) for r < 0.43×10 9 cm indicates g-mode waves. The boundary dynamics are those discussed in Arnett, et al. (2015). The need for a turning flow requires that the radial velocity (red) approach zero faster than the horizontal velocities (blue), which have a peak near the boundary. There is a well-mixed region in which this braking and turning occur, and which is subadiabatic. Here the Brunt-Väisälä frequency is real; waves are generated by convective fluctuations and propagate from this region (see §5.3, and Fig. 14, Woodward et al. 2015). The asymmetric behavior of convective boundaries, similar to that discussed in Gabriel, et al. (2014); Montalbán, et al. (2013), can naturally arise in this braking region (0.43 × 10 9 < r < 0.44 × 10 9 cm in Fig. 4). It is a consequence of carefully joining turbulent flow to a radiative region.
The sigmoid shape of the boundary is strikingly similar to that of the averaged mean flow of Garaud, Gagnier & Verhoeven (2017), obtained with DNS for a shear flow setup; see their Figure 3 (our use of the bottom boundary requires a flip in orientation). Arnett, et al. (2015) showed that a turning plume would have such a shearing interface, even without differential rotation. The boundary conditions for differential rotation and convection seem to have a deep family resemblance. Cristini, et al. (2017Cristini, et al. ( , 2018) obtain a similar "sigmoid-like" profile for a carbon-burning shell. This shape from simulations has striking similarity to that inferred from asteroseismology (Arnett & Moravveji 2017) for the outer boundaries of convective cores of red giants. At first sight this is puzzling; radiative diffusion has a very minor role in carbon and oxygen burning, but is more important in red giants (Viallet, et al. 2013). It suggests that the shape of the composition gradient may be insensitive to the effects of radiative diffusion. A clue appears in Fig. 4, which shows a significant shear due to the horizontal velocities (blue curve), which may drive Kelvin-Helmholtz (KH) instabilities . In the low Mach-number limit, this is a kinetic instability, driven just by the shape of the horizontal velocity in the radial direction (Drazin (2002), see Fig. 8.3), and would be insensitive to radiative diffusion. See §3.3 in Woodward et al. (2015) for further discussion.
Consider a global perspective: as nuclear burning proceeds the composition gradient steepens until KH instabilities induce some mixing to flatten it again. The gradient in velocity profile evolves toward neutral stability to KH, and lies near the composition profile that it sculpts, as shown in Fig. 4.
An increase in entropy inside the convective zone tends to increase the size of the turbulent region. However, the amount of increase may be sensitively dependent upon the gradients of thermal and compositional entropy in the boundary regions, as well as the turbulent velocity. A decrease in entropy gives a tendency to shed turbulent layers, which dissipate when no longer driven. Again there may be a sensitive dependence upon the gradients Figure 4. Comparison of the 16 O composition profiles (black), averaged over angle, for different resolutions ("med-res" (384 × 256 2 , dots) and "hi-res" cases (768 × 512 2 , heavy dots) of (Viallet, et al. 2013), and a "very-hi-res" case (1536 × 1024 2 , heavier dots, the "Perth" simulation, see Arnett, et al. (2015)), at the bottom boundary. This shape has striking similarity to that inferred from asteroseismology (Arnett & Moravveji 2017) for convective cores of red giants. See also discussion of the C+C shell in Cristini, et al. (2017Cristini, et al. ( , 2018. Overdrawn are the rms velocity profiles averaged over a spherical shell; red is radial and blue is for the two horizontal directions. Instantaneous profiles are not exactly spherical or constant but dynamic. The velocities do not go to zero outside the convective zone (r < 0.43 × 10 9 cm ) because of significant but small wave motion (see Fig 1).
of thermal and compositional entropy in the boundary regions, as well as the turbulent velocity. This is essentially the Richardson criterion for stability for a layered system; (Turner 1973, §10.2

.3),
But what determines the velocity field at the interface? There is a similarity to the previous discussion ( §3.1) of the turbulent cascade and ILES. Approaching the convective boundary resembles approaching the Kolmogorov scale in velocity space; near the boundary only the small scales can "fit". Turbulence implies a shear layer at the interface.

CONCLUSION
We have systematically re-examined the unphysical aspects of MLT in 1D stellar evolution, which were summarized in (Renzini 1987). We have (1) performed 3D turbulent implicit large eddy simulations (ILES, Grinstein, Magolin & Rider 2007) of the classical Böhm-Vitense problem of stellar convection (Böhm-Vitense 1958), (2) applied Reynoldsaveraging (RA) to these numerical results, and (3) developed several analytic approximations to illustrate the physics involved. The RA-ILES combination is "exact" to the word length of the computer used, unlike the Reynolds-averaged Navier-Stokes equations (RANS) which are not closed (Tritton 1988). This advantage is a consequence of the sub-grid behavior of RA-ILES, which mimics a Kolmogorov turbulent cascade, allowing both stellar size scales and turbulent scales to be described.
The evolution of the turbulent kinetic energy is chaotic, involving a few dominant modes (unstable 3D rolls, Lorenz 1963). Driving and damping of turbulent convection are out of phase; this makes interpretation of time aspects of numerical convergence a challenge, but provides insights into the dynamics of turbulence and stellar variability. We show an example of such turbulent pulses which are not significantly damped in 30 turnover times; stellar convection can attain a dynamic steadystate, with chaotic fluctuations.
The RA-ILES simulations automatically provide the velocity structure at boundaries, including composition profile, structure of deceleration region, and the basis for prediction of wave generation. ILES allows us to describe a turbulent boundary layer, which uses the KH instability to sculpt a composition profile similar to that inferred from astereoseismology (Arnett & Moravveji 2017). Entrainment of stable matter occurs at these dynamic boundary layers (Meakin & Arnett 2007b). The process of turning the flow at the boundary also insures that waves are generated, connecting the convective motion to the wave generation.
We will quantify the resolution errors due to zoning in the companion paper (Arnett, et al. 2019), and predict the turbulent dissipation length (the 'mixing length') from the simulations.