EDP Sciences
Home All issues Volume 593 (September 2016) A&A, 593 (2016) A121 Full HTML
Free Access
Issue
A&A
Volume 593, September 2016
Article Number A121
Number of page(s) 12
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/201628296
Published online 04 October 2016

© 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 one-dimensional 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öhm-Vitense 1958; Abbett et al. 1997; Trampedach 2010; Brandenburg 2015). Nevertheless, non-local processes can significantly influence stellar convection (e.g. Spruit 1997; Gough & Weiss 1976; Canuto & Dubovikov 1997). Thus non-linear 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.

Whole-star 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 two-dimensional or three-dimensional, 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 spherical-shell geometry and boundaries on the dynamo. The goal of this work is to produce a similar study of the impact of the spherical-shell geometry on compressible convection without the fundamental influence of magnetic fields.

Although the impact of the spherical-shell 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 two-layer 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 two-dimensional 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, small-scale motions feed back on large-scale 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 non-local 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 near-surface 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 three-dimensional 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 near-surface layers. We study a prototypical model of a young, low-mass 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 pre-main 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 × 1011 cm of the total radius of 2.13 × 1011 cm. This large convective envelope allows us to study deep stellar convection, far from the physically complicated near-surface 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 non-uniform and are not linearly dependent on the radius. We define the near-surface layers as the portion of the star where the pressure scale height, hp = −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 Rayleigh-Bénard convection, which is driven by a uniform linear gradient of temperature, carefully controlled in a laboratory environment. In Rayleigh-Bé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 non-uniform and vary non-linearly 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 non-uniform gradients, regular patterns cannot be identified. Instead the diagnostics we use to evaluate the impact of the spherical-shell geometry and boundaries are velocity statistics gathered over long periods of steady convection.

thumbnail Fig. 1

Radial structure of our young sun model, showing the profiles of the natural log of the pressure scale height hp, 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σT3/ 3κρ is defined using the Stefan-Boltzmann 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 one-dimensional 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 low-mass 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 one-dimensional 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 one-dimensional stellar evolution models by studying critical short phases using stellar hydrodynamics in two and three dimensions.

Table 1

Parameters for two-dimensional 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 Courant-Friedrichs-Lewy (CFL) limit for time-explicit methods. Convergence of an implicit scheme can be computationally demanding. In the present study, the system of equations is discretized in time using a Crank-Nicolson scheme. To integrate the compressible hydrodynamic equations in time, a Jacobian free Newton-Krylov (JFNK) solver (Knoll & Keyes 2004) is employed. Instead of storing a Jacobian, a JFNK method uses matrix-vector products that can be estimated efficiently with finite difference methods. In MUSIC, a physics-based 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 semi-implicit 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 on-going 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 well-known 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 subgrid-scale 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 Rayleigh-Taylor instability, an ideal Kelvin-Helmholtz instability, and the Taylor-Green vortex. For these standard hydrodynamic problems MUSIC produces the expected results, and exhibits the expected numerical convergence properties.

3. Spherical-shell simulation volumes

3.1. Resolution

The compressible hydrodynamic Eqs. (1)−(3)are solved in a two-dimensional spherical shell using spherical coordinates: radius r and colatitude θ. In studying the impact of the spherical shell on compressible convection we examine only two-dimensional convection. Two-dimensional convection is known to produce higher velocity structures than three-dimensional 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 two-dimensional convection is the extreme case, and our results will also be relevant for three-dimensional convection simulations.

We consider the suite of two-dimensional 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 low-resolution simulations (Low1-8) use radial grid spacing Δr/R ≈ 2.8 × 10-3. The high-resolution simulations (Hi1-4) use radial grid spacing Δr/R ≈ 1.4 × 10-3. The extra-high-resolution 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 near-surface 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 near-surface layers, for r/R> 0.94, a non-uniform grid that is decreasing radially toward the surface is spliced on top of the uniform grid. The non-uniform portion of the grid is defined by a geometric sequence Δri = Δri + 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 Low1-8 and Hi1-4, respectively. The spliced grid of simulation Hi5 allows for a resolution of Δr/hp(r) ≈ 0.36 at the surface; for comparison, the uniform grid in simulation Hi4, has Δr/hp(r) ≈ 0.67 at the surface. For either the spliced grid or the uniform grid, our grid spacing in the near-surface 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 large-scale 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 high-resolution simulations for this benchmarking study of large-scale convection in a spherical-shell geometry, we are able to simulate over long times. Because the one-dimensional 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 near-surface layers is under-resolved. 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 non-dimensional 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.

thumbnail Fig. 2

Time-averaged radial profile, produced from a volume-weighted 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 one-dimensional 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 non-penetrative and stress-free boundary conditions on the radial boundaries.

When the near-surface 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 Stefan-Boltzmann constant and Ts(θ,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 one-dimensional 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 one-dimensional stellar evolution data.

To treat the density, we use the hydrostatic equilibrium boundary condition described by Grimm-Strele 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 spherical-shell geometry and the type of star.

4. Results

4.1. Convective turnover time

The convective turnover time is a fundamental parameter produced in one-dimensional 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 hp 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 volume-weighted and V(r,θ) is a volume element. The integration covers the convection zone between rmin/R = 0.47 and rmax/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 root-mean-square (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 time-averaged 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 volume-weighted average in θ are shown in Fig. 3, and are compared with a result from the one-dimensional Lyon stellar evolution code. The result from the one-dimensional calculation uses a velocity produced by mixing length theory: ⟨ τloc ⟩ t | 1D ≡ hp/vMLT. Two-dimensional hydrodynamic calculations produce a radial profile with a shape similar to the one-dimensional 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 two-dimensional 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 two-dimensional simulations produce higher velocities than three-dimensional simulations, we expect that the convective turnover time from three-dimensional simulations will lie between the one- and two-dimensional values.

thumbnail 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 one-dimensional Lyon stellar evolution code after dividing by a factor of 5. The shaded areas indicate one standard deviation above and below each of these time-averaged 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 steady-state convection for which each simulation has been followed, is also indicated in units of the convective turnover time. The values of τconv for simulations Low1-9 are slightly higher than for simulations Hi1-5. This reflects that lower velocities result from lower resolution simulations where numerical dissipation is larger. For both sets of simulations Low1-9 and Hi1-5, τconv is on the order of 106 s. The one-dimensional Lyon stellar evolution code produces a longer convective turnover time on the order of 5 × 106 s.

To understand this difference between the time-averaged values of τconv, it is useful to analyze the spread of τglobal sampled at a fixed time interval on the order of τconv/ 103 during our simulations. Figure 4 shows the full spread of this data as a box plot. In the left panel of this figure, simulations Low1-9 are shown on the vertical axis, while the spread of τglobal data is shown on the horizontal axis. The high value produced by the one-dimensional Lyon stellar evolution code is only occasionally reached. In the right panel of this figure, similar box plots are shown for higher resolution simulations Hi1-5. 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.

thumbnail Fig. 4

Standard Tukey box plot of the global time scale τglobal sampled at a fixed time interval throughout the simulations (left) Low1-9 and (right) Hi1-5. 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 × 106 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 three-dimensional 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 near-surface layers and more of the radiative zone than simulation Hi1. From this figure it is clear that Hi2 has developed small-scale convection in the near-surface 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.

thumbnail Fig. 5

Instantaneous radial velocity in simulations of the young sun. The left panel shows a typical snapshot from simulation Hi1, which excludes the near-surface layers and the lower radiative zone. The right panel is typical of simulation Hi2, which includes a large extent of the radiative zone. Small-scale convective motions are clearly visible in the near-surface 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 ⟨ vrms ⟩ 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 time-averaged 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 near-surface 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 Low1-9. 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 ⟨ vrms ⟩ 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 ⟨ vrms ⟩ 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 one-dimensional Lyon stellar evolution code. Above this boundary, ⟨ vrms ⟩ t smoothly rises to a maximum in each simulation. This maximum in the time-averaged 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 time-averaged rms velocity remains relatively flat. The near-surface layers are omitted from these simulations and are not shown on this plot. An area containing one standard deviation above and below the time-average 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.

thumbnail Fig. 6

Left: time-averaged radial profile of rms velocity in four low-resolution simulations that exclude the near-surface layers but include a different extent of the radiative zone. Right: time-averaged radial profile of rms velocity in simulations Low7 and Low8, that include the near-surface layers up to the surface at r/R = 1.0. The shaded areas indicate one standard deviation above and below the time-averaged 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 one-dimensional stellar evolution calculation and the Schwarzschild criterion. Gray vertical lines mark the inner radial boundary of the spherical-shell simulation volume for each simulation shown.

Open with DEXTER

The right panel of Fig. 6 shows the radial profile of ⟨ vrms ⟩ t for simulations Low7 and Low8. These two simulations include the near-surface 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 time-averaged 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 build-up.

There is a difference of approximately a factor of 5 between the time-averaged rms velocity produced in our two-dimensional hydrodynamics simulations and the mixing length theory velocity used in the one-dimensional 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 two-dimensional simulations over the times we examine.

4.2.2. Coupling to the near-surface layers: impact on velocity

The left panel of Fig. 7 compares the amplitude of time-averaged rms velocities in the stellar interior for two simulations that include the near-surface layers, Low8 and Low9. These simulations use identical simulation volumes, but differ in how well the near-surface layers are resolved; simulation Low9 uses the spliced grid where grid spacing decreases in the near-surface 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 non-spherically symmetric way, as discussed in Sect. 3.2. The different treatment of the near-surface layers has a significant impact on convective velocities throughout the convection zone, in line with the predictions of Spruit (1997). The maximum of ⟨ vrms ⟩ 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 ⟨ vrms ⟩ 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 ⟨ vrms ⟩ 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).

thumbnail Fig. 7

Time-averaged radial profile of rms velocity in simulations Low8 and Low9 (left) which include the near-surface 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 time-averaged 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 spherical-shell simulation volume for each simulation shown.

Open with DEXTER

The right panel of Fig. 7 compares the radial profile of time-averaged 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 near-surface layers in order to allow energy to radiate as blackbody radiation. Simulation Hi4 is identical to simulation Low8, but uses double the grid-size. Although simulation Hi4 experiences a faster rise in rms velocity in the upper convection zone than simulation Hi5, a larger magnitude of time-averaged 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.

thumbnail Fig. 8

Left: overshooting layer in simulation Hi1 represented by the shaded area between the radii of ro,top and ro,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 one-dimensional 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 one-dimensional model of the young sun.
Open with DEXTER

4.3. Overshooting layer width

Large-scale 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 one-dimensional 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 one-dimensional 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 one-dimensional stellar evolution calculations. The width of the overshooting layer can only be defined from a long-time average of the dynamics, which quantifies the interaction between the convection and radiative zones.

To evaluate the sensitivity of the overshooting layer to the spherical-shell 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 ro,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, ro,top is listed for each simulation. The bottom of the overshooting layer, ro,bot, is defined as the radius in the radiative zone where the enthalpy flux becomes negligible. In practice ro,bot is also the point where waves in the radiative zone begin to impact the radial profile of time-averaged 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 ro,top and ro,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 so-called 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 one-dimensional 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 two-dimensional fluid effects and the low resolution of the fluid simulations compared to the one-dimensional stellar evolution calculation.

The overshooting layer width, wo, for each simulation is given in Table 1 as a percentage of the total stellar radius R. For simulations Low1-9, the overshooting layer width is typically about 4% of the young sun’s radius, or 0.21hp at the boundary to the radiative zone; for the higher resolution simulations Hi1-5, it is approximately 14% of the young sun’s radius, or 0.76hp. 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 high-resolution 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 spherical-shell 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 Low1-9 we measure ro,bot/R ≈ 0.41, while for simulations Hi1-5 this value is ro,bot/R ≈ 0.39. We observe that the inner radial boundary of simulation Low1 is identical to ro,bot; this likely affects the unusually large overshooting layer width in this simulation.

thumbnail Fig. 9

Overshooting layer width defined in Eq. (8)vs. the maximum of the time-averaged 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 two-dimensional 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 near-surface layers is used to analyze the impact of an open boundary condition on the outer radius of the convection zone. In the non-local 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 small-scale convection in modeled near-surface layers also affects the amplitude of convection rolls throughout the convection zone. The combined effect of spherical-shell 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 near-surface layers are not well resolved, we do observe a non-local effect in agreement with Spruit (1997). To resolve the flows in the near-surface 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 three-dimensional 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 near-surface 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, shear-flows in a deep interior, or overshooting at the lower boundary of the convection zone.

This work has focused on analyzing the spherical-shell 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 two-dimensional and focus exclusively on large-scale stellar flows. We therefore refrain from physical conclusions broader than the impact of the spherical-shell geometry and boundary conditions. A focused, more complete study of convection and convective overshooting in the young sun, including three-dimensional 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 ERC-AdG No. 341157-COCO2CASA. 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 E-Infrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. DiRAC is part of the National E-Infrastructure. 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

All Tables

Table 1

Parameters for two-dimensional hydrodynamic simulations of the young sun.

In the text

All Figures

thumbnail Fig. 1

Radial structure of our young sun model, showing the profiles of the natural log of the pressure scale height hp, 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
thumbnail Fig. 2

Time-averaged radial profile, produced from a volume-weighted 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 one-dimensional 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
thumbnail 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 one-dimensional Lyon stellar evolution code after dividing by a factor of 5. The shaded areas indicate one standard deviation above and below each of these time-averaged lines; the shaded areas are nearly identical for these two simulations.
Open with DEXTER
In the text
thumbnail Fig. 4

Standard Tukey box plot of the global time scale τglobal sampled at a fixed time interval throughout the simulations (left) Low1-9 and (right) Hi1-5. 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 × 106 s marks the value of τconv that we obtain from mixing length theory.
Open with DEXTER
In the text
thumbnail Fig. 5

Instantaneous radial velocity in simulations of the young sun. The left panel shows a typical snapshot from simulation Hi1, which excludes the near-surface layers and the lower radiative zone. The right panel is typical of simulation Hi2, which includes a large extent of the radiative zone. Small-scale convective motions are clearly visible in the near-surface 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
thumbnail Fig. 6

Left: time-averaged radial profile of rms velocity in four low-resolution simulations that exclude the near-surface layers but include a different extent of the radiative zone. Right: time-averaged radial profile of rms velocity in simulations Low7 and Low8, that include the near-surface layers up to the surface at r/R = 1.0. The shaded areas indicate one standard deviation above and below the time-averaged 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 one-dimensional stellar evolution calculation and the Schwarzschild criterion. Gray vertical lines mark the inner radial boundary of the spherical-shell simulation volume for each simulation shown.

Open with DEXTER
In the text
thumbnail Fig. 7

Time-averaged radial profile of rms velocity in simulations Low8 and Low9 (left) which include the near-surface 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 time-averaged 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 spherical-shell simulation volume for each simulation shown.

Open with DEXTER
In the text
thumbnail Fig. 8

Left: overshooting layer in simulation Hi1 represented by the shaded area between the radii of ro,top and ro,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 one-dimensional 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 one-dimensional model of the young sun.
Open with DEXTER
In the text
thumbnail Fig. 9

Overshooting layer width defined in Eq. (8)vs. the maximum of the time-averaged 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

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.