Sphericalshell boundaries for twodimensional compressible convection in a 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} MaxPlanckInstitut für Astrophysik, Karl Schwarzschild Strasse 1, 85741 Garching, Germany
Received: 11 February 2016
Accepted: 19 June 2016
Context. Studies of stellar convection typically use a sphericalshell geometry. The radial extent of the shell and the boundary conditions applied are based on the model of the star investigated. We study the impact of different twodimensional spherical shells on compressible convection. Realistic profiles for density and temperature from an established onedimensional stellar evolution code are used to produce a model of a large stellar convection zone representative of a young lowmass star, like our sun at 10^{6} years of age.
Aims. We analyze how the radial extent of the spherical shell changes the convective dynamics that result in the deep interior of the young sun model, far from the surface. In the nearsurface layers, simple smallscale convection develops from the profiles of temperature and density. A central radiative zone below the convection zone provides a lower boundary on the convection zone. The inclusion of either of these physically distinct layers in the spherical shell can potentially affect the characteristics of deep convection.
Methods. We perform hydrodynamic implicit large eddy simulations of compressible convection using the MUltidimensional Stellar Implicit Code (MUSIC). Because MUSIC has been designed to use realistic stellar models produced from onedimensional stellar evolution calculations, MUSIC simulations are capable of seamlessly modeling a whole star. Simulations in twodimensional spherical shells that have different radial extents are performed over tens or even hundreds of convective turnover times, permitting the collection of wellconverged statistics.
Results. To measure the impact of the sphericalshell geometry and our treatment of boundaries, we evaluate basic statistics of the convective turnover time, the convective velocity, and the overshooting layer. These quantities are selected for their relevance to onedimensional stellar evolution calculations, so that our results are focused toward studies exploiting the socalled 321D link. We find that the inclusion in the spherical shell of the boundary between the radiative and convection zones decreases the amplitude of convective velocities in the convection zone. The inclusion of nearsurface layers in the spherical shell can increase the amplitude of convective velocities, although the radial structure of the velocity profile established by deep convection is unchanged. The impact of including the nearsurface layers depends on the speed and structure of smallscale convection in the nearsurface layers. Larger convective velocities in the convection zone result in a commensurate increase in the overshooting layer width and a decrease in the convective turnover time. These results provide support for nonlocal aspects of convection.
Key words: methods: numerical / convection / stars: interiors / stars: lowmass / stars: evolution
© ESO, 2016
1. Introduction
Because convection underlies the fundamental processes of heat transport, mixing, shear, and the stellar dynamo, it also influences stellar evolution over long times. Studies of stellar evolution typically use onedimensional calculations that evolve physical quantities dependent on the radial position interior to a star. These studies model the impact of convection on other physical quantities using one of the several variations of stellar mixing length theory that depend on the local temperature gradient (e.g. Vitense 1953; BöhmVitense 1958; Abbett et al. 1997; Trampedach 2010; Brandenburg 2015). Nevertheless, nonlocal processes can significantly influence stellar convection (e.g. Spruit 1997; Gough & Weiss 1976; Canuto & Dubovikov 1997). Thus nonlinear hydrodynamic convection simulations that include a whole star are a key step to improving our understanding of stellar convection and eventually to improving models of stellar evolution.
Wholestar simulations present considerable numerical challenges. For that reason, numerical studies of hydrodynamic and magnetohydrodynamic convection are often performed in a simulation volume that is a spherical shell containing a portion of a star (e.g. Gilman 1983; Cole et al. 2014; Yadav et al. 2013; Käpylä et al. 2011a,b; Miesch et al. 2000; Grote & Busse 2001; Quataert & Gruzinov 2000; Nelson et al. 2011; Simitev et al. 2011). A spherical shell refers simply to any portion of a spherical domain that is limited both in angular and radial extent. A spherical shell can be twodimensional or threedimensional, and sometimes is called a spherical wedge. Several earlier works (Cole et al. 2016; Mitra et al. 2009; Heimpel et al. 2005) have studied the impact of the sphericalshell geometry and boundaries on the dynamo. The goal of this work is to produce a similar study of the impact of the sphericalshell geometry on compressible convection without the fundamental influence of magnetic fields.
Although the impact of the sphericalshell geometry on compressible convection has not been studied, much fundamental work has been done to study how far convective plumes overshoot the boundary into the stable region. Earlier work on overshooting during compressible convection has concentrated primarily on direct numerical simulation in a local domain with Cartesian geometry in two (Roxburgh & Simmons 1993; Hurlburt et al. 1994) and three dimensions (Ziegler & Rüdiger 2003; Brummell et al. 2002, 2010). Typically, a convectively unstable layer is stacked on top of a stable layer. The structure of this model is formed using a piecewise continuous twolayer polytropic stratification in the density, temperature, and thermal diffusivity. Direct numerical simulations using the anelastic approximation rather than compressible convection have also been performed in a twodimensional cylindrical geometry (Rogers & Glatzmaier 2005; Rogers et al. 2006); these works used a realistic stellar stratification rather than a polytropic stratification, but the thermal diffusivity was artificially increased. Several other early studies have produced large eddy simulations, typically using a standard Smagorinsky subgrid scale model, in a local domain with Cartesian geometry (Xie & Toomre 1993; Singh et al. 1995, 1998; Saikia et al. 2000; Pal et al. 2007). Direct numerical simulation is a sensible choice when the focus is dynamo action; when a magnetic field is present, smallscale motions feed back on largescale motions in ways that are not completely understood (Miesch et al. 2015). However when hydrodynamic convection is the focus, the largest scales can potentially be the most important.
Because nonlocal effects of convection and coupling between physically different regions in a star may play a significant role in the basic character of convective processes (e.g. Brun et al. 2011; Latour et al. 1981; Spruit 1997), the necessary radial extent of the spherical shell is an open question. For practical reasons, simulations of stellar convection in a spherical shell often include a minimal radial extent of a modeled convection zone; neighboring, physically distinct zones of the star are neglected, and closed boundary conditions are applied at the top and bottom of this convectively unstable region of the star. In this work we compare results obtained when the interaction between a central radiative zone and an interior convection zone, and between an interior convection zone and nearsurface layers is allowed. We allow or disallow this coupling either by including these extra zones of the star in the spherical shell, or by excluding them and applying a closed boundary condition. Our motivation is to establish a physically reasonable position for the spherical shell and treatment of the boundaries that can serve as a basis for future studies of two and threedimensional convection at the bottom of the stellar convection zone, and overshooting of the boundary to the radiative zone.
In contrast to other solar and stellar simulations, the physical model of a star that we study in this work only allows for convection; the possibility of studying additional physical effects such as rotation, the effect of a tachocline, chemical mixing, or magnetic fields is intentionally omitted from the current study. This simplifies the boundaries between the lower radiative zone, the convection zone, and the nearsurface layers. We study a prototypical model of a young, lowmass star: the young sun. 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 of the star. Figure 1 sketches the radial structure of the young sun model. A central radiative zone below a large convection zone is expected based on the radial entropy profile and an evaluation of the Schwarzschild criterion. An overshooting layer, where mixing of physically different flows from the radiative and convection zones can take place, can only be defined dynamically; this is discussed in detail in Sect. 4.3. More importantly for this work, the young sun model has a large convection zone; it is convectively unstable over 1.2 × 10^{11} cm of the total radius of 2.13 × 10^{11} cm. This large convective envelope allows us to study deep stellar convection, far from the physically complicated nearsurface layers. We term the type of convection analyzed in this work “stellar convection” because both the stratification in the density and the temperature gradient that drives the convection are nonuniform and are not linearly dependent on the radius. We define the nearsurface layers as the portion of the star where the pressure scale height, h_{p} = −p/ (∂p/∂r), drops dramatically. Further details of the young sun model are discussed in Sect. 2.
The most thoroughly studied type of convection is RayleighBénard convection, which is driven by a uniform linear gradient of temperature, carefully controlled in a laboratory environment. In RayleighBénard convection, patterns can be identified (Bodenschatz et al. 2000; Weiß et al. 2012; Emran & Schumacher 2015). In a star, the gradients of temperature and density are generally nonuniform and vary nonlinearly over the full stellar radius. The radial variation of temperature and density for our young sun model is shown in Fig. 1. Although coherent structures such as plumes and convection rolls form during convection driven by nonuniform gradients, regular patterns cannot be identified. Instead the diagnostics we use to evaluate the impact of the sphericalshell geometry and boundaries are velocity statistics gathered over long periods of steady convection.
Fig. 1 Radial structure of our young sun model, showing the profiles of the natural log of the pressure scale height h_{p}, the temperature T, and density ρ. The pressure scale height and temperature are scaled by arbitrary values in order to appear on the scale of this diagram. R is the total radius of the young sun. Average radial widths shown for physical zones are predicted by the young sun model calculated with the Lyon stellar evolution code. 

Open with DEXTER 
This work is structured as follows. In Sect. 2 we discuss the simulation framework by outlining the physical and numerical models used by the MUltidimensional Stellar Implicit Code (MUSIC). In Sect. 3 we discuss the details of several spherical shells that we formulate. In Sect. 4 we compare simulations performed in different spherical shells with different radial extents on the basis of three different statistics related to the convective velocities. In Sect. 5 we discuss the implications of these results for our planned future work on deep stellar convection and convective overshooting.
2. Simulations
Convection is typically modeled using one of three approaches: the Boussinesq hydrodynamic equations, the anelastic hydrodynamic equations, or the compressible hydrodynamic equations (for a concise summary see Glatzmaier 2014). While the Boussinesq and anelastic approximations are simpler and numerically more efficient, there is reason to believe that these physical approximations could affect the realism of basic results in simulations of a whole star (Lecoanet et al. 2014; Verhoeven et al. 2015). In this work we use the compressible hydrodynamic equations.
The MUSIC code solves the inviscid compressible hydrodynamic equations for density ρ, momentum ρu, and internal energy ρe:
Here the thermal conductivity χ = 16σT^{3}/ 3κρ is defined using the StefanBoltzmann constant σ and the Rosseland mean opacity κ. The thermal conductivity is thus realistic to the young sun model, and has not been enhanced (for examples of ways that the thermal conductivity may be enhanced, see Browning et al. 2004; Rogers et al. 2006; Tian et al. 2009; Strugarek et al. 2011). A major feature of our simulations is the use of an equation of state and realistic opacities that are standard for onedimensional stellar evolution calculations. Opacities are interpolated from the OPAL (Iglesias & Rogers 1996) and Ferguson et al. (2005) tables, which cover a range of temperatures suitable for the description of the entire structure of a lowmass star. The compressible hydrodynamic Eqs. (1)−(3)are closed by determining the gas 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 from 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 spherically symmetric vector identical to that used in the Lyon stellar evolution code, and not evolved by our simulations. Thus MUSIC simulation results should contribute to the 321D link (Arnett 2014; Arnett & Meakin 2009), i.e. the effort to improve onedimensional stellar evolution models by studying critical short phases using stellar hydrodynamics in two and three dimensions.
Parameters for twodimensional hydrodynamic simulations of the young sun.
Early work to develop the MUSIC code was reported in detail in Viallet et al. (2011, 2013). Time integration is implicit in order to permit time steps larger than the CourantFriedrichsLewy (CFL) limit for timeexplicit methods. Convergence of an implicit scheme can be computationally demanding. In the present study, the system of equations is discretized in time using a CrankNicolson scheme. To integrate the compressible hydrodynamic equations in time, a Jacobian free NewtonKrylov (JFNK) solver (Knoll & Keyes 2004) is employed. Instead of storing a Jacobian, a JFNK method uses matrixvector products that can be estimated efficiently with finite difference methods. In MUSIC, a physicsbased preconditioner (Park et al. 2009) targets the physical processes that cause the system to become stiff. In practice, this preconditioning matrix is calculated as a semiimplicit solution of the system (Viallet et al. 2016). In all simulations in this work, we require a general CFL number to be less than 10, while a CFL number based on simple advection is restricted to 0.5. The time step calculated from these constraints produces good convergence of relevant basic quantities, such as average kinetic energy. The development of efficient implicit solvers is the subject of ongoing research (e.g. Chacón 2008; Main & Farhat 2014; Alvarez Laguna et al. 2016).
The spatial discretization of Eqs. (1)−(3)is accomplished using a staggered grid and a finite volume approach. Physical quantities are interpolated to the grid using an upwind limited interpolation similar to the Monotone Upwind Schemes for Conservation Laws (MUSCL) method (Thornber et al. 2008). This method employs the wellknown van Leer flux limiter (Roe 1986). MUSIC is designed as a large eddy simulation (LES). In the present work no additional viscosity is applied either through fixed coefficients or subgridscale modeling. This is a common tactic for astrophysical hydrodynamics, in order to obtain the minimum possible dissipation (for a similar discussion see Miesch et al. 2015). The MUSIC code has been benchmarked (Goffrey et al. 2016) against standard hydrodynamic problems that isolate fundamental physical mechanisms relevant to stellar hydrodynamics, including an ideal RayleighTaylor instability, an ideal KelvinHelmholtz instability, and the TaylorGreen vortex. For these standard hydrodynamic problems MUSIC produces the expected results, and exhibits the expected numerical convergence properties.
3. Sphericalshell simulation volumes
3.1. Resolution
The compressible hydrodynamic Eqs. (1)−(3)are solved in a twodimensional spherical shell using spherical coordinates: radius r and colatitude θ. In studying the impact of the spherical shell on compressible convection we examine only twodimensional convection. Twodimensional convection is known to produce higher velocity structures than threedimensional convection (see for example Meakin & Arnett 2007). Boundary effects are also larger in two dimensions than in three. Thus a study on the impact of boundary conditions and spherical shells on twodimensional convection is the extreme case, and our results will also be relevant for threedimensional convection simulations.
We consider the suite of twodimensional simulations of spherical shell convection summarized in Table 1. In this table, the radial position and extent of the spherical shell are shown in units of the total stellar radius R. Most of the simulations use a uniform grid characterized by a fixed radial spacing Δr. The lowresolution simulations (Low18) use radial grid spacing Δr/R ≈ 2.8 × 10^{3}. The highresolution simulations (Hi14) use radial grid spacing Δr/R ≈ 1.4 × 10^{3}. The extrahighresolution simulation ExH2 uses radial grid spacing Δr/R = 7 × 10^{4}. Each of these simulations that uses a uniform grid also uses a constant energy flux across the upper surface.
To explore the impact of certain boundary conditions on the energy, described in detail in Sect. 3.2, a higher resolution of the nearsurface layers is required. To accomplish this, the remaining two simulations (Low9 and Hi5) use a spliced grid. Our spliced grid is composed of a uniform grid with fixed grid spacing Δr for r/R ≤ 0.94. In the nearsurface layers, for r/R> 0.94, a nonuniform grid that is decreasing radially toward the surface is spliced on top of the uniform grid. The nonuniform portion of the grid is defined by a geometric sequence Δr_{i} = Δr_{i + 1}/ 1.05 (Geroux et al. 2016). The grid spacing in the uniform portion of the grid in simulations Low9, and Hi5 is identical to that in simulations Low18 and Hi14, respectively. The spliced grid of simulation Hi5 allows for a resolution of Δr/h_{p}(r) ≈ 0.36 at the surface; for comparison, the uniform grid in simulation Hi4, has Δr/h_{p}(r) ≈ 0.67 at the surface. For either the spliced grid or the uniform grid, our grid spacing in the nearsurface layers is too low to resolve the complex physics in these layers with any precision. Our motivation for modeling these layers is solely to provide a simple physically motivated open boundary condition (e.g. Käpylä et al. 2010; Cameron et al. 2011) on the interior convection zone. Such an open boundary condition allows the exchange of momentum, density, and thermal fluctuations with other zones of the star; open boundary conditions are highly desirable when dealing with largescale flows.
The resolutions we employ are comparable to earlier LES studies (e.g. Brun & Palacios 2009); we note that simulation Hi2 has a total grid size of 608 × 512. By not pursuing extremely highresolution simulations for this benchmarking study of largescale convection in a sphericalshell geometry, we are able to simulate over long times. Because the onedimensional stellar evolution model of the young sun stops at the photosphere, defined by an optical depth τ = 2 / 3, our simulations do not include the usual increase in entropy in the stable atmospheric layers (e.g. Magic, Z. 2016; Trampedach et al. 2014; Abbett et al. 1997). The preceding drop in entropy in the nearsurface layers is underresolved. However, the entropy jump at the bottom of the convection zone is resolved in our simulations. The entropy profile for simulations Hi4 and Hi5 is shown in Fig. 2. Because the characteristic length scales, velocities, and thermal diffusivity vary throughout the radius of a star, the Rayleigh number and Reynolds number are not specified in a general sense for these simulations. We note that such nondimensional parameters can potentially affect convective heat transport (Ahlers & Xu 2001; Kerr & Herring 2000) and dynamo action (Simitev & Busse 2005). In MUSIC, numerical truncation errors contribute to both the viscosity and thermal diffusivity. In addition MUSIC simulations include an explicit thermal diffusivity related to the thermal conductivity in Eq. (3)so that the Prandtl number is everywhere less than one. The simulations possible with a code like MUSIC are still many orders of magnitude away in parameter space from the highly turbulent conditions likely to be found in realistic stellar convection zones. The LES results should therefore be viewed merely as indicators of the properties of realistic stellar flows.
Fig. 2 Timeaveraged radial profile, produced from a volumeweighted average in the angular direction, of specific entropy in simulations Hi4 and Hi5. A heavy vertical black line marks the boundary between the radiative and convection zones, determined from the radial profile of entropy and the Schwarzschild criterion produced by the onedimensional stellar model. A gray vertical line marks the inner radial boundary of the spherical shell at 0.21 of the stellar radius. 

Open with DEXTER 
3.2. Boundary conditions
Because the placement of boundaries and the choice of boundary conditions affect the physical outcome of a hydrodynamic simulation, we use boundary conditions that are targeted to maintaining the original radial profiles of density and temperature. Each simulation volume begins at 20° from the north pole, and ends 20° before the south pole. We impose periodicity 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.
When the nearsurface layers are included in the spherical shell and resolved using a spliced grid, in simulations Low9 and Hi5, a blackbody radiation law is used on the outer radial boundary. To radiate as a blackbody, the energy flux is allowed to vary as , where σ is the StefanBoltzmann constant and T_{s}(θ,t) is the temperature along the surface. This physically realistic boundary condition can only be effectively used when the steep temperature gradient near the surface is sufficiently resolved; otherwise, it results in artificially high cooling rates. In this work when a uniform grid is used, the energy flux and surface luminosity are held constant at the outer radial boundary at the correct value of the energy flux at that radius in the onedimensional stellar evolution calculation; the surface luminosity of our model for the young sun is 2.32 times the luminosity of our current sun. On the inner radial boundary the energy flux is always held constant at the value indicated by the onedimensional stellar evolution data.
To treat the density, we use the hydrostatic equilibrium boundary condition described by GrimmStrele et al. (2015). Integration of the equation for hydrostatic equilibrium produces two possibilities for the density stratification: Here j is a grid index in the boundary cells, with j = 0 indicating the last radial cell inside the spherical shell. Equation (4)assumes constant internal energy and constant radial acceleration due to gravity in the boundary cells. Equation (5)assumes a constant pressure scale height in the boundary cells. This type of boundary condition assumes an ideal gas at the boundaries.
When the boundary condition does not closely match the density stratification of the stellar model, an artificial boundary layer develops near the simulation boundary. When the inner radial boundary of the spherical shell is in the radiative zone, we impose a constant radial derivative on the density at this boundary. A constant derivative is a comparatively simple boundary condition, but for the young sun we find that it more accurately and smoothly maintains the density stratification than using the hydrostatic equilibrium boundary condition on the inner radial boundary. When the inner radial boundary of the spherical shell is in the convection zone, we impose the hydrostatic equilibrium boundary condition on the density using Eq. (5). We find that in the lower convection zone, this hydrostatic equilibrium boundary condition minimizes the size of the boundary layer that would otherwise develop.
At the outer radial boundary we find that the hydrostatic equilibrium boundary conditions produce results indistinguishable from simply assuming a constant radial derivative for the density. We apply a boundary condition on the density that maintains hydrostatic equilibrium using Eq. (4). The accuracy of a particular boundary condition depends on the placement of the boundary and the physics of the fluid flow near the boundary, which are dictated by the sphericalshell geometry and the type of star.
4. Results
4.1. Convective turnover time
The convective turnover time is a fundamental parameter produced in onedimensional stellar evolution studies and used in further modeling (e.g. Kim & Demarque 1996; Matt et al. 2015). It can be defined in a variety of nearly equivalent ways (e.g. Landin et al. 2010; Meakin & Arnett 2007). In this work we define a local convective turnover time at any position in our convection zone as (6)In this expression h_{p} is the pressure scale height as a function of position and time;  v  is simply the velocity magnitude. Our stellar hydrodynamics model does not include rotation or mean flows, and no mean flows result in our simulations; the velocity magnitude in the convection zone is a purely convective velocity. A global time scale τ_{global} can thus be defined from the local convective turnover time by averaging over the entire convection zone: (7)Here the integration is volumeweighted and V(r,θ) is a volume element. The integration covers the convection zone between r_{min}/R = 0.47 and r_{max}/R = 0.85, and the full angular extent of our spherical shells. This lower limit is chosen as the deepest point within the convection zone of the young sun model that appears to be unaffected by the changing dynamics near the boundary of the radiative zone. The upper limit is chosen as a point in the upper convection zone that is contained in all of our spherical shells; in the upper convection zone the contribution to the convective turnover time becomes small. In practice the variation of the pressure scale height in the angular direction and in time is small; an expression based on an average of the pressure scale height divided by the rootmeansquare (rms) of the radial velocity yields similar results to our definition in Eq. (7). However, capturing small variations is an objective of our simulations and this motivates the form of Eq. (7). Ultimately, the convective turnover time is timeaveraged over the full time span of the simulation (see Table 1) to produce a single number representative of convection in the star, i.e. τ_{conv} = ⟨ τ_{global} ⟩ _{t}, where the brackets ⟨ ... ⟩ _{t} indicate a time average.
The radial profiles of ⟨ τ_{loc} ⟩ _{t} from simulations Low2 and Hi2 calculated using a volumeweighted average in θ are shown in Fig. 3, and are compared with a result from the onedimensional Lyon stellar evolution code. The result from the onedimensional calculation uses a velocity produced by mixing length theory: ⟨ τ_{loc} ⟩ _{t}  _{1D} ≡ h_{p}/v_{MLT}. Twodimensional hydrodynamic calculations produce a radial profile with a shape similar to the onedimensional stellar evolution calculation. However, at the bottom of the convection zone ⟨ τ_{loc} ⟩ _{t} initially decreases more sharply than ⟨ τ_{loc} ⟩ _{t}  _{1D}. In the upper convection zone, the twodimensional hydrodynamic simulations display a radial profile of ⟨ τ_{loc} ⟩ _{t} that is flatter and has a higher value. The magnitude of ⟨ τ_{loc} ⟩ _{t} in the hydrodynamic calculations is smaller by approximately a factor of 5 than ⟨ τ_{loc} ⟩ _{t}  _{1D}. Because twodimensional simulations produce higher velocities than threedimensional simulations, we expect that the convective turnover time from threedimensional simulations will lie between the one and twodimensional values.
Fig. 3 Radial profiles of the local convective turnover time ⟨ τ_{loc} ⟩ _{t} averaged over the full simulation time in simulations Low2 and Hi2, compared with the value obtained from the onedimensional Lyon stellar evolution code after dividing by a factor of 5. The shaded areas indicate one standard deviation above and below each of these timeaveraged lines; the shaded areas are nearly identical for these two simulations. 

Open with DEXTER 
Table 1 summarizes the average value τ_{conv} = ⟨ τ_{global} ⟩ _{t}, as well as the standard deviation of τ_{global}(t). These statistics are calculated from data sampled at a fixed time interval throughout our statistically steady convection simulations. The standard deviation is used in this context as an error estimate on the mean. In this table the total simulation time, i.e. the span of time in steadystate convection for which each simulation has been followed, is also indicated in units of the convective turnover time. The values of τ_{conv} for simulations Low19 are slightly higher than for simulations Hi15. This reflects that lower velocities result from lower resolution simulations where numerical dissipation is larger. For both sets of simulations Low19 and Hi15, τ_{conv} is on the order of 10^{6} s. The onedimensional Lyon stellar evolution code produces a longer convective turnover time on the order of 5 × 10^{6} s.
To understand this difference between the timeaveraged values of τ_{conv}, it is useful to analyze the spread of τ_{global} sampled at a fixed time interval on the order of τ_{conv}/ 10^{3} during our simulations. Figure 4 shows the full spread of this data as a box plot. In the left panel of this figure, simulations Low19 are shown on the vertical axis, while the spread of τ_{global} data is shown on the horizontal axis. The high value produced by the onedimensional Lyon stellar evolution code is only occasionally reached. In the right panel of this figure, similar box plots are shown for higher resolution simulations Hi15. By comparing simulations with identical grid spacing, we eliminated the effect of resolution. We are then able to observe that, although different spherical shells of the star are covered, the data sets of τ_{global} strongly overlap. We also observe that simulation Low1, which includes the boundary at the bottom of the convection zone, but does not include the radiative zone, experiences a larger variation in time for τ_{global} than do simulations that include the radiative zone.
Fig. 4 Standard Tukey box plot of the global time scale τ_{global} sampled at a fixed time interval throughout the simulations (left) Low19 and (right) Hi15. A box plot (e.g. Spitzer et al. 2014; McGill et al. 1978) is designed to show characteristics of the spread of the data. The line in the middle of the box is the median. The box itself represents the middle 50% of the data, so that the edges are the 25th and 75th percentiles. The whiskers, i.e. error bars, represent the extent of the data when outliers are discounted. Outliers are represented by circles. For a Tukey box plot, outliers are defined as data that lie farther than 1.5 times the interquartile range from the box. The vertical black line at 5 × 10^{6} s marks the value of τ_{conv} that we obtain from mixing length theory. 

Open with DEXTER 
4.2. Radial structure and amplitude of convective velocities in the young sun
In our model of the young sun, steady convection is characterized by wide regions of upflow velocity, and thinner regions with faster downflows. In contrast to threedimensional convection, in two dimensions the convection rolls are all aligned. Figure 5 illustrates this structure in the radial velocity for simulations Hi1 and Hi2. Simulation Hi2 includes more of the nearsurface layers and more of the radiative zone than simulation Hi1. From this figure it is clear that Hi2 has developed smallscale convection in the nearsurface layers that is absent in simulation Hi1. The additional extent of the radiative zone permits the growth of more waves in the radiative zone of simulation Hi2. However, the structure of flows in the bulk of the convection zone are not visibly distinguishable between simulations Hi1 and Hi2.
Fig. 5 Instantaneous radial velocity in simulations of the young sun. The left panel shows a typical snapshot from simulation Hi1, which excludes the nearsurface layers and the lower radiative zone. The right panel is typical of simulation Hi2, which includes a large extent of the radiative zone. Smallscale convective motions are clearly visible in the nearsurface layers of simulation Hi2. These simulations have identical temperature and density profiles, grid spacing, and boundary conditions. They are both visualized during steady convection. The sole difference is the position of the simulation boundaries in the radial direction. The color schemes are identical: blue indicates an inward flow, while red indicates an outward flow. 

Open with DEXTER 
Radial positions of convection rolls are time dependent. Therefore it is meaningful to examine radial profiles of the time average of the rms velocity ⟨ v_{rms} ⟩ _{t} which can be conveniently linked with helioseismology data and predictions from mixing length theory (e.g. Miesch et al. 2012). These radial profiles of velocity can provide details about the dynamic radial structure of the young sun.
4.2.1. Coupling to a central radiative zone: impact on velocity
The left panel of Fig. 6 compares the timeaveraged rms velocity as a function of radius in the deep stellar interior for four simulations that include different extents of the radiative zone, but all exclude the nearsurface layers. In the portion of the radiative zone we consider, which spans approximately 0.1 <r/R< 0.42, discrete large amplitude waves are destabilized. These waves propagate in the angular direction and can survive for long times; similar waves are identified by Alvan et al. (2015) as possible gravity waves. Simulations Low4 and Low6, which include a large extent of the radiative zone, display waves of particularly large amplitude in the radiative zone. However the amplitude of these waves does not appear to consistently increase with the extent of the radiative zone that is included in the spherical shell for simulations Low19. In addition, waves of similarly large amplitude are not present in the radiative zone of simulations Hi2 or ExH2, which are higher resolution but otherwise identical to simulation Low4. The length of time of the simulation appears to be critical; the tens of convective turnover times of simulations Hi2 or ExH2 may not be a sufficient amount of time to observe the excitation of large amplitude waves. We speculate that the large amplitude of these waves may be produced during intermittent events. We also observe no clear correlation in our simulations between the amplitudes of the waves in the radiative zone and the amplitude of ⟨ v_{rms} ⟩ _{t} in the lower convection zone.
Although the extent of the radiative zone included in the simulations in the left panel of Fig. 6 is different, the radial profiles of ⟨ v_{rms} ⟩ _{t} are remarkably similar in the convection zone. The bottom of the convection zone, indicated by a heavy vertical line, is calculated from the entropy profile produced by the onedimensional Lyon stellar evolution code. Above this boundary, ⟨ v_{rms} ⟩ _{t} smoothly rises to a maximum in each simulation. This maximum in the timeaveraged rms velocity lies approximately in the range 200 m/s to 300 m/s for these four simulations. Throughout the rest of the convection zone, the timeaveraged rms velocity remains relatively flat. The nearsurface layers are omitted from these simulations and are not shown on this plot. An area containing one standard deviation above and below the timeaverage for simulation Low2 is shaded. Because the average values lie within a standard deviation of each other, we consider the dynamics of these four simulations to be indistinguishable. We conclude that, as long as the boundary at the bottom of the convection zone is included, the convective dynamics are largely unaffected by the extent of the radiative zone included in the spherical shell.
Fig. 6 Left: timeaveraged radial profile of rms velocity in four lowresolution simulations that exclude the nearsurface layers but include a different extent of the radiative zone. Right: timeaveraged radial profile of rms velocity in simulations Low7 and Low8, that include the nearsurface layers up to the surface at r/R = 1.0. The shaded areas indicate one standard deviation above and below the timeaveraged lines; in the left panel the shaded area for simulation Low2 is shown. A heavy vertical black line marks the boundary between the radiative and convection zones, calculated from the radial profile of entropy produced by the onedimensional stellar evolution calculation and the Schwarzschild criterion. Gray vertical lines mark the inner radial boundary of the sphericalshell simulation volume for each simulation shown. 

Open with DEXTER 
The right panel of Fig. 6 shows the radial profile of ⟨ v_{rms} ⟩ _{t} for simulations Low7 and Low8. These two simulations include the nearsurface layers and are identical except that Low7 does not include the boundary at the bottom of the convection zone. Eliminating this natural boundary clearly results in a significantly elevated profile of the timeaveraged rms velocity in the lower convection zone. We observe that when some extent of the radiative zone is included in the simulation, the radiative zone appears to act as an energetic buffer. Energy from convective motions that reach the radiative boundary may contribute to waves in the radiative zone. In simulation Low7, where the radiative boundary is not included, all ballistic motions are simply reflected back into the convection zone. Higher velocity structures in the convection zone result.
We conjecture that this result could vary for different stellar models. A peculiarity of our young sun model is that the energy flux changes with the radius. Because Low7 and Low8 use boundary conditions consistent with the radial profile of energy flux in the young sun model, the energy flux through their lower boundaries is not identical. Placing the boundaries in different physical zones could lead to different energy buildup.
There is a difference of approximately a factor of 5 between the timeaveraged rms velocity produced in our twodimensional hydrodynamics simulations and the mixing length theory velocity used in the onedimensional calculation. This observation is not independent of the factor of 5 found in our comparison of convective turnover times in Sect. 4.1. The pressure scale height is used to produce the convective turnover time from the rms velocity; however the pressure scale height is indistinguishable between one and twodimensional simulations over the times we examine.
4.2.2. Coupling to the nearsurface layers: impact on velocity
The left panel of Fig. 7 compares the amplitude of timeaveraged rms velocities in the stellar interior for two simulations that include the nearsurface layers, Low8 and Low9. These simulations use identical simulation volumes, but differ in how well the nearsurface layers are resolved; simulation Low9 uses the spliced grid where grid spacing decreases in the nearsurface layers, with the result that the grid spacing near the surface is approximately half the size of that in simulation Low8. In addition, simulation Low9 is allowed to radiate heat dependent on the local temperature on the boundary, i.e. in a nonspherically symmetric way, as discussed in Sect. 3.2. The different treatment of the nearsurface layers has a significant impact on convective velocities throughout the convection zone, in line with the predictions of Spruit (1997). The maximum of ⟨ v_{rms} ⟩ _{t} at the bottom of the convection zone is almost twice as large in simulation Low9 as in Low8. The shaded areas that indicate one standard deviation above and below the time averages for Low8 and Low9 do not overlap. Aside from the difference in velocity magnitude throughout the convection zone, the profiles of ⟨ v_{rms} ⟩ _{t} exhibit a similar shape in simulations Low8 and Low9, indicating that the radial structure of stellar convection is preserved. Although distinctly different stars are produced by simulations Low8 and Low9, the approximate shape of the radial profile of ⟨ v_{rms} ⟩ _{t} in the upper convection zone is similar. This trend of increasing rms velocity with radius can also be compared to results produced for the current sun by the ASH, MURaM, and STAGGER codes (see Fig. 4 of Miesch et al. 2012).
Fig. 7 Timeaveraged radial profile of rms velocity in simulations Low8 and Low9 (left) which include the nearsurface layers up to the surface at r/R = 1.0, and simulations Hi4 and Hi5 (right). The shaded areas indicate one standard deviation above and below each timeaveraged line. A heavy vertical black line marks the boundary between the radiative and convection zones, calculated from the radial profile of entropy and the Schwarzschild criterion. Gray vertical lines mark the inner radial boundary of the sphericalshell simulation volume for each simulation shown. 

Open with DEXTER 
The right panel of Fig. 7 compares the radial profile of timeaveraged rms velocity in simulations Hi4 and Hi5. Simulation Hi5 has the same treatment of the surface layers as simulation Low9, using a spliced grid to better resolve the temperature gradient in the nearsurface layers in order to allow energy to radiate as blackbody radiation. Simulation Hi4 is identical to simulation Low8, but uses double the gridsize. Although simulation Hi4 experiences a faster rise in rms velocity in the upper convection zone than simulation Hi5, a larger magnitude of timeaveraged rms velocity is observed throughout the lower convection zone for simulation Hi5 than for simulation Hi4. In comparison with the left panel of Fig. 7, the simulations in the right panel show that, in the lower convection zone, a higher velocity is generally produced when the grid spacing is reduced. This occurs because we use only numerical dissipation of the velocity, which is reduced with higher resolution. Aside from this difference in amplitude, we find that properties of the flow in the lower convection zone produce a comparable form for the radial structure of the overshoot region, for the resolutions we are able to investigate. This is reassuring for future efforts to numerically study overshooting at the lower boundary of the convection zone. The differences observed in the upper convection zone are understandable because the treatment of the surface clearly affects the character of surface convection.
Fig. 8 Left: overshooting layer in simulation Hi1 represented by the shaded area between the radii of r_{o,top} and r_{o,bot}. The heavy vertical line shows the position of the boundary between the radiative and convection zones, indicated by the entropy profile produced by the onedimensional stellar evolution calculation. The full line shows the enthalpy flux, scaled by the magnitude of the large negative peak. The dashed line shows the Schwarzschild discriminant S(r), scaled to appear on this graph. Right: Schwarzschild discriminant in simulations Hi4 and Hi5, compared with the value obtained from the onedimensional model of the young sun. 

Open with DEXTER 
4.3. Overshooting layer width
Largescale convective motions behave ballistically, terminating at a range of radial points near the boundary with the convectively stable radiative zone. This situation can be classified as either convective overshooting or convective penetration, a distinction based on the impact that the convective flow has on the physics of the radiative zone (Brummell et al. 2002; Maeder 2009; Zahn 1991). In our simulations, we use a realistic thermal diffusivity identical to that used in the onedimensional stellar evolution calculation and that has not been artificially enhanced. Because of this, our simulations are in the large Péclet number regime in the deep stellar interior, and penetrative convection can be expected at the lower boundary of the convection zone. Simulations at large Péclet numbers in deep stellar interiors have been reported by Meakin & Arnett (2007) and Viallet et al. (2015).
An aspect of the dynamic structure of a star that results from the dynamic coupling across the boundary between the radiative zone and convection zones is the width of the overshooting layer. The width of the overshooting layer is sometimes called an overshooting length, and is related to an overshooting parameter in onedimensional simulations (e.g. Schröder et al. 1997; Zhang & Li 2012; Renzini 1987). Overshooting establishes a layer between the stable radiative zone and the convection zone where convective motions mix into a more quiescent fluid. The physics of this layer has interesting implications for the stellar chemistry, stellar evolution, and magnetic field generation (Marik & Petrovay 2002; Marques et al. 2006; Tian et al. 2009). We do not address the physics of the overshooting layer in this work aside from evaluating its width. The overshooting layer width that we calculate here is a convenient measure that can be benchmarked for simulations in different spherical shells; it is not the same as the overshooting length used in onedimensional stellar evolution calculations. The width of the overshooting layer can only be defined from a longtime average of the dynamics, which quantifies the interaction between the convection and radiative zones.
To evaluate the sensitivity of the overshooting layer to the sphericalshell geometry and boundary conditions, we adopt one possible method for calculating a width for the overshooting layer described in Brun et al. (2011) and Browning et al. (2004). We define the overshooting layer width as (8)where r_{o,top} is the radius at the top of the overshooting layer, defined as the point where the time average of the enthalpy flux first changes sign. In Table 1, r_{o,top} is listed for each simulation. The bottom of the overshooting layer, r_{o,bot}, is defined as the radius in the radiative zone where the enthalpy flux becomes negligible. In practice r_{o,bot} is also the point where waves in the radiative zone begin to impact the radial profile of timeaveraged enthalpy flux. We define the enthalpy flux following Freytag et al. (1996)(9)where the enthalpy H = e + P/ρ is calculated in the standard way in terms of the pressure P, density ρ, and internal energy e. The second term in Eq. (9)subtracts the effect of any bulk mass flow, which is typically small. The enthalpy flux in the region surrounding the overshooting layer is illustrated in the left panel of Fig. 8 for simulation Hi1. Both r_{o,top} and r_{o,bot} are labeled in this figure. The Schwarzschild discriminant S(r) = ∇_{ad}−∇ is also shown in Fig. 8 for simulation Hi1, as well as for simulations Hi4 and Hi5. In the Schwarzschild discriminant, the socalled adiabatic gradient ∇_{ad} = ∂log T/∂log P_{ad} is calculated using the equation of state. The local gradient ∇ = ∂log T/∂log P is calculated from our hydrodynamic simulations and appropriately averaged. The Schwarzschild criterion amounts to the statement that the Schwarzschild discriminant, S(r), must be greater than zero for stability against convection (Osterbrock & Schwarzschild 1958; Lebovitz 1965). We see small deviations from the onedimensional profile in the Schwarzschild discriminant in these figures. We do not expect that these small deviations are due to convective penetration because of the limited time that the simulations are observed in comparison with the time scale for thermal evolution at this depth. The small differences possibly stem from twodimensional fluid effects and the low resolution of the fluid simulations compared to the onedimensional stellar evolution calculation.
The overshooting layer width, w_{o}, for each simulation is given in Table 1 as a percentage of the total stellar radius R. For simulations Low19, the overshooting layer width is typically about 4% of the young sun’s radius, or 0.21h_{p} at the boundary to the radiative zone; for the higher resolution simulations Hi15, it is approximately 14% of the young sun’s radius, or 0.76h_{p}. Intuitively, the width of the overshooting layer is linked to the velocity amplitude. The difference in the overshooting layer width that we observe in simulations of different resolution clearly reflects the different velocity of convection rolls in the convection zone, and indicates a higher level of interaction between these zones when higher velocities are present. This relationship is shown in Fig. 9; the degree of certainty for the largest overshooting layer widths recorded in this figure are low, because they correspond to highresolution simulations that were performed for comparatively short times. Beyond this link to local velocity amplitudes, the overshooting layer width appears to be independent of the sphericalshell geometry and boundary conditions.
When the overshooting layer width is larger, its growth is not centered around the boundary where the Schwarzschild criterion indicates stability. It begins substantially higher in the convection zone and expands proportionately less into the radiative zone. For Low19 we measure r_{o,bot}/R ≈ 0.41, while for simulations Hi15 this value is r_{o,bot}/R ≈ 0.39. We observe that the inner radial boundary of simulation Low1 is identical to r_{o,bot}; this likely affects the unusually large overshooting layer width in this simulation.
Fig. 9 Overshooting layer width defined in Eq. (8)vs. the maximum of the timeaveraged rms velocity above the radiative zone boundary in the overshooting layer. A local polynomial fit to the data is also shown. 

Open with DEXTER 
5. Summary and discussion
We have examined how the extent of a twodimensional spherical shell, and therefore the boundary conditions imposed on an interior convection zone, affect the characteristics of compressible convection. A spherical shell that extends deeper into the central radiative zone is used to analyze the impact of a physically motivated open boundary condition on the inner radius of the convection zone. A spherical shell that extends into the nearsurface layers is used to analyze the impact of an open boundary condition on the outer radius of the convection zone. In the nonlocal theory of stellar convection, coupling between different zones of a star is expected to strongly influence stellar dynamics.
The results of our numerical experiments show that indeed this kind of coupling can affect the amplitude of the convective velocities that evolve throughout the convection zone. The inclusion of the radiative zone as an inner radiative boundary on the convection zone produces a lower velocity of convection rolls throughout the convection zone. The presence of smallscale convection in modeled nearsurface layers also affects the amplitude of convection rolls throughout the convection zone. The combined effect of sphericalshell geometry and boundary treatment can result in a difference greater than a factor of two in our velocity data. Thus, even though the flows in our modeled nearsurface layers are not well resolved, we do observe a nonlocal effect in agreement with Spruit (1997). To resolve the flows in the nearsurface layers well and to allow more accurate coupling with the convection zone, adaptive mesh refinement is a promising tactic; this is a development direction that we are pursuing. Along with an increase in the amplitude of convective velocities, a predictable increase in the width of the overshooting layer, and a decrease in the convective turnover time is generally observed.
The salient aspect of these results for our ongoing studies of threedimensional convective overshoot is that, although coupling between different zones changes the amplitude of convective velocities, the radial structure of the convective dynamics in the young sun maintains a strikingly similar profile, in particular in the overshoot region. Thus, a simulation that comprises accurate modeling of the nearsurface layers of the young sun will be able to produce results for deep convective processes physically similar to simpler models. These results will be relevant to a lower velocity regime, but can still contribute in critical ways to the 321D link. In general the results from simulations may be more realistic if the region of interest is far from the surface, such as fluid and chemical mixing, shearflows in a deep interior, or overshooting at the lower boundary of the convection zone.
This work has focused on analyzing the sphericalshell geometry and boundary conditions. In summary, we obtain two helpful observations: (1) that coupling between different zones in the young sun model changes the amplitude of dynamic quantities; but (2) the dynamic structure of the star is not affected. The simulations we analyze in this work are twodimensional and focus exclusively on largescale stellar flows. We therefore refrain from physical conclusions broader than the impact of the sphericalshell geometry and boundary conditions. A focused, more complete study of convection and convective overshooting in the young sun, including threedimensional simulations, is planned based on these results. A broader physical discussion of convection may await those calculations.
Acknowledgments
The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework (FP7/2007−2013)/ERC grant agreement No. 320478. M. Viallet is funded by the European Research Council though grant ERCAdG No. 341157COCO2CASA. 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 Supercomputer, a DiRAC Facility jointly funded by STFC, 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]
 Ahlers, G., & Xu, X. 2001, Phys. Rev. Lett., 86, 3320 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Alvan, L., Strugarek, A., Brun, A., Mathis, S., & Garcia, R. 2015, A&A, 581, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Alvarez Laguna, A., Lani, A., Deconinck, H., Mansour, N., & Poedts, S. 2016, J. Comput. Phys., 318, 252 [NASA ADS] [CrossRef] [Google Scholar]
 Arnett, D. W. 2014, Proc. IAU, 9, 459 [CrossRef] [Google Scholar]
 Arnett, W. D., & Meakin, C. 2009, Proc. IAU, 5, 106 [CrossRef] [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]
 Bodenschatz, E., Pesch, W., & Ahlers, G. 2000, Annu. Rev. Fluid Mech., 32, 709 [NASA ADS] [CrossRef] [Google Scholar]
 BöhmVitense, E. 1958, Z. Astrophys., 46, 108 [Google Scholar]
 Brandenburg, A. 2015, ArXiv eprints [arXiv:1504.03189] [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 [NASA ADS] [CrossRef] [Google Scholar]
 Brummell, N. H., Tobias, S. M., & Cattaneo, F. 2010, Geophys. Astrophys. Fluid Dyn., 104, 565 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A., & Palacios, A. 2009, ApJ, 702, 1078 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., Miesch, M. S., & Toomre, J. 2011, ApJ, 742, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Cameron, R., Voegler, A., & Schuessler, M. 2011, A&A, 533, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Canuto, V., & Dubovikov, M. 1997, ApJ, 484, L161 [NASA ADS] [CrossRef] [Google Scholar]
 Chacón, L. 2008, Phys. Plasmas, 15, 056103 [NASA ADS] [CrossRef] [Google Scholar]
 Cole, E., Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2014, ApJ, 780, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Cole, E., Brandenburg, A., Käpylä, P. J., & Käpylä, M. J. 2016, A&A, in press, DOI: 10.1051/00046361/201628165 [Google Scholar]
 Emran, M. S., & Schumacher, J. 2015, J. Fluid Mech., 776, 96 [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]
 Geroux, C., Baraffe, I., Viallet, M., et al. 2016, A&A, 588, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gilman, P. 1983, ApJS, 53, 243 [NASA ADS] [CrossRef] [Google Scholar]
 Glatzmaier, G. 2014, Introduction to Modeling Convection in Planets and Stars: Magnetic Field, Density Stratification, Rotation (Princeton University Press) [Google Scholar]
 Goffrey, T., Pratt, J., Viallet, M., et al. 2016, A&A, submitted [Google Scholar]
 Gough, D., & Weiss, N. 1976, MNRAS, 176, 589 [NASA ADS] [CrossRef] [Google Scholar]
 GrimmStrele, H., Kupka, F., LöwBaselli, B., et al. 2015, New Astron., 34, 278 [NASA ADS] [CrossRef] [Google Scholar]
 Grote, E., & Busse, F. 2001, Fluid Dyn. Res., 28, 349 [NASA ADS] [CrossRef] [Google Scholar]
 Heimpel, M., Aurnou, J., AlShamali, F., & Perez, N. G. 2005, Earth Planet. Sci. Lett., 236, 542 [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]
 Käpylä, P., Korpi, M., & Brandenburg, A. 2010, A&A, 518, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P., Mantere, M., Guerrero, G., Brandenburg, A., & Chatterjee, P. 2011a, A&A, 531, A162 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J., Mantere, M., & Brandenburg, A. 2011b, Astron. Nachr., 332, 883 [NASA ADS] [CrossRef] [Google Scholar]
 Kerr, R. M., & Herring, J. R. 2000, J. Fluid Mech., 419, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, Y.C., & Demarque, P. 1996, ApJ, 457, 340 [NASA ADS] [CrossRef] [Google Scholar]
 Knoll, D. A., & Keyes, D. E. 2004, J. Comput. Phys., 193, 357 [NASA ADS] [CrossRef] [Google Scholar]
 Landin, N., Mendes, L., & Vaz, L. 2010, A&A, 510, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Latour, J., Toomre, J., & Zahn, J.P. 1981, ApJ, 248, 1081 [NASA ADS] [CrossRef] [Google Scholar]
 Lebovitz, N. 1965, ApJ, 142, 229 [NASA ADS] [CrossRef] [Google Scholar]
 Lecoanet, D., Brown, B. P., Zweibel, E. G., et al. 2014, ApJ, 797, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Maeder, A. 2009, in Physics, Formation and Evolution of Rotating Stars (Heidelberg: Springer Berlin), Astron. Astrophys. Library 109 [Google Scholar]
 Magic, Z. 2016, A&A, 586, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Main, A., & Farhat, C. 2014, J. Comp. Phys., 258, 613 [NASA ADS] [CrossRef] [Google Scholar]
 Marik, D., & Petrovay, K. 2002, A&A, 396, 1011 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marques, J. P., Monteiro, M. J., & Fernandes, J. 2006, MNRAS, 371, 293 [NASA ADS] [CrossRef] [Google Scholar]
 Matt, S. P., Brun, A. S., Baraffe, I., Bouvier, J., & Chabrier, G. 2015, ApJ, 799, L23 [NASA ADS] [CrossRef] [Google Scholar]
 McGill, R., Tukey, J. W., & Larsen, W. A. 1978, Am. Stat., 32, 12 [Google Scholar]
 Meakin, C. A., & Arnett, D. 2007, ApJ, 667, 448 [NASA ADS] [CrossRef] [Google Scholar]
 Miesch, M. S., Elliott, J. R., Toomre, J., et al. 2000, ApJ, 532, 593 [NASA ADS] [CrossRef] [Google Scholar]
 Miesch, M. S., Featherstone, N. A., Rempel, M., & Trampedach, R. 2012, ApJ, 757, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Miesch, M., Matthaeus, W., Brandenburg, A., et al. 2015, Space Sci. Rev., 1 [Google Scholar]
 Mitra, D., Tavakol, R., Brandenburg, A., & Moss, D. 2009, ApJ, 697, 923 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, N. J., Brown, B. P., Brun, A. S., Miesch, M. S., & Toomre, J. 2011, ApJ, 739, L38 [NASA ADS] [CrossRef] [Google Scholar]
 Osterbrock, D. E., & Schwarzschild, M. 1958, Structure and Evolution of the Stars (Princeton, NJ: Princeton University Press) [Google Scholar]
 Pal, P. S., Singh, H. P., Chan, K. L., & Srivastava, M. 2007, Astrophys. Space Sci., 307, 399 [NASA ADS] [CrossRef] [Google Scholar]
 Park, H., Nourgaliev, R. R., Martineau, R. C., & Knoll, D. A. 2009, J. Comput. Phys., 228, 9131 [NASA ADS] [CrossRef] [Google Scholar]
 Quataert, E., & Gruzinov, A. 2000, ApJ, 539, 809 [NASA ADS] [CrossRef] [Google Scholar]
 Renzini, A. 1987, A&A, 188, 49 [NASA ADS] [Google Scholar]
 Roe, P. 1986, Annu. Rev. Fluid Mech., 18, 337 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, T. M., & Glatzmaier, G. A. 2005, 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]
 Roxburgh, L., & Simmons, J. 1993, A&A, 277, 93 [NASA ADS] [Google Scholar]
 Saikia, E., Singh, H. P., Chan, K., Roxburgh, I., & Srivastava, M. 2000, ApJ, 529, 402 [NASA ADS] [CrossRef] [Google Scholar]
 Schröder, K.P., Pols, O. R., & Eggleton, P. P. 1997, MNRAS, 285, 696 [NASA ADS] [CrossRef] [Google Scholar]
 Simitev, R., & Busse, F. 2005, J. Fluid Mech., 532, 365 [NASA ADS] [CrossRef] [Google Scholar]
 Simitev, R., Busse, F., & Kosovichev, A. 2011, ArXiv eprints [arXiv:1102.1092] [Google Scholar]
 Singh, H., Roxburgh, I., & Chan, K. 1995, A&A, 295, 703 [NASA ADS] [Google Scholar]
 Singh, H. P., Roxburgh, I. W., & Chan, K. L. 1998, A&A, 340, 178 [NASA ADS] [Google Scholar]
 Spitzer, M., Wildenhain, J., Rappsilber, J., & Tyers, M. 2014, Nature Methods, 11, 121 [CrossRef] [Google Scholar]
 Spruit, H. 1997, Mem. Soc. Astron. It., 68, 397 [NASA ADS] [Google Scholar]
 Strugarek, A., Brun, A., & Zahn, J.P. 2011, A&A, 532, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Thornber, B., Mosedale, A., Drikakis, D., Youngs, D., & Williams, R. 2008, J. Comput. Phys., 227, 4873 [NASA ADS] [CrossRef] [Google Scholar]
 Tian, C.L., Deng, L.C., & Chan, K.L. 2009, MNRAS, 398, 1011 [NASA ADS] [CrossRef] [Google Scholar]
 Trampedach, R. 2010, Ap&SS, 328, 213 [NASA ADS] [CrossRef] [Google Scholar]
 Trampedach, R., Stein, R. F., ChristensenDalsgaard, J., Nordlund, Å., & Asplund, M. 2014, MNRAS, 445, 4366 [NASA ADS] [CrossRef] [Google Scholar]
 Verhoeven, J., Wiesehöfer, T., & Stellmach, S. 2015, ApJ, 805, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Viallet, M., Baraffe, I., & Walder, R. 2011, A&A, 531, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Viallet, M., Baraffe, I., & Walder, R. 2013, A&A, 555, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Viallet, M., Meakin, C., Prat, V., & Arnett, D. 2015, A&A, 580, A61 [NASA ADS] [CrossRef] [EDP Sciences] [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, ZAp, 32, 135 [NASA ADS] [Google Scholar]
 Weiß, S., Seiden, G., & Bodenschatz, E. 2012, New J. Phys., 14, 053010 [NASA ADS] [CrossRef] [Google Scholar]
 Xie, X., & Toomre, J. 1993, ApJ, 405, 747 [NASA ADS] [CrossRef] [Google Scholar]
 Yadav, R. K., Gastine, T., Christensen, U. R., & Duarte, L. D. 2013, ApJ, 774, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Zahn, J.P. 1991, A&A, 252, 179 [NASA ADS] [Google Scholar]
 Zhang, Q., & Li, Y. 2012, ApJ, 746, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Ziegler, U., & Rüdiger, G. 2003, A&A, 401, 433 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
All Figures
Fig. 1 Radial structure of our young sun model, showing the profiles of the natural log of the pressure scale height h_{p}, the temperature T, and density ρ. The pressure scale height and temperature are scaled by arbitrary values in order to appear on the scale of this diagram. R is the total radius of the young sun. Average radial widths shown for physical zones are predicted by the young sun model calculated with the Lyon stellar evolution code. 

Open with DEXTER  
In the text 
Fig. 2 Timeaveraged radial profile, produced from a volumeweighted average in the angular direction, of specific entropy in simulations Hi4 and Hi5. A heavy vertical black line marks the boundary between the radiative and convection zones, determined from the radial profile of entropy and the Schwarzschild criterion produced by the onedimensional stellar model. A gray vertical line marks the inner radial boundary of the spherical shell at 0.21 of the stellar radius. 

Open with DEXTER  
In the text 
Fig. 3 Radial profiles of the local convective turnover time ⟨ τ_{loc} ⟩ _{t} averaged over the full simulation time in simulations Low2 and Hi2, compared with the value obtained from the onedimensional Lyon stellar evolution code after dividing by a factor of 5. The shaded areas indicate one standard deviation above and below each of these timeaveraged lines; the shaded areas are nearly identical for these two simulations. 

Open with DEXTER  
In the text 
Fig. 4 Standard Tukey box plot of the global time scale τ_{global} sampled at a fixed time interval throughout the simulations (left) Low19 and (right) Hi15. A box plot (e.g. Spitzer et al. 2014; McGill et al. 1978) is designed to show characteristics of the spread of the data. The line in the middle of the box is the median. The box itself represents the middle 50% of the data, so that the edges are the 25th and 75th percentiles. The whiskers, i.e. error bars, represent the extent of the data when outliers are discounted. Outliers are represented by circles. For a Tukey box plot, outliers are defined as data that lie farther than 1.5 times the interquartile range from the box. The vertical black line at 5 × 10^{6} s marks the value of τ_{conv} that we obtain from mixing length theory. 

Open with DEXTER  
In the text 
Fig. 5 Instantaneous radial velocity in simulations of the young sun. The left panel shows a typical snapshot from simulation Hi1, which excludes the nearsurface layers and the lower radiative zone. The right panel is typical of simulation Hi2, which includes a large extent of the radiative zone. Smallscale convective motions are clearly visible in the nearsurface layers of simulation Hi2. These simulations have identical temperature and density profiles, grid spacing, and boundary conditions. They are both visualized during steady convection. The sole difference is the position of the simulation boundaries in the radial direction. The color schemes are identical: blue indicates an inward flow, while red indicates an outward flow. 

Open with DEXTER  
In the text 
Fig. 6 Left: timeaveraged radial profile of rms velocity in four lowresolution simulations that exclude the nearsurface layers but include a different extent of the radiative zone. Right: timeaveraged radial profile of rms velocity in simulations Low7 and Low8, that include the nearsurface layers up to the surface at r/R = 1.0. The shaded areas indicate one standard deviation above and below the timeaveraged lines; in the left panel the shaded area for simulation Low2 is shown. A heavy vertical black line marks the boundary between the radiative and convection zones, calculated from the radial profile of entropy produced by the onedimensional stellar evolution calculation and the Schwarzschild criterion. Gray vertical lines mark the inner radial boundary of the sphericalshell simulation volume for each simulation shown. 

Open with DEXTER  
In the text 
Fig. 7 Timeaveraged radial profile of rms velocity in simulations Low8 and Low9 (left) which include the nearsurface layers up to the surface at r/R = 1.0, and simulations Hi4 and Hi5 (right). The shaded areas indicate one standard deviation above and below each timeaveraged line. A heavy vertical black line marks the boundary between the radiative and convection zones, calculated from the radial profile of entropy and the Schwarzschild criterion. Gray vertical lines mark the inner radial boundary of the sphericalshell simulation volume for each simulation shown. 

Open with DEXTER  
In the text 
Fig. 8 Left: overshooting layer in simulation Hi1 represented by the shaded area between the radii of r_{o,top} and r_{o,bot}. The heavy vertical line shows the position of the boundary between the radiative and convection zones, indicated by the entropy profile produced by the onedimensional stellar evolution calculation. The full line shows the enthalpy flux, scaled by the magnitude of the large negative peak. The dashed line shows the Schwarzschild discriminant S(r), scaled to appear on this graph. Right: Schwarzschild discriminant in simulations Hi4 and Hi5, compared with the value obtained from the onedimensional model of the young sun. 

Open with DEXTER  
In the text 
Fig. 9 Overshooting layer width defined in Eq. (8)vs. the maximum of the timeaveraged rms velocity above the radiative zone boundary in the overshooting layer. A local polynomial fit to the data is also shown. 

Open with DEXTER  
In the text 