Issue 
A&A
Volume 604, August 2017



Article Number  A125  
Number of page(s)  15  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201630362  
Published online  25 August 2017 
Extreme value statistics for twodimensional convective penetration in a premain sequence star
^{1} Astrophysics, College of Engineering, Mathematics and Physical Sciences, University of Exeter, EX4 4QL, Exeter, UK
email: janepratt@utexas.edu
^{2} École Normale Supérieure de Lyon, CRAL (UMR CNRS 5574), Université de Lyon 1, 69007 Lyon, France
^{3} Monash Centre for Astrophysics (MoCA) and School of Physics & Astronomy, Monash University, Clayton Vic 3800, Australia
^{4} MaxPlanckInstitut für Astrophysik, Karl Schwarzschild Strasse 1, 85741 Garching, Germany
Received: 26 December 2016
Accepted: 12 June 2017
Context. In the interior of stars, a convectively unstable zone typically borders a zone that is stable to convection. Convective motions can penetrate the boundary between these zones, creating a layer characterized by intermittent convective mixing, and gradual erosion of the density and temperature stratification.
Aims. We examine a penetration layer formed between a central radiative zone and a large convection zone in the deep interior of a young lowmass star. Using the Multidimensional Stellar Implicit Code (MUSIC) to simulate twodimensional compressible stellar convection in a spherical geometry over long times, we produce statistics that characterize the extent and impact of convective penetration in this layer.
Methods. We apply extreme value theory to the maximal extent of convective penetration at any time. We compare statistical results from simulations which treat nonlocal convection, throughout a large portion of the stellar radius, with simulations designed to treat local convection in a small region surrounding the penetration layer. For each of these situations, we compare simulations of different resolution, which have different velocity magnitudes. We also compare statistical results between simulations that radiate energy at a constant rate to those that allow energy to radiate from the stellar surface according to the local surface temperature.
Results. Based on the frequency and depth of penetrating convective structures, we observe two distinct layers that form between the convection zone and the stable radiative zone. We show that the probability density function of the maximal depth of convective penetration at any time corresponds closely in space with the radial position where internal waves are excited. We find that the maximal penetration depth can be modeled by a Weibull distribution with a small shape parameter. Using these results, and building on established scalings for diffusion enhanced by largescale convective motions, we propose a new form for the diffusion coefficient that may be used for onedimensional stellar evolution calculations in the large Péclet number regime. These results should contribute to the 321D link.
Key words: methods: numerical / convection / stars: interiors / stars: lowmass / stars: evolution
© ESO, 2017
1. Introduction
Convection is a fundamental stellar process which underlies heat transport, mixing, shear, and the dynamo. During many phases of stellar evolution, a star’s interior is characterized by layers that are convectively unstable. Outside of a convection zone are layers that are stable to convection. Convective fluid motions can cross the boundary between a convectively unstable region and a neighboring stable region. When the Péclet number is large, these convective motions slowly erode the density and thermal stratification of the stable region. When this occurs, the process is termed convective penetration. We study convective motions that may contribute to convective penetration below a convection zone, and thus can have an impact on the structure and evolution of a star.
Studies of stellar evolution utilize onedimensional calculations that evolve physical quantities as a function of the radial position interior to a star. The effects of convection on other physical quantities are typically modeled using stellar mixing length theory, which depends on the local temperature gradient (e.g., Vitense 1953; BöhmVitense 1958; Abbett et al. 1997; Trampedach 2010; Brandenburg 2016). The phenomenon of convective penetration has been included in mixing length theory in two competing ways, based on different nonlocal convection theories. One is an overshooting length, expressed as a fraction of the pressure scale height at the convective boundary, which defines a region outside of the convection zone where the stellar model assumes full convective mixing. The other is a region of enhanced onedimensional diffusivity (Noels et al. 2010; Zhang 2013), where the intensity of convective mixing changes with radius. Using statistical observations, this work addresses both of these methods and thereby contributes to the 321D link.
Stellar convection is a nonlinear and nonlocal phenomenon. The extent to which convective motions can extend beyond the theoretical boundary of a convection zone is intuitively linked to the local velocity, density, and temperature of the fluid at the boundary^{1}. In this work we will refer to convective flow structures that possess radial velocity as plumes, without invoking geometry or physics specific to plume studies. Plumes frequently move a short distance beyond the convective boundary. Less frequently, plumes penetrate the stable zone to a larger extent. These plumes that extend deeper clearly have a greater effect on the star, because they penetrate with larger momentum and affect the stratification deeper inside the stable zone. Convective penetration is therefore naturally a process in which the intermittency is significant. A theoretical description of such flows should be built upon statistics gathered by observing convective penetration over long times. The Multidimensional Stellar Implicit Code (MUSIC), which uses implicit time integration, allows us to cover efficiently the long times necessary for converged statistics on convective penetration. This work thus builds on earlier explorations (e.g., Hurlburt et al. 1986, 1994; Brummell et al. 2002; Rogers & Glatzmaier 2005b; Rogers et al. 2006), which simulate convective overshooting or penetration over shorter times, typically between two and five convective turnover times.
In addition to observing the longtime behavior of convective penetration, a statistical treatment is necessary to model meaningfully instances of large convective penetration. Extreme value theory (EVT) is a statistical theory that has been widely used to predict events that have large impacts. Examples are record rainfall, floods, dangerously high winds, earthquakes, and avalanches. For the case of convective flows penetrating a stable zone, an analysis using extreme value theory is natural.
This work is structured as follows. In Sect. 2 we discuss the young Sun model and the numerical framework of our simulations. In Sect. 3 we examine convective penetration based on simulation data that covers a range of times, and a range of different spatial resolutions and extents. Motivated by these results, in Sect. 4 we introduce the framework of extreme value theory, and compare statistical results for simulations in a truncated domain that model penetration due to local convection, and simulations in a large domain that model penetration due to nonlocal convection. In Sect. 5 we build on these statistical results to formulate a onedimensional radial diffusion coefficient that is relevant in the large Péclet number limit. In Sect. 6 we summarize and discuss the implications of these results.
2. Simulations
In this work we examine a model of a prototypical premain sequence star, the young Sun, using the numerical setup benchmarked in Goffrey et al. (2017). In Pratt et al. (2016) the structure and properties of the stellar model of the young Sun are described in detail. Here we briefly summarize the salient features of our simulations. We perform twodimensional implicit large eddy simulations (LES) of the young Sun using the MUSIC code. The stellar radius of our young Sun model is approximately three times larger than the current Sun, it is one solar mass, and has homogeneous chemical composition. The radial profiles of density and temperature for the young Sun model are typical for a premain sequence star that is no longer accreting and is gradually contracting. The luminosity is increasing with the interior radius. Based on the radial entropy profile and an evaluation of the Schwarzschild criterion, a central radiative zone is expected below a large convection zone. This convection zone spans the outer 1.2 × 10^{11} cm of the total radius of 2.13 × 10^{11} cm. A layer exists between the convection zone and the radiative zone, where mixing of convective flows into the radiative zone can take place. This layer is the focus of our work, and is located far from the physically complicated nearsurface layers. In this layer, the stiffness (e.g., Hurlburt et al. 1994; Brummell et al. 2002; Rogers et al. 2006) is approximately S ≈ 13. Our simulations do not include additional stellar physics such as rotation, a tachocline, chemical mixing, or magnetic fields, that we intentionally omit from our statistical study of convective penetration. This simplifies the dynamics at the boundary between the interior radiative zone and the convection zone, where we observe penetrative convective motions. These additional effects may influence the results we obtain. In three dimensions the shape of plumes may be different from two dimensions (Schmalzl et al. 2004; van der Poel et al. 2015), and this could affect mixing properties. Rotation and shear may tilt plumes, modifying the angle of penetration (e.g. Brun et al. 2017). A magnetic field may shift the Schwarzschild criteria (Gough & Tayler 1966) and fundamentally change the mechanisms of plume braking (Hughes & Proctor 1988). These additional physical effects will be explored in future work.
Parameters for twodimensional compressible hydrodynamic simulations of the young Sun.
The MUSIC code solves the inviscid compressible hydrodynamic equations for density ρ, momentum ρu, and internal energy ρe: using a finitevolume method, a MUSCL scheme for interpolation (Van Leer 1977), and a van Leer flux limiter (Van Leer 1974; Roe 1986). Time integration in the MUSIC code is fully implicit, and uses a Jacobian free NewtonKrylov (JFNK) solver (Knoll & Keyes 2004) with a variant of physicsbased preconditioning (Viallet et al. 2016). For all simulations examined in this work, an adaptive timestep maintains a general CFL number of ten, while a CFL number based on simple advection is restricted to be 0.5. We find that this produces good convergence of relevant basic quantities, such as the average kinetic energy.
MUSIC simulations are designed to contribute to the 321D link (Arnett 2014; Arnett & Meakin 2009), the effort to improve onedimensional stellar evolution models by studying critical short phases using stellar hydrodynamics in two and three dimensions. One aspect that makes our simulations relevant to the 321D link is the use of an equation of state and realistic opacities standard in onedimensional stellar evolution calculations. Opacities are interpolated from the OPAL (Iglesias & Rogers 1996) and Ferguson et al. (2005) tables, which cover a range in temperature suitable for the description of the entire structure of a lowmass star. The compressible hydrodynamic equations, given by Eqs. (1)−(3), are closed by determining the pressure p(ρ,e) and temperature T(ρ,e) from a tabulated equation of state for a solar composition mixture. This equation of state accounts for partial ionization of atomic species by solving the Saha equation, and neglects partial degeneracy of electrons; it is suitable for the description of our solar model at a young age. The initial state for MUSIC simulations is produced using data extracted from a onedimensional model calculated using the Lyon stellar evolution code (Baraffe & El Eid 1991; Baraffe et al. 1997, 1998), which uses the same opacities and equation of state implemented in MUSIC. In Eq. (2), g is the gravitational acceleration, a sphericallysymmetric vector identical to that used in the stellar evolution code, and not evolved in our simulations.
The thermal diffusivity in MUSIC is also realistic for the young Sun model, and has not been enhanced. In Eq. (3)the thermal conductivity is defined using the StefanBoltzmann constant σ and the Rosseland mean opacity . An artificially high thermal diffusivity is a key issue that has been found to cause differences between mixing length theory and numerical results (Rempel 2004), and can result in structural changes for a star. Despite these consequences, enhancing the thermal diffusivity has become a conventional step in numerical studies of convective overshooting (Browning et al. 2004; Rogers et al. 2006; Tian et al. 2009; Strugarek et al. 2011), because it allows the star to more rapidly achieve energy balance by shortening the KelvinHelmholtz time scale. For the young Sun model considered here, the KelvinHelmholtz time scale is several orders of magnitude greater than its convective turnover time. Thus although the young Sun model has not achieved energy balance, during the time that our simulations span, there is no measured change in the thermal profile, and the convection may be considered to be in steady state.
When a realistically small thermal diffusivity is used, a high resolution is typically required to resolve the corresponding smallscale temperature fluctuations. No attempt is made to fully resolve the small scales related to the realistically small thermal diffusivity that we use; we rely on the implicit large eddy simulation paradigm to approximate the effect of this smallscale dissipation. Our use of the realistic thermal diffusivity for the young Sun model allows us to study convective penetration because it produces a large Péclet number; this is discussed fully in Sect. 2.2. In MUSIC, numerical truncation errors contribute to both the viscosity and thermal diffusivity. Because MUSIC simulations also use an explicit thermal diffusivity related to the thermal conductivity in Eq. (3), the Prandtl number is expected to be everywhere less than one. The characteristic length scales, velocities, viscosity, and thermal diffusivity vary throughout the radius of a star. The viscosity and thermal diffusivity may also be dependent on local dynamics, the wavenumber, and the estimation method (Domaradzki et al. 2003; Zhou et al. 2014; Radice et al. 2015; Schranner et al. 2015; Zhou & Thornber 2016). For these reasons, the Rayleigh number and Reynolds number are not specified in a general sense for our simulations.
2.1. Sphericalshell geometry and boundary conditions
Equations (1)−(3)are solved in a twodimensional spherical shell using spherical coordinates: radius r and colatitude θ. A twodimensional geometry has been chosen in order to allow us to produce a long time series of data that also has a satisfactory resolution near the bottom of the convection zone. Convective penetration is the result of largescale convection, and our simulations do not model the significantly threedimensional effects of turbulence, rotation, or magnetic fields present in a realistic star. Therefore a twodimensional simulation provides an approximation to threedimensional stellar convection. It has long been recognized that twodimensional stellar simulations result in higher velocity magnitudes than threedimensional simulations (Muthsam et al. 1995; Meakin & Arnett 2007). However, because our simulations do not approach the high velocities realistic to a star, the difference between two and threedimensional velocities do not seriously detract from a study of the form of convective penetration.
Fig. 1 Typical time snapshots of radial velocity in simulation YS2 (left), YS5 (middle), and YS6 (right). The zero point is colored black, while inward flows are gray and outward flows are red. 
The placement of boundaries and the choice of boundary conditions affect the physical outcome of a hydrodynamic simulation. We use boundary conditions targeted to maintain the original radial profiles of density and temperature of the onedimensional stellar evolution model of the young Sun. Each simulation volume begins at 20° from the north pole, and ends 20° before the south pole, allowing for a wide convection zone with multiple descending and ascending plumes. Periodicity is imposed on all physical quantities at the boundaries in θ. Radial boundaries require more sophisticated treatment. In velocity, we impose nonpenetrative and stressfree boundary conditions on the radial boundaries. The energy flux and luminosity are held constant at the inner radial boundary at the value of the energy flux at that radius in the onedimensional stellar evolution calculation. Typically we hold the energy flux and surface luminosity constant on the outer radial boundary as well; the surface luminosity of our model for the young Sun is 2.32 times the luminosity of the current Sun. An alternative is to allow the surface to radiate energy with the local surface temperature. In this case the energy flux varies as where σ is the StefanBoltzmann constant and T_{s}(θ,t) is the temperature along the surface. This boundary condition can only be effectively used when the temperature gradient near the surface is sufficiently resolved; otherwise it results in artificially high cooling rates. On the inner radial boundary of the spherical shell, a constant radial derivative on the density is imposed, as discussed in Pratt et al. (2016). At the outer radial boundary a hydrostatic equilibrium boundary condition (GrimmStrele et al. 2015) that maintains hydrostatic equilibrium by assuming constant internal energy and constant radial acceleration due to gravity in the boundary cells is applied to the density.
Consider the suite of twodimensional simulations of sphericalshell convection summarized in Table 1. For all of these simulations, the table lists radial grid spacing and the radial extent of the spherical shell in units of the total stellar radius R. Rogers et al. (2006) report using a grid of N_{r} × N_{θ} = 2048 × 1500, with 620 grid cells in the radiative zone to study convective overshooting in the current Sun. By comparison, in simulation YS5 we use a grid of N_{r} × N_{θ} = 2432 × 2048 with 922 grid cells in the radiative zone to study convective penetration in the young Sun. The large section of the radiative zone simulated at high resolution is useful to precisely resolve convective penetration in the young Sun, which has a larger convection zone than the current Sun.
Table 1 provides the parameters of simulations YS06, which consist of three different sets of simulations. Simulations YS02 model convection in the local approximation, using a spherical shell of limited radial extent around the penetration layer. These simulations are designed to examine penetration due to convection that arises and is dynamically limited to the layer around the bottom of the convection zone, similar to a boxinastar approach. Simulations YS35 model convection in the nonlocal approximation, using a spherical shell that includes nearly the whole star, terminating in the nearsurface layers. These simulations are designed to examine penetration due to convection that is dynamically accurate to the large convection zone of the young Sun model, where convective dynamics may interact nonlocally across more than half of the stellar radius. Among each of these sets of simulations, three different resolutions are examined for otherwise identical physical setups and for boundary conditions that include a constant energy flux at the outer surface. In MUSIC simulations, which do not include an explicit viscosity, higher resolution produces higher characteristic velocities in the convection zone. Thus these different resolution simulations are performed to show the robustness of our statistical results with respect to convective velocities. An additional simulation YS6 is performed in order to isolate the effect of allowing the energy flux at the stellar surface to vary according to the local temperature. A typical snapshot of the radial velocity in each of these simulation volumes is shown in Fig. 1. The effect of the placement of the simulation boundaries and the boundary conditions imposed on the properties of convection was examined in Pratt et al. (2016). In this work we examine exclusively the penetration of plumes below the convection zone.
For each simulation, Table 1 provides a time span in units of the convective turnover time of that simulation. This time span indicates a period after steady convection has developed, during which the global kinetic energy of the system is approximately constant. The average radial structure of the star’s density and temperature does not evolve over this period. All timeaverages and other statistics are calculated over this time span.
2.2. Characterization based on the Péclet number
The Péclet number is a nondimensional ratio that measures the relative importance of two physical effects that determine the deceleration of a fluid element entering a layer that is stable to convection: advection and thermal diffusion. Therefore it is standardly used to distinguish between convective overshooting and convective penetration. Low Péclet numbers (Pe ≤ 1) indicate that thermal diffusion is the dominant process and convective motions overshoot the convection zone without significantly changing the stratification of the stable layer; high Péclet numbers indicate penetration can take place, a process where convective fluid motions modify the thermal and density stratification in the stable layer outside the convection zone (e.g., Zahn 1991). The Péclet number is typically defined Here the pressure scale height h_{p} represents a characteristic length scale, and the rootmeansquare velocity v_{RMS} represents a characteristic convective velocity. The thermal diffusivity α is produced from our tabulated equation of state, using the definition of the thermal conductivity χ in Eq. (3), and the specific heat capacity at constant pressure c_{P}. Numerical truncation errors due to the finite resolution of a simulation can increase this thermal diffusivity.
In our simulations, radial profiles of temperature and thermal diffusivity are taken directly from the model of the young Sun produced by a onedimensional stellar evolution calculation. Because of this, and because we use a realistically small thermal diffusivity, the Péclet number is much greater than one in the area of interest at the bottom of the convection zone, and convective penetration is expected. Figure 2 shows the average radial profile of the Péclet number in simulation YS4. Neglecting the contributions to thermal diffusivity from numerical truncation, the Péclet number is – in the layer at the bottom of the convection zone for all of the simulations in Table 1. By perturbing the background state we are able to measure an approximate effective thermal diffusivity. We find that while numerical truncation errors increase the effective thermal diffusivity beyond the value implied by the explicit physical thermal diffusivity, our simulations are still in the Pe ≫ 1 regime. A Péclet number larger than one is expected; Péclet numbers are generally higher in low and medium mass stars, like the young Sun model, than for higher mass stars (e.g. Meynet & Maeder 2000). Péclet numbers are also expected to be much higher in stellar interiors than in the nearsurface layers. In addition to these physical points, twodimensional hydrodynamic simulations are known for producing somewhat higher Péclet numbers because they exhibit higher velocities than three dimensional simulations (Muthsam et al. 1995; Meakin & Arnett 2007).
Fig. 2 Radial profile of the timeaveraged and horizontallyvolumeaveraged Péclet number, as defined by Eq. (4), from simulation YS4. The contribution to the thermal diffusivity by numerical truncations are not included in this profile. The heavy vertical black line marks the boundary between the stable radiative zone and the convection zone, determined from the radial profile of entropy and the Schwarzschild criterion produced by the onedimensional stellar model. 
Based on the Péclet number, we define convective motions that cross beyond the boundary of the convection zone to be convective penetration. Indeed in test simulations where the KelvinHelmholtz time scale of the young Sun is artificially lowered by the use of an enhanced thermal diffusivity and corresponding increased luminosity, we observe penetrative motions that change the density and thermal stratification of the star. However in the simulations presented in this work, the KelvinHelmholtz time scale is many orders of magnitude larger than the characteristic time scale associated with convective motions^{2}. Therefore we do not observe erosion of the stellar structure due to convective penetration; the temperature and density profiles are approximately constant during our simulations. We do study penetrative flows that have the potential to change stellar structure. However our study may also have implications for convective overshooting motions in the low Péclet number regime.
Over the last decade it has become common to refer to the mixing at the boundary of a convection zone as turbulent entrainment (e.g. Meakin & Arnett 2007). We note that convective penetration has long been considered a type of entrainment in which mixing is due to largescale convective motions (e.g. Tennekes & Driedonks 1981), and that in a physical star, mixing at a boundary may be due to a combination of fluid effects including convection, turbulence, and also shear (e.g. Jonker et al. 2013). Generally in LES of stellar convection, the range of scales that correspond to turbulent mixing include scales considerably smaller than the grid spacing of the simulation. Indeed for all of the simulations considered in this work, flows are expected to be in the laminar regime with moderate Reynolds numbers, and the effect of turbulent entrainment is not targeted.
3. Results: the penetration depth
3.1. Determination of the extent of convective penetration
Historically several different definitions have been used for the depth of convective penetration. Broadly speaking, a physical quantity is sought that either drops to a low level, or changes sign where convective motions cease. To measure the penetration depth using a quantity that drops to a low level, an arbitrary percentage of the quantity is chosen. An example of such a quantity is the kinetic energy density . In several earlier works, the point where convective motions cease has been defined as the point where the kinetic energy density reaches 1% (Brummell et al. 2002), or 5% (Rogers et al. 2006) of its peak value near the bottom of the convection zone.
In our simulations, difficulty with this type of measurement arises after convective plumes interact with the underlying radiative region. Convective penetration and overshooting generate a spectrum of waves in the radiative zone that propagate in the angular direction and feed back on the convective plumes. These waves have been identified as internal gravity waves in many cases (Hurlburt et al. 1986; Andersen 1994; Dintrans et al. 2003; Rogers & Glatzmaier 2005a; Rogers et al. 2006; Lecoanet & Quataert 2013; Shiode et al. 2013; Brun et al. 2013; Alvan et al. 2014, 2015; Pinçon et al. 2016). We observe waves in the radiative zone evolving over tens of convective turnover times, and they spoil an evaluation of a percentage of the peak convection zone kinetic energy density. This problem is clear from the longtimeaverage of the kinetic energy density shown in Fig. 3 for simulation YS3. In this longtime average, these waves appear as a small peak in the kinetic energy density at a radius 0.35 <r/R< 0.4, precisely in the region where plumes that have a large penetration depth terminate. A second issue arises when we observe steady convection over tens of convective turnover times: the peak value of the kinetic energy density near the bottom of the convection zone evolves. Thus a criterion based on a percentage of the peak value changes during the long time span of our simulations.
Fig. 3 Radial profile of the timeaveraged and horizontally volumeaveraged kinetic energy, , scaled to its maximum, from simulation YS3. The heavy vertical black line marks the convective boundary. 
For these reasons it is advantageous to base the penetration depth on a quantity that changes sign at the point in the radiative zone where convective motions cease. Two such quantities have been proposed (as discussed in, e.g. Hurlburt et al. 1986, 1994; Ziegler & Rüdiger 2003; Rogers et al. 2006; Tian et al. 2009; Chan et al. 2010): the vertical kinetic energy flux and the vertical heat flux. We examine both of these measures, and compare the results. The vertical kinetic energy flux F_{k} is defined as (6)where ρ is the density, and v is the velocity. The penetration layer can be defined by the region of positive vertical kinetic energy flux. A radial profile of the vertical kinetic energy flux, calculated by averaging both in time and using a volume weightedaverage in the horizontal (θ) direction, is shown in Fig. 4 for simulation YS4. Previous works have used this type of averaged profile to determine the penetration depth. In Fig. 4, the vertical kinetic energy flux exhibits a smooth positive peak; the width of the layer where convective mixing takes place is well defined, but small. In any timesnapshot, larger penetration may be identified at particular angles.
Fig. 4 Radial profile of the timeaveraged and horizontally volumeaveraged vertical kinetic energy flux F_{k} defined in Eq. (6), and vertical heat flux F_{c} defined in Eq. (7)for simulation YS4. Both of these fluxes have been normalized to their absolute maximum value. The heavy vertical black line marks the convective boundary. 
The second criteria used to measure the penetration depth is the vertical heat flux: (7)where the temperature perturbation δT is defined relative to a time and horizontally volumeaveraged temperature profile. Although the time and horizontally volumeaveraged temperature profile is allowed to evolve during our simulations, it does not change substantially; this simplifies the calculation of the temperature perturbation. A layer defined by convective penetration is associated with negative vertical heat flux. Figure 4 also shows the timeaveraged and horizontally volumeaveraged vertical heat flux for simulation YS4. The penetration layer is clearly identifiable as the negative peak of vertical heat flux surrounding the bottom of the convection zone. There is significant discrepancy between the penetration depth defined from the averaged vertical kinetic energy flux and the averaged vertical heat flux. Because of this, we examine both measures in detail. We note that the use of Lagrangian coherent structures (e.g. Hadjighasem et al. 2016) to define the shape of convective plumes offers an alternative way to understand this discrepancy.
Fig. 5 Angular structure of the penetration layer at three arbitrary times spread throughout simulation YS5. The penetration depth in this illustration is determined by zeros of the vertical kinetic energy flux. The boundary between the convection zone and the stable radiative zone is indicated by a solid black line. The vertical axis is in units of the pressure scale height h_{p} at this boundary. A dashed black line indicates the average penetration depth at this time. 
Fig. 6 PDF of the penetration position, r_{o}, based on plumes penetrating at all angles and sampled at a fixed time interval throughout the full time span of each simulation indicated in Table 1. a) Penetration depth calculated from the vertical kinetic energy flux for the local convection simulations YS02. b) Penetration depth calculated from the vertical kinetic energy flux in the nonlocal convection simulations YS35. c) Penetration depth calculated from the vertical heat flux in the local convection simulations YS02. d) Penetration depth calculated from the vertical heat flux in the nonlocal convection simulations YS35. The heavy vertical black line marks the convective boundary. 
Horizontal and time averages have been used historically to define the width of the penetration layer. However the full statistics of convective penetration have not been previously examined, nor has the use of a simple average to represent them been specifically justified. Straightforward averages can mask contributions from the intermittent convective plumes that penetrate deeper into the stable radiative zone than the average plume. The distribution of penetration lengths may not be symmetric, or may have other relevant features that are not captured by a mean. Just such intermittent plumes can have a potentially larger impact on the physics of the penetration layer. Indeed ChristensenDalsgaard et al. (2011) find that consideration of temporal and spatial inhomogeneity of overshooting could lead to a thermal profile that agrees with observations from helioseismology. We therefore approach the full statistics of penetrating plumes with the goal of quantifying the contributions of intermittency.
3.2. Shape of the penetration layer
Using the first zero of the vertical kinetic energy flux and the vertical heat flux in the stable zone as our criteria, we examine the depth of convective penetration for simulations YS05, described in Table 1. In hydrodynamic simulations without rotation, the penetration is statistically independent of the colatitude θ. However at any instant in time the depth that plumes penetrate beyond the convection zone varies in θ, defining the structure of the penetration layer. Three typical snapshots that illustrate the shape of plumes penetrating into the radiative zone are shown in Fig. 5. In this figure, the bottom of the convection zone, established by the entropy profile and the Schwarzschild criterion, is shown as a thick horizontal line. The depth of convective penetration at any angular position in the spherical shell, based on the zero point of the vertical kinetic energy flux, is shown as a shaded area below this line. At the angular resolution of simulation YS5, penetrating plumes often fill several grid spaces. A general penetration depth calculated by a straightforward average in angle is indicated by the dashed horizontal line in each case; indeed we find that in Fig. 5a, for example, an average will hide the physically important more deeply penetrating plumes at approximately 40° and approximately 150°. This figure is calculated from the vertical kinetic energy flux; a qualitatively similar picture is produced from the vertical heat flux. We see similar pictures for each of the simulations considered in this work, independent of the resolution or boundary conditions used. Comparing with threedimensional hydrodynamic simulations during steady convection, we also observe similar pictures, although statistics based on longtimes are not available for threedimensional stellar convection.
To build a statistical model of convective penetration, we first consider the extent that plumes penetrate into the radiative zone at each angle defined in our simulation grid, and we sample the simulation data at a fixed time interval identical for each simulation. This fixed time interval is on the order of τ_{conv}/ 10^{3}, a value selected to capture the fastest moving plumes entering the penetration layer. We define the penetration position r_{o} as the position where convective motions cease based on either of our criteria. Figure 6 shows the probability density function (PDF) of this penetration position for each set of simulations. Each PDF is calculated using bins of equal size, which are identical for each set of simulations, and is properly normalized. From these PDFs, two important observations can be made. First, the PDFs possess either a multimodal shape or a significant shoulder. This identifies two distinct layers within the region below the bottom of the convection zone. Immediately below the convection zone is a shallow layer where weak plumes frequently penetrate the radiative zone. Below this layer is a second layer where stronger plumes penetrate, and more dramatically affect the stable radiative zone. This strong penetration occurs with a typical probability between approximately 1 / 2 and 1 / 7th of the shallow layer. The second important observation arises from the probability of larger penetration: these PDFs have a heavy tail that is nonGaussian in appearance. Both of these observations suggest that a straightforward average does not accurately describe the data.
Unlike the averaged profiles, the PDFs obtained based on the zeroes of the vertical kinetic energy flux and the zeroes of the vertical heat flux in Fig. 6 yield characteristically similar results. The PDFs calculated from the vertical heat flux are slightly more sharply peaked for small penetration extent than those calculated from the vertical kinetic energy flux. However, the peaks and shoulders are present and their approximate positions are unchanged regardless of the physical criterium used.
Meaningful differences between the PDFs calculated for the local simulations YS02 and the nonlocal simulations YS35 are clear from Fig. 6. The PDF of the nonlocal simulations show that penetration in simulations YS35 reaches deeper into the radiative zone. In addition, the multimodel behavior of the PDFs and the heavy tail is more significant. This may be directly related to the higher velocities that arise for identicallyresolved simulations when a large portion of the star is simulated, as reported by Pratt et al. (2016). Additional effects, such as the larger range of spatial scales of plumes that appear in the nonlocal simulations, may also be responsible, but are more difficult to quantify.
4. Extreme value statistics of maximal convective penetration
4.1. Formulation
The heavy tails evident in Fig. 6 motivate a more detailed statistical examination to characterize the depth of convective penetration that is meaningful for modeling stellar evolution. To pursue this, we define a maximal penetration depth to be the lowest position in the radiative zone that is reached at a given time: (8)Here we select the minimum over data at different angles, at the same time. In Fig. 7 the PDF of this maximal penetration depth is compared with the PDF of the penetration depth of all plumes (as shown in Fig. 6) and the time and horizontally volumeaveraged Péclet number. The Péclet number reaches a minimum near the bottom of the convection zone, above r/R = 0.4. Between 0.35 <r/R< 0.4, however, the Péclet number raises again; this is the signature of waves excited by convective penetration, that form in a narrow band in radius below the convection zone in the young Sun. The rms velocity profile and the rms radial velocity profile are impacted by the velocity of the waves in this radial range. The PDF of the maximal penetration depth correlates with the position where the Péclet number has a maximum due to these waves. The quantity r_{max} thus pinpoints an important physical consequence of convective penetration, while the PDF of the penetration depth for all plumes r_{o} does not.
Fig. 7 PDF of the convective penetration depth for all plumes r_{o} (dashed line), PDF of the maximal penetration depth at any time r_{max} (solid line), and the time and horizontally volumeaveraged Péclet number (dotted line) for simulation YS4. Each of these quantities is vertically shifted and scaled for comparison. The convective boundary is delineated by a heavy vertical line. 
Fig. 8 Cumulative distribution function F of the maximal penetration length, Δr_{max} defined in Eq. (9). a) Maximal penetration length calculated from the vertical kinetic energy flux for the local convection simulations YS02. b) Maximal penetration length calculated from the vertical kinetic energy flux for the nonlocal convection simulations YS35. c) Maximal penetration length calculated from the vertical heat flux for the local convection simulations YS02. d) Maximal penetration length calculated from the vertical heat flux for the nonlocal convection simulations YS35. Data is shown as point symbols, while a fit to the GEVD form is shown as a line. 
Because the PDF of the maximal penetration depth targets the waves that are the most obvious dynamical effect of penetration, it is reasonable to suggest that it should determine the relevant penetration length (or overshooting length) for stellar evolution calculations. We therefore define a maximal penetration length to be the difference between the maximal penetration depth and the convective boundary: (9)Here r_{B} is the radial position of the convective boundary and the maximum is taken over all angles θ at a single time t. Although in this work we treat penetration below a convection zone, the absolute value allows us to refer to penetration above or below a convective boundary. This construction allows us to focus on the larger amount of convective penetration that is likely to influence stellar structure.
The maximal penetration length Δr_{max} could be used straightforwardly as an overshooting length in simulations of stellar evolution. A simple average of Eq. (9)produces an overshooting length ℓ_{ov} ~ 0.3h_{p} in our simulations, where h_{p} is the pressure scale height at the convective boundary. A calculation of an overshooting length based on a time and horizontally averaged profile of the vertical kinetic energy flux, as in Fig. 4, produces a much smaller overshooting length, ℓ_{ov} ~ 0.1h_{p} from our simulations. In recent years stellar calculations have adopted overshooting lengths in a wide range 0.05h_{p} ≤ ℓ_{ov} ≤ 0.4h_{p} (Basu 1997; Chen & Han 2002; Brummell et al. 2002; Rogers et al. 2006; Politano et al. 2010; Liu et al. 2012; Montalbán et al. 2013; Jin et al. 2015), which includes both of these values. Both of these values are consistent with previously reported results for simulations with a similar stiffness to the young Sun (Brummell et al. 2002; Rogers et al. 2006). Those earlier results demonstrate that the overshooting length is dependent on the Prandtl number, the Rayleigh number, and the resolution of the simulation, as well as the stiffness. However our analysis suggests that larger values of the overshooting length than previously predicted would better reproduce the physics of large Péclet number convective penetration in the interior of a young lowmass star.
4.2. Application of extreme value theory
The maximal penetration length in Eq. (9)suggests concepts from extreme value theory (Castillo et al. 2005). In extreme value theory, unlike the central limit theorem that underlies much of probability theory, the interest is not the central values of an underlying distribution, but in accurately modeling the tails. The centerpiece of extreme value theory is the generalized extreme value distribution (GEVD), which can be used to model the probability of maximal events. The cumulative distribution function (CDF) identified with the GEVD has the form: (10)Here κ is commonly called the shape parameter, μ is the location parameter, and λ is the scale parameter (see e.g. Sect. 9.1.1 Eqs. (9.3) and (9.4) of Castillo et al. 2005; CharrasGarrido & Lezaud 2013; Gomes & Guillou 2015). Historically the welldocumented Frechét distribution corresponds to a positive shape parameter κ> 0, the Weibull distribution corresponds to κ< 0, and the Gumbel distribution corresponds to κ = 0. The shape parameter is related directly to the heaviness of the tail of the underlying distribution. For this reason κ is also sometimes called the extreme value index (EVI). For the Gumbel case of κ = 0 the GEVD may be expressed in the simpler form: (11)Thus the graph of the natural log of the negative natural log of F will be a straight line if the distribution is of the Gumbel form. It will curve downward if the distribution is of the Weibull form, and upward if the distribution is of the Frechét form. When a strict upper limit exists for the variable considered, the Weibull case is predicted. Our maximal penetration length is defined so that 0 < Δr_{max}< 0.43R, meaning that convective plumes cannot penetrate beyond the center of the star. Thus a Weibull distribution is expected.
The CDFs of the maximum penetration length for our simulations are shown in Fig. 8. The forms of the CDFs in Fig. 8 indeed appear to have slight downward curvature, consistent with the Weibull form. For all of these CDFs, most data lie at lower values of Δr_{max}, and thus the distribution is clearest in this region; a lack of data at high values of Δr_{max}, due to the comparative rarity of highly penetrating plumes, may contribute to the downward curvature in this region. However simulations with different resolution produce remarkably similar CDFs, although different resolutions produce different convective velocities, and different spans of time are examined. Increased resolution of a simulation leads to an improved characterization of the physical flows. In a direct sense, it also leads to a larger sample size over which our maximum is taken to define the length Δr_{max}. However in a fluid simulation the data sampled is not strictly independent, so that it is not clear to what extent a higher resolution is related to an improvement in the use of extreme value theory, where taking a maximum over a large sample of independent data is assumed. We do find that identical distributions are produced when our data is sampled ten times less frequently in time, on the order of the autocorrelation time of the flow at the bottom of the convection zone.
To precisely determine the parameters of the GEVD for each of our simulations, we use the package evd (Stephenson 2002; Penalva et al. 2013) publicly available for R (The R Project for Statistical Computing)^{3}. Independently, we confirm these parameters using a nonlinear least squares fit to the GEVD. Although in our fluid simulations the data is not strictly independent, as assumed by extreme value theory, we find that the GEVD provides an excellent fit. The shape, location, and scale parameters resulting from our fit are summarized in Table 2. In the table a value of κ = 0 is indicated when the best fit is of the Gumbel form. The value of the shape parameter κ is either small and negative, or zero for all of our simulations. The available timeseries of data for our highest resolution simulations is shorter than for lower resolutions. Despite this, a statistically significant trend is observed in which the value of κ is explicitly zero in the high resolution, high velocity simulations YS4 and YS5. Based on this observation, we predict that the Weibull characterization of our distributions should converge toward the simpler Gumbel distribution at high resolution, realistic stellar velocities, and when stellar convection is observed over long times.
Parameters for the generalized extreme value distribution of maximal convective penetration length Δr_{max}.
4.3. Radiation with the local surface temperature
To understand whether the treatment of the surface produces different results for convective penetration, we perform an additional simulation YS6. This simulation has identical resolution to simulations YS0 (local convection) and YS3 (nonlocal convection), but in simulation YS6 the nearsurface layers are included and the surface of the star is allowed to radiate energy with the local surface temperature. The PDF of the penetration depth of all plumes, r_{o}, for YS6 is highly similar to simulation YS3; this comparison is shown in Fig. 9. The penetration depth calculated from the vertical kinetic energy flux in Fig. 9a reveals a slightly larger probability of penetration between 0.39 <r< 0.42 than in simulation YS3. The penetration depth calculated from the vertical heat flux in Fig. 9b reveals a slightly larger probability of penetration between 0.42 <r< 0.43 than in simulation YS3. Although each of these simulations is observed over more than a hundred convective turnover times, the small difference in the statistics may still result from the limited length of time that these simulations are observed, particularly because convective penetration is intermittent.
Fig. 9 PDF of the penetration extent of all plumes, r_{o}. a) Penetration depth calculated from the vertical kinetic energy flux. b) Penetration depth calculated from the vertical heat flux. The convective boundary is delineated by a heavy vertical line. 
The parameters for the GEVD of the maximal penetration length in simulation YS6 are included in Table 2. For both the vertical kinetic energy flux and the vertical heat flux criteria, the GEVD for simulation YS6 is best fit by a Gumbel distribution. The treatment of the surface in simulation YS6 is physically more accurate than in YS3. That these two different treatments of the surface radiation produce similar results may be due to the fact that both treat the surface only approximately and neglect important physical effects. However the variable surface radiation in YS6 leads to higher velocity throughout the convection zone compared with either simulation YS3 or YS0; this was observed in Pratt et al. (2016). This is a significant observation, because it reinforces the finding that the overshooting length is not solely linked to the velocity of structures entering the penetration layer (Schmitt et al. 1984; Hurlburt et al. 1994; Saikia et al. 2000; Zahn 2000).
5. Determination of a diffusion coefficient for convection
5.1. Background
Convective overshooting and penetration have long been modeled by analyzing the enhanced diffusion in the penetration layer due to intermittent convective flows (e.g. as described by Freytag et al. 2010; Noels et al. 2010; Zhang 2013). Such models are based on the definition of a onedimensional diffusion coefficient, as in (12)Here A is any scalar quantity, and D is a diffusion coefficient. Based on dimensional arguments, a diffusion coefficient has typically been estimated for astrophysical applications by either a characteristic length scale multiplied by a characteristic velocity scale (e.g. Van Ballegooijen 1982; Andrássy & Spruit 2015), or a characteristic squared velocity multiplied by a characteristic time scale.
Equation (12)typically describes processes like molecular diffusion, and does not explicitly account for additional mixing due to largescale convective motions, which enhance diffusion. In the setting of steady twodimensional largescale convection^{4}, the relation between the diffusion coefficient and the Péclet number of a convective flow has been studied in both the small and large Péclet number regimes (Moffatt 1983; Rosenbluth et al. 1987; Shraiman 1987; Haynes & Vanneste 2014). In the small Péclet number regime this relation has been analytically and numerically shown to be: (13)Here D is the enhanced diffusion coefficient and D_{0} is the diffusion coefficient based on smallscale diffusive processes, which may be molecular diffusion in a liquid, or smallscale turbulent diffusion in a fluid. The constant c_{1} is related to the aspect ratio of the convection rolls. In the large Péclet number regime the corresponding relation is: (14)where c_{2} is a constant based on the aspect ratio of the convection cells. Rosenbluth et al. (1987) evaluate c_{2} for a wide range of possible convection roll aspect ratios, and find that with only small variation. The enhanced diffusion coefficients in Eqs. (13)and (14)are derived for an ideal setting rather than stellar convection, but are suggestive of diffusive scalings relevant to different Péclet number regimes.
Near the surface of a star, the small Péclet number limit is most appropriate; the Péclet number in simulations of this layer has been reported to be of order one. In this setting, Eq. (13)indicates that an enhancement of the diffusion coefficient dependent on the squared Péclet number is reasonable. If the thermal diffusivity and characteristic length scale are assumed to be approximately constant, the Péclet number is proportional to the velocity. This situation can be related directly to the enhanced diffusivity proposed by Freytag et al. (1996). In their radiation hydrodynamics simulations, Freytag et al. (1996) observe exponential decay of the rms radial velocity beyond the boundary of the convection zone. They use this result to define a piecewisecontinuous diffusion coefficient that is constant in the convection zone. Beyond the convection boundary, the piecewisecontinuous diffusion coefficient is defined using the Ansatz that it should decay exponentially as the velocity: (15)Here v_{B} and r_{B} are the rms radial velocity at the convective boundary, and the position of the convective boundary, respectively. The pressure scaleheight h_{p} at the convective boundary is used to define the characteristic time scale τ = h_{p}/v_{B} and the velocity scaleheight at the convective boundary h_{v} = fh_{p}, where f is a constant of proportionality. The exponentially decaying diffusion coefficient of Eq. (15)has been incorporated in several stellar evolution codes, including MESA (Paxton et al. 2013), ATON 3.1 (Ventura et al. 2008), and GARSTEC (Weiss & Schlattl 2008).
5.2. Diffusion coefficient in the large Péclet number regime
Fig. 10 Timeaveraged and horizontally volumeaveraged rms radial velocity near the bottom of the convection zone in simulation YS4. The timeaveraged and horizontally volumeaveraged squareroot of the Péclet number is shown, scaled and shifted, for comparison. The heavy vertical black line marks the convective boundary. 
In contrast with the early simulations examined by Freytag et al. (1996), in our comparatively highresolution young Sun simulations the lower convective boundary is in the large Péclet number regime. For this physically and numerically distinct case, we observe that the rms radial velocity does not decay exponentially beyond the convective boundary. The decay of the natural log of rms radial velocity is shown in Fig. 10, and it does not decay linearly with radius in the penetration layer. Motivated by Eq. (14), we examine the squareroot of the Péclet number in this figure in addition to the rms radial velocity, and find that predictably it exhibits a trend of decay outside the convection zone similar to the rms radial velocity. The observed departure from exponential decay motivates our proposal of a new form for the diffusion coefficient for stellar interiors in the large Péclet number regime.
Fig. 11 Time and horizontally volumeaveraged squareroot of the Péclet number compared with the data for the cumulative distribution function F(r_{max}/R) = 1 − F(Δr_{max}/R) calculated from the vertical kinetic energy flux. The best fit to this data from simulation YS4 is a Gumbel distribution with the parameters indicated in Table 2. The Péclet number has been normalized and shifted vertically. The dashed line indicates a Péclet number that is also shifted in radius so that the common shape is clear. The heavy vertical black line marks the convective boundary. 
Based on our statistical analysis in Sect. 4, we predict that the maximal penetration length should match the profile of the Gumbel distribution: exp( − exp( − x)). Figure 11 shows that the radial profile of the squareroot of the Péclet number closely matches the form of the CDF of maximal penetration length. We propose a diffusion coefficient of the form: (16)Here penetration beneath a convection zone is assumed, and Pe_{B} is a characteristic Péclet number in the convection zone. In constructing this diffusion coefficient, we have used the identity: (17)to relate the cumulative distribution functions of the maximal penetration depth r_{max} and the maximal penetration length Δr_{max}. This identity allows us to express the diffusion coefficient in the same form as Freytag et al. (1996), as a difference between the radial coordinate and the position of the convective boundary indicated by r_{B}. The prefactor D_{0}Pe_{B}^{1 / 2} in Eq. (16)is representative of diffusion in the convection zone, and may be estimated by comparison with stellar evolution models for convective diffusion. Using mixing length theory, diffusion in the convection zone is characterized by D_{MLT} = 1 / 3 L_{MLT}v_{MLT}, where the length scale is proportional to the pressure scale height^{5}. The values for μ and λ in Eq. (16)may be estimated from simulations YS35 in Table 2. In units of the pressure scale height at the base of the convection zone, h_{p} ≈ 0.178 R, these parameters give a range of We note that μ is approximately equivalent to the maximal penetration length estimated in Sect. 4. As with D_{F96}, the diffusion coefficient in Eq. (16)should be applied only in the layer of the star affected by convective penetration, so that it does not introduce infinitesimally small amounts of mixing into deeper layers of the star. This proposed diffusion coefficient is shown in Fig. 12, calibrated using the mixinglengththeory diffusion coefficient, D_{MLT}, used to produce the young Sun model. Both D_{MLT}, and D_{F96} are shown for comparison. The mixinglengththeory diffusion coefficient does not include any diffusion in the penetration layer. A detailed comparison of the exponentially decaying diffusion coefficient D_{F96} with our diffusion coefficient in Eq. (16)naturally depends on the parameters that define the exponential decay. These are commonly defined by treating the constant of proportionality between the velocity scaleheight and pressure scaleheight as a free parameter f ≤ 0.1 (e.g. Angelou et al. 2015). For these small values of f our diffusion coefficient models a higher level of diffusion throughout the penetration layer than the exponentially decaying diffusion coefficient of D_{F96}. Because of the higher level of diffusion predicted, and the different form of its decay below the convective boundary, it is worth exploring the effect of this new diffusion coefficient on the abundance of light elements like lithium and beryllium in solartype stars (Boesgaard 1976; Pinsonneault 1994), and to characterize its observational signatures. In simulations that include magnetic fields, this differing amount of diffusion due to convective penetration may affect magnetic flux in the penetration layer, and thus have consequences for the stellar dynamo (e.g. Nordlund et al. 1992; Rüdiger & Brandenburg 1995; Tobias et al. 2001; Rempel 2003; Brummell et al. 2010).
Fig. 12 Diffusion coefficient based on EVT proposed in Eq. (16)compared with a diffusion coefficient defined from mixing length theory D_{MLT}, and the decaying exponential diffusion coefficient D_{F96} of Freytag et al. (1996) defined in Eq. (15). Here the peak of the mixing length theory diffusion coefficient has been used to calibrate diffusion in the convection zone. 
6. Summary and discussion
We have studied flows that penetrate below a convection zone due to largescale twodimensional stellar convection in a prototypical young lowmass star. We have analyzed and compared convective penetration in two fundamental types of simulations: (1) those that model local convection by truncating the spherical shell to a minimal radial layer around the convective boundary, and (2) those that model nonlocal convection by simulating nearly the whole star. For each of these situations we examine stellar convection at three different resolutions, allowing us to observe penetration relating to a range of characteristic convective velocities. To compare with our simulations that hold the energy flux constant at the outer boundary, we also consider a nonlocal convection simulation where the surface radiation is allowed to vary with the surface temperature. This simulation is characterized by higher velocities in the convection zone. Our simulations do not account for rotation, shear flows, chemical mixing, or magnetic fields. These processes may change the properties of the flow at the convective boundary and alter the quantitative results. Analysis of these effects will be explored in future work.
Historically, a time and horizontally volumeaveraged radial profile of the vertical kinetic energy flux has often been used to determine the depth of the penetration layer. Different results are produced from a comparison of the time and horizontally volumeaveraged radial profiles of the vertical kinetic energy flux with the vertical heat flux, a second way to define the penetration layer. Here we show that when plumes at all angles and sampled at a fixed time interval are considered, PDFs of convective penetration depth calculated from these two different fluxes are characteristically similar. The PDFs are nonGaussian, with a heavy tail. Based on these PDFs, we define two distinct layers that form between the convection zone and the stable radiative zone: a shallow layer where convective plumes frequently penetrate, and a deeper layer where convective plumes penetrate intermittently.
Motivated by the nonGaussian character of convective penetration, we examine the maximal penetration depth at any time. The PDF of this quantity is centered closely around the position where waves are excited by convective penetration in the radiative zone. This suggests that the maximal penetration depth is a physically significant length scale to define the penetration layer. Based on a statistical analysis of our simulation data, we determine that a reasonable estimate of the overshooting length in the young Sun is ℓ_{ov} ~ 0.3h_{p}. We observe that the cumulative distribution function of the maximal penetration length can be modeled accurately by a Weibull distribution with a small shape parameter. As both higher resolution and higher velocities are examined, this Weibull distribution appears to converge toward a Gumbel distribution. This approach toward measuring convective penetration is new. Further work is required to demonstrate that it accurately describes convective penetration when rotation, shear flows or magnetic fields are considered. However, it provides a promising avenue to quantify the mixing due to overshooting or penetration in stars. Analysis of the properties of the waves generated by convective penetration into the radiative zone will be explored in future work.
Building on these statistical results, as well as scalings for twodimensional convection, we propose a new form for the convectionenhanced diffusion coefficient suitable for onedimensional stellar evolution calculations. Unlike the D_{F96} diffusion coefficient (Freytag et al. 1996), our proposed diffusion coefficient is targeted for the large Péclet number flows characteristic of stellar interiors. This represents a step forward for stellar evolution modeling, because the D_{F96} diffusion coefficient has been broadly applied to both small and large Péclet number flows, outside of its range of demonstrated physical validity. Followup studies to explore the effect of our diffusion coefficient in onedimensional stellar evolution calculations are now required to analyze possible differences that may result (Baraffe et al., in prep.). Isolating observable signatures to discern the structure of overshooting and penetration layers is currently a challenge, but may also shed light on the differences between these two forms of the diffusion coefficient.
In this work we have focused on identifying the amount of mixing due to convective penetration, which can be modeled by a simple diffusion coefficient. However convective penetration also has a significant impact on the thermal profile of a penetration layer. A complete model of penetration therefore should include a treatment of heat transport relevant for stellar evolution calculations. Such a model is currently being pursued based on the simulations presented.
The local interplay of velocity, density, and temperature have been related to buoyancy breaking of plumes (e.g., Latour et al. 1981; Hughes & Proctor 1988; Spruit et al. 1990).
The R Project for Statistical Computing: https://cran.rproject.org/
Acknowledgments
The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework (FP7/20072013)/ERC grant agreement No. 320478.
I.B. thanks MOCA at Monash University for warm hospitality.
This work used the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment is funded by BIS National EInfrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. DiRAC is part of the National EInfrastructure.
This work also used the University of Exeter Supercomputers, ZEN and ISCA, funded by the STFC/DiRAC, the Large Facilities Capital Fund of BIS and the University of Exeter.
References
 Abbett, W. P., Beaver, M., Davids, B., et al. 1997, ApJ, 480, 395 [NASA ADS] [CrossRef] [Google Scholar]
 Alvan, L., Brun, A. S., & Mathis, S. 2014, A&A, 565, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Alvan, L., Strugarek, A., Brun, A., Mathis, S., & Garcia, R. 2015, A&A, 581, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Andersen, B. 1994, Sol. Phys., 152, 241 [NASA ADS] [CrossRef] [Google Scholar]
 Andrássy, R., & Spruit, H. 2015, A&A, 578, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Angelou, G. C., D’Orazi, V., Constantino, T. N., et al. 2015, MNRAS, 450, 2423 [NASA ADS] [CrossRef] [Google Scholar]
 Arnett, D. W. 2014, IAU Proc., 9, 459 [CrossRef] [Google Scholar]
 Arnett, W. D., & Meakin, C. 2009, IAU Proc., 5, 106 [Google Scholar]
 Baraffe, I., & El Eid, M. F. 1991, A&A, 245, 548 [NASA ADS] [Google Scholar]
 Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. 1997, A&A, 327, 1054 [NASA ADS] [Google Scholar]
 Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. 1998, A&A, 337, 403 [NASA ADS] [Google Scholar]
 Basu, S. 1997, MNRAS, 288, 572 [NASA ADS] [Google Scholar]
 Boesgaard, A. 1976, PASP, 88, 353 [NASA ADS] [CrossRef] [Google Scholar]
 BöhmVitense, E. 1958, ZAp, 46, 108 [NASA ADS] [Google Scholar]
 Brandenburg, A. 2016, ApJ, 832, 6 [Google Scholar]
 Browning, M. K., Brun, A. S., & Toomre, J. 2004, ApJ, 601, 512 [NASA ADS] [CrossRef] [Google Scholar]
 Brummell, N. H., Clune, T. L., & Toomre, J. 2002, ApJ, 570, 825 [Google Scholar]
 Brummell, N. H., Tobias, S. M., & Cattaneo, F. 2010, Geophys. Astrophys. Fluid Dyn., 104, 565 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., Alvan, L., Strugarek, A., Mathis, S., & Garcia, R. A. 2013, JPCS, 440, 012043 [Google Scholar]
 Brun, A. S., Strugarek, A., Varela, J., et al. 2017, ApJ, 836, 192 [NASA ADS] [CrossRef] [Google Scholar]
 Castillo, E., Hadi, A. S., Balakrishnan, N., & Sarabia, J.M. 2005, Extreme value and related models with applications in engineering and science (NJ: Wiley Hoboken) [Google Scholar]
 Chan, K. L., Cai, T., & Singh, H. P. 2010, IAU Proc., 6, 317 [CrossRef] [Google Scholar]
 CharrasGarrido, M., & Lezaud, P. 2013, Journal de la Société Française de Statistique, 154 [Google Scholar]
 Chen, X., & Han, Z. 2002, MNRAS, 335, 948 [NASA ADS] [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J., Monteiro, M., Rempel, M., & Thompson, M. J. 2011, MNRAS, 414, 1158 [NASA ADS] [CrossRef] [Google Scholar]
 Dintrans, B., Brandenburg, A., Nordlund, Å., & Stein, R. F. 2003, in Asteroseismology Across the HR Diagram (Springer), 237 [Google Scholar]
 Domaradzki, J. A., Xiao, Z., & Smolarkiewicz, P. K. 2003, Phys. Fluids, 15, 3890 [NASA ADS] [CrossRef] [Google Scholar]
 Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585 [NASA ADS] [CrossRef] [Google Scholar]
 Freytag, B., Ludwig, H.G., & Steffen, M. 1996, A&A, 313, 497 [Google Scholar]
 Freytag, B., Allard, F., Ludwig, H.G., Homeier, D., & Steffen, M. 2010, A&A, 513, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goffrey, T., Pratt, J., Viallet, M., et al. 2017, A&A, 600, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gomes, M. I., & Guillou, A. 2015, Int. Stat. Rev., 83, 263 [CrossRef] [Google Scholar]
 Gough, D., & Tayler, R. 1966, MNRAS, 133, 85 [NASA ADS] [CrossRef] [Google Scholar]
 GrimmStrele, H., Kupka, F., LöwBaselli, B., et al. 2015, New Astron., 34, 278 [NASA ADS] [CrossRef] [Google Scholar]
 Hadjighasem, A., Karrasch, D., Teramoto, H., & Haller, G. 2016, Phys. Rev. E, 93, 063107 [NASA ADS] [CrossRef] [Google Scholar]
 Haynes, P., & Vanneste, J. 2014, J. Fluid Mech., 745, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Hughes, D., & Proctor, M. 1988, Annu. Rev. Fluid Mech., 20, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Hurlburt, N. E., Toomre, J., & Massaguer, J. M. 1986, ApJ, 311, 563 [NASA ADS] [CrossRef] [Google Scholar]
 Hurlburt, N. E., Toomre, J., Massaguer, J. M., & Zahn, J.P. 1994, ApJ, 421, 245 [NASA ADS] [CrossRef] [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Jin, J., Zhu, C., & Lü, G. 2015, PASJ, 67, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Jonker, H. J., van Reeuwijk, M., Sullivan, P. P., & Patton, E. G. 2013, J. Fluid Mech., 732, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Knoll, D. A., & Keyes, D. E. 2004, J. Comput. Phys., 193, 357 [NASA ADS] [CrossRef] [Google Scholar]
 Latour, J., Toomre, J., & Zahn, J.P. 1981, ApJ, 248, 1081 [Google Scholar]
 Lecoanet, D., & Quataert, E. 2013, MNRAS, 430, 2363 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, Z., Pakmor, R., Röpke, F., et al. 2012, A&A, 548, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Meakin, C. A., & Arnett, D. 2007, ApJ, 667, 448 [NASA ADS] [CrossRef] [Google Scholar]
 Meynet, G., & Maeder, A. 2000, A&A, 361, 101 [NASA ADS] [Google Scholar]
 Moffatt, H. 1983, Rep. Prog. Phys., 46, 621 [NASA ADS] [CrossRef] [Google Scholar]
 Montalbán, J., Miglio, A., Noels, A., et al. 2013, ApJ, 766, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Muthsam, H., Göb, W., Kupka, F., Liebich, W., & Zöchling, J. 1995, A&A, 293, 127 [NASA ADS] [Google Scholar]
 Noels, A., Montalban, J., Miglio, A., Godart, M., & Ventura, P. 2010, Astrophys. Space Sci., 328, 227 [NASA ADS] [CrossRef] [Google Scholar]
 Nordlund, A., Brandenburg, A., Jennings, R. L., et al. 1992, ApJ, 392, 647 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Penalva, H., Neves, M., & Nunes, S. 2013, Metodoloski Zvezki, 10, 17 [Google Scholar]
 Pinçon, C., Belkacem, K., & Goupil, M. 2016, A&A, 588, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pinsonneault, M. 1994, in Cool Stars, Stellar Systems, and the Sun, ASP Conf. Proc., 64, 254 [NASA ADS] [Google Scholar]
 Politano, M., Van der Sluys, M., Taam, R. E., & Willems, B. 2010, ApJ, 720, 1752 [NASA ADS] [CrossRef] [Google Scholar]
 Pratt, J., Baraffe, I., Goffrey, T., et al. 2016, A&A, 593, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Radice, D., Couch, S. M., & Ott, C. D. 2015, Comput. Astrophys. Cosmol., 2, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Rempel, M. 2003, A&A, 397, 1097 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rempel, M. 2004, ApJ, 607, 1046 [NASA ADS] [CrossRef] [Google Scholar]
 Roe, P. 1986, Annu. Rev. Fluid Mech., 18, 337 [Google Scholar]
 Rogers, T. M., & Glatzmaier, G. A. 2005a, MNRAS, 364, 1135 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, T. M., & Glatzmaier, G. A. 2005b, ApJ, 620, 432 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, T. M., Glatzmaier, G. A., & Jones, C. 2006, ApJ, 653, 765 [NASA ADS] [CrossRef] [Google Scholar]
 Rosenbluth, M., Berk, H., Doxas, I., & Horton, W. 1987, Phys. Fluids, 30, 2636 [NASA ADS] [CrossRef] [Google Scholar]
 Rüdiger, G., & Brandenburg, A. 1995, A&A, 296, 557 [NASA ADS] [Google Scholar]
 Saikia, E., Singh, H. P., Chan, K., Roxburgh, I., & Srivastava, M. 2000, ApJ, 529, 402 [NASA ADS] [CrossRef] [Google Scholar]
 Schmalzl, J., Breuer, M., & Hansen, U. 2004, EPL, 67, 390 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schmitt, J., Rosner, R., & Bohn, H. 1984, ApJ, 282, 316 [NASA ADS] [CrossRef] [Google Scholar]
 Schranner, F. S., Domaradzki, J. A., Hickel, S., & Adams, N. A. 2015, Comput. Fluids, 114, 84 [CrossRef] [Google Scholar]
 Shiode, J. H., Quataert, E., Cantiello, M., & Bildsten, L. 2013, MNRAS, 430, 1736 [NASA ADS] [CrossRef] [Google Scholar]
 Shraiman, B. I. 1987, Phys. Rev. A, 36, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Spruit, H. C., Nordlund, A., & Title, A. 1990, ARA&A, 28, 263 [NASA ADS] [CrossRef] [Google Scholar]
 Stephenson, A. G. 2002, R News, 2, [Google Scholar]
 Strugarek, A., Brun, A., & Zahn, J.P. 2011, A&A, 532, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tennekes, H., & Driedonks, A. 1981, BoundaryLayer Meteorology, 20, 515 [NASA ADS] [CrossRef] [Google Scholar]
 Tian, C.L., Deng, L.C., & Chan, K.L. 2009, MNRAS, 398, 1011 [NASA ADS] [CrossRef] [Google Scholar]
 Tobias, S. M., Brummell, N. H., Clune, T. L., & Toomre, J. 2001, ApJ, 549, 1183 [NASA ADS] [CrossRef] [Google Scholar]
 Trampedach, R. 2010, Ap&SS, 328, 213 [Google Scholar]
 Van Ballegooijen, A. 1982, A&A, 113, 99 [NASA ADS] [Google Scholar]
 van der Poel, E. P., Stevens, R. J., & Lohse, D. 2015, J. Fluid Mech., 736, 177 [Google Scholar]
 Van Leer, B. 1974, J. Comput. Phys., 14, 361 [NASA ADS] [CrossRef] [Google Scholar]
 Van Leer, B. 1977, J. Comput. Phys., 23, 276 [Google Scholar]
 Ventura, P., D’Antona, F., & Mazzitelli, I. 2008, Astrophys. Space Sci., 316, 93 [Google Scholar]
 Viallet, M., Goffrey, T., Baraffe, I., et al. 2016, A&A, 586, A153 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vitense, E. 1953, Z. Astrophys., 32, 135 [NASA ADS] [Google Scholar]
 Weiss, A., & Schlattl, H. 2008, Astrophys. Space Sci., 316, 99 [Google Scholar]
 Zahn, J.P. 1991, A&A, 252, 179 [NASA ADS] [Google Scholar]
 Zahn, J.P. 2000, Annals of the New York Academy of Sciences, 898, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, Q. 2013, ApJS, 205, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Zhou, Y., & Thornber, B. 2016, J. Fluids Eng., 138, 070905 [CrossRef] [Google Scholar]
 Zhou, Y., Grinstein, F. F., Wachtor, A. J., & Haines, B. M. 2014, Phys. Rev. E, 89, 013303 [NASA ADS] [CrossRef] [Google Scholar]
 Ziegler, U., & Rüdiger, G. 2003, A&A, 401, 433 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
Parameters for twodimensional compressible hydrodynamic simulations of the young Sun.
Parameters for the generalized extreme value distribution of maximal convective penetration length Δr_{max}.
All Figures
Fig. 1 Typical time snapshots of radial velocity in simulation YS2 (left), YS5 (middle), and YS6 (right). The zero point is colored black, while inward flows are gray and outward flows are red. 

In the text 
Fig. 2 Radial profile of the timeaveraged and horizontallyvolumeaveraged Péclet number, as defined by Eq. (4), from simulation YS4. The contribution to the thermal diffusivity by numerical truncations are not included in this profile. The heavy vertical black line marks the boundary between the stable radiative zone and the convection zone, determined from the radial profile of entropy and the Schwarzschild criterion produced by the onedimensional stellar model. 

In the text 
Fig. 3 Radial profile of the timeaveraged and horizontally volumeaveraged kinetic energy, , scaled to its maximum, from simulation YS3. The heavy vertical black line marks the convective boundary. 

In the text 
Fig. 4 Radial profile of the timeaveraged and horizontally volumeaveraged vertical kinetic energy flux F_{k} defined in Eq. (6), and vertical heat flux F_{c} defined in Eq. (7)for simulation YS4. Both of these fluxes have been normalized to their absolute maximum value. The heavy vertical black line marks the convective boundary. 

In the text 
Fig. 5 Angular structure of the penetration layer at three arbitrary times spread throughout simulation YS5. The penetration depth in this illustration is determined by zeros of the vertical kinetic energy flux. The boundary between the convection zone and the stable radiative zone is indicated by a solid black line. The vertical axis is in units of the pressure scale height h_{p} at this boundary. A dashed black line indicates the average penetration depth at this time. 

In the text 
Fig. 6 PDF of the penetration position, r_{o}, based on plumes penetrating at all angles and sampled at a fixed time interval throughout the full time span of each simulation indicated in Table 1. a) Penetration depth calculated from the vertical kinetic energy flux for the local convection simulations YS02. b) Penetration depth calculated from the vertical kinetic energy flux in the nonlocal convection simulations YS35. c) Penetration depth calculated from the vertical heat flux in the local convection simulations YS02. d) Penetration depth calculated from the vertical heat flux in the nonlocal convection simulations YS35. The heavy vertical black line marks the convective boundary. 

In the text 
Fig. 7 PDF of the convective penetration depth for all plumes r_{o} (dashed line), PDF of the maximal penetration depth at any time r_{max} (solid line), and the time and horizontally volumeaveraged Péclet number (dotted line) for simulation YS4. Each of these quantities is vertically shifted and scaled for comparison. The convective boundary is delineated by a heavy vertical line. 

In the text 
Fig. 8 Cumulative distribution function F of the maximal penetration length, Δr_{max} defined in Eq. (9). a) Maximal penetration length calculated from the vertical kinetic energy flux for the local convection simulations YS02. b) Maximal penetration length calculated from the vertical kinetic energy flux for the nonlocal convection simulations YS35. c) Maximal penetration length calculated from the vertical heat flux for the local convection simulations YS02. d) Maximal penetration length calculated from the vertical heat flux for the nonlocal convection simulations YS35. Data is shown as point symbols, while a fit to the GEVD form is shown as a line. 

In the text 
Fig. 9 PDF of the penetration extent of all plumes, r_{o}. a) Penetration depth calculated from the vertical kinetic energy flux. b) Penetration depth calculated from the vertical heat flux. The convective boundary is delineated by a heavy vertical line. 

In the text 
Fig. 10 Timeaveraged and horizontally volumeaveraged rms radial velocity near the bottom of the convection zone in simulation YS4. The timeaveraged and horizontally volumeaveraged squareroot of the Péclet number is shown, scaled and shifted, for comparison. The heavy vertical black line marks the convective boundary. 

In the text 
Fig. 11 Time and horizontally volumeaveraged squareroot of the Péclet number compared with the data for the cumulative distribution function F(r_{max}/R) = 1 − F(Δr_{max}/R) calculated from the vertical kinetic energy flux. The best fit to this data from simulation YS4 is a Gumbel distribution with the parameters indicated in Table 2. The Péclet number has been normalized and shifted vertically. The dashed line indicates a Péclet number that is also shifted in radius so that the common shape is clear. The heavy vertical black line marks the convective boundary. 

In the text 
Fig. 12 Diffusion coefficient based on EVT proposed in Eq. (16)compared with a diffusion coefficient defined from mixing length theory D_{MLT}, and the decaying exponential diffusion coefficient D_{F96} of Freytag et al. (1996) defined in Eq. (15). Here the peak of the mixing length theory diffusion coefficient has been used to calibrate diffusion in the convection zone. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.