© ESO, 2013

1. Introduction

Cool main-sequence stars such as the Sun have thick envelopes in which convective heat transport dominates the other energy transport mechanisms and determines the temperature structure. Although the atmosphere sitting on top of the convective envelope is convectively stable, the motions and the thermal inhomogeneity of the layers below extend into the optically thin layers and affect their 3D structure (see Nordlund et al. 2009, for a review). The convective motions are correlated with the local temperature and pressure fluctuations, which affect the line opacities and thus the spectral information obtainable from a star. Near-surface convection and the resulting spatial intensity and velocity patterns, which are not directly observable, thus have an impact on stellar spectra. This is in addition to their role in generating and structuring magnetic fields and in chemically mixing the envelope material. The knowledge of the 3D structure is therefore essential if one aims at determining stellar parameters (such as abundances, surface gravity, effective temperature, magnetic field strength and geometry, etc.) from spectra using inversion techniques (see, e.g., Collet et al. 2007; Wende et al. 2009).

Comprehensive three-dimensional (magneto-)hydrodynamic simulations have become an important tool to study solar and stellar convection. For the solar case, where spatially resolved observations are available, the simulations were found to be in excellent agreement with the measurements (e.g. Stein & Nordlund 1998; Nordlund et al. 2009; Pereira et al. 2013). The application of comprehensive 3D convection simulations to stars other than the Sun was pioneered by Nordlund & Dravins (1990a,b,c). In more recent years, 3D simulations of stellar convection were used for studies of abundances (e.g. Asplund et al. 1999; Collet et al. 2009), dust formation in very-cool stars (e.g. Freytag et al. 2010), later evolutionary stages (e.g. Freytag et al. 2002; Collet et al. 2007; Chiavassa et al. 2010), and the calibration of 1D methods to include 3D effects such as line broadening by micro- and macroturbulence (e.g. Wende et al. 2009) or the mixing-length parameter (e.g. Ludwig et al. 1999; Trampedach 2007). A grid of 37 stellar 3D simulations of main-sequence and giant stars was recently presented by Trampedach et al. (2013). First simulations including magnetic fields in stellar near-surface convection were performed by Beeck et al. (2011) and Wedemeyer et al. (2013).

In this paper, we consider MURaM simulations of non-magnetic convection in six main-sequence stars for a systematic parameter study of near-surface convection in cool main-sequence stars. The simulation setup was chosen as similar as possible for the six stars, in order to facilitate a comparison between the simulation results. Here, we mainly focus on the general characteristics of the non-magnetic convection and atmospheric structure, such as flow geometries, energy fluxes, and mean profiles of temperature, density, and other quantities. In subsequent papers of this series, we will study the morphology of the granulation patterns and its influence on line shapes and the effects of magnetic fields.

2. Simulations

2.1. The MURaM code

The MURaM (MPS/University of Chicago Radiative MHD) code was developed by the MHD simulation groups at the Max Planck Institute for Solar System Research (MPS) in Katlenburg-Lindau and at the University of Chicago (Vögler 2003; Vögler et al. 2005). It solves the three-dimensional, time-dependent MHD equations for a compressible and partially ionised plasma. The non-grey radiation transport included in the code is based on the short-characteristics scheme (Kunasz & Auer 1988) and uses opacity binning (Nordlund 1982) based on opacity distribution functions computed with ATLAS9 (Kurucz 1993). For the simulation of the G2V star, the opacity bins previously obtained for the solar MURaM simulations (see Vögler et al. 2004) were applied. The opacity bins for the G2V simulation were used in low horizontal resolution simulations to obtain the approximate stratification for the K0V and F3V stars. The τ-sorting method was applied to horizontal averages of these atmospheres to produce a better estimate for the opacity bins. Analogously, the bins of the K0V star were then used as a starting point for determining opacity binning for the K5V and M0V stars, and the opacity bins thus obtained for the M0V star were used as the initial choice for the M2V star. The procedure for each star could be iterated, however we found that the changes caused by a second iteration were minimal. For a detailed description of the radiation transfer in the MURaM code, see Vögler et al. (2004).

The OPAL equation-of-state (Rogers 1994; Rogers et al. 1996) for a plasma with a solar composition (Anders & Grevesse 1989) was used to generate look up tables for temperature, pressure, and entropy as a function of density and internal energy density. The MURaM code uses discretised spatial derivatives (4th-order centred differences) on a three-dimensional Cartesian grid and explicit time stepping (4th-order Runge-Kutta solver). In order to stabilise the scheme, artificial diffusivities are introduced. The code version used for the simulations discussed here employed a Minmod slope limiter which sets the numerical diffusivity of the scheme (see Rempel et al. 2009). The bottom boundary was open to vertical inflows and arbitrarily inclined outflows. The pressure at the bottom boundary was adjusted to maintain a fixed mass in the computational domain. For the simulations dicussed here, the entropy density of the material entering the domain through the bottom boundary was varied in an iterative process until the desired value of the effective temperature was (approximately) reached. In all six simulations, the top boundary was closed for flows. Test simulations with open top boundary were also run and yielded no significantly different results. The code uses periodic side boundary conditions.

2.2. Stellar parameters

For the near-surface layers and atmosphere of a cool star, the governing parameters are the gravitational acceleration, g, the effective temperature, T_eff, and the chemical composition. We use solar abundances in all cases. The effective temperature and the gravitational acceleration were chosen to match the conditions in cool main-sequence stars. We have carried out simulations corresponding to the following spectral types: F3V, G2V, K0V, K5V, M0V, and M2V (stellar parameters given in Table 1). Figure 1 shows the location of the six models in the log g-log T_eff plane along with three isochrones marking the position of the main sequence (Bressan et al. 2012).

While gravitational acceleration and chemical composition explicitely enter the simulations as parameters, the effective temperature is indirectly specified through the bottom boundary condition of the code, see Sect. 2.1. For the analysis presented here, the simulations have been run long enough with fixed inflowing entropy density for any transients to dissappear, however T_eff varies slightly due to oscillations and granulation. The standard deviation of the temporal fluctuations of T_eff is given in Table 1.

2.3. Simulation setup

The dimensions of the computational domain (“local box”) were adapted according to the stellar parameters: the height of the box was chosen to contain about 13 to 15 pressure scale heights (at least six below and six above the optical surface). The vertical cell size (height resolution Δz) was set sufficiently small to resolve steep temperature gradients and to maintain Δz < H_p/5 at every depth, where H_p is the local pressure scale height. For most of the models, 300 or less cells in the vertical direction were sufficient to meet these criteria. Only the F3V model with strongly varying pressure scale height and a very steep local temperature gradient, required 800 cells in the vertical direction. The two equal horizontal box dimensions were scaled to the expected granule (convection cell) size. In order to reduce the effects of the periodic boundary conditions and obtain good statistics while maintaining sufficient spatial grid resolution, the boxes were chosen big enough to contain 30 to 50 granules at any given time. The horizontal dimensions were resolved into 512 × 512 grid cells. The largest computational domain of the simulations thus comprised 512 × 512 × 800 ≈ 2.1 × 10⁸ cells. Table 2 gives a summary of the computational box dimensions and grid resolutions.

3. Results

3.1. General morphology of near-surface convection

Figure 2 gives maps of the bolometric intensity emerging vertically from the simulated stellar surfaces for single snapshots of the time-dependent simulations. All simulations show intensity patterns reminiscent of solar granulation. The typical size of the granules varies from ~5 Mm for F3V to ~0.3 Mm for M2V. The rms bolometric intensity contrast (denoted by σ_I in Table 1) decreases from about 20% for F3V to less than 3% in the M2V simulation, reflecting decreasing temperature fluctuations on surfaces of constant optical depth (see Sect. 3.3).

There are qualitative changes in the visual appearance of the surface convection along the sequence of simulated stars. For instance, the granulation pattern of the F3V model appears “rough” and irregular owing to numerous shock waves at the optical surface. Shocks are rarer and weaker in the near-surface layers of the cooler stars since the typical convective velocities are lower (also in relation to the sound speed; cf. Fig. 6). At the cool end of our model sequence, the M-dwarf granules, which are sustained by the slowest convective flows, have more irregular shapes but less brightness substructure than their counterparts on the simulated G- and K-type stars. As we report quantitatively in Paper II, their dark intergranular lanes are thinner (with respect to the granule size) and vary more strongly in intensity and width than those of the other stars (see also Ludwig et al. 2002).

	Fig. 1 Stellar parameters of the six models along with three isochrones by Bressan et al. (2012), solid line: zero-age main sequence (ZAMS), dashed line: age of 1 Ga, dotted line: age of 4.5 Ga (approximate solar age) on the log g–log Teff plane.
Open with DEXTER

	Fig. 2 Maps of bolometric intensity emerging vertically from single snapshots of the six simulations. The grey scale of each image is saturated at ±2 σI, where σI is the rms contrast (cf. Table 1). The significant difference in the length scales of the images is illustrated by the inset in the upper left panel, which shows all other images on the same scale as the one from the F3V simulation.
Open with DEXTER

Ludwig et al. (2006) found “dark knots” associated with strong downflows and vortex motion in simulations of convection in M-type main- and pre-main-sequence objects. Our simulations show knots of high vorticity associated with strong downflows in all models (some examples in Fig. 2 are: G2V, (x,y) = (8.7 Mm,4.4 Mm); K5V, (x,y) = (0.36 Mm,0.52 Mm); M0V, (x,y) = (0.45 Mm,0.6 Mm)). They become increasingly stable and prominent at lower effective temperatures. In our models, some of these vortices are evacuated strongly enough by the effect of the centrifugal force to become brighter than their surroundings (cf. vortices in solar simulations studied by Moll et al. 2011, 2012). Most frequently these bright vortex structures occur in our two K-type simulations.

A more detailed analysis of the granulation properties and their effects on spectral lines is given in Beeck et al. (2013, Paper II, hereafter).

3.2. Velocity field

As the visible granulation pattern is created by convective flows, it is strongly correlated to the vertical velocities at the optical surface, υ_z(z = 0). Figure 3 shows υ_z(z = 0) for four of the six simulations. The snapshots are taken at the same time as in Fig. 2. The colour scale of the images saturate at 2 υ_z,rms(z₀) with z₀: = ⟨z⟩_τR = 1, values of which are given in Table 1. The granules visible in Fig. 2 correspond to upflows, while the dark intergranular lanes correspond to downdrafts. In the G-, K-, and M-type simulations, an anti-correlation between size and mean upflow velocity of the granules is indicated: while most of the small convection cells appear (almost) saturated in Fig. 3, meaning their velocity reaches 2 υ_z,rms(z₀), the larger granules appear paler, meaning their upflow speed is lower. In the F3V simulation, this effect is not visible, due to a strong large-scale modulation of the vertical velocity at the optical surface. This large scale pattern might hint to a strong mesogranulation in this spectral type. Unfortunately, the length scale of this modulation is the horizontal box size, which raises the question whether this effect is produced, enhanced, or modified by the periodic boundary condition. A test simulation with a box twice as large has shown similar but weaker large-scale modulation but probably is still strongly influenced by the periodic boundary condition.

Figure 4 shows maps of the vertical velocity at a depth of 4.6 pressure scale heights below the optical surface, where the average pressure is 100 times the average pressure at the optical surface, p₀: = ⟨p⟩_τR = 1. The typical size of the convection cells is significantly larger at this depth than at the surface. A rough estimate based on mass conservation and stationarity predicts a proportionality between the horizontal scale of the vertical velocity pattern at a given depth, D_hor, and the local density scale height, H_ϱ (Nordlund et al. 2009), (1)where υ_hor and υ_ver are the horizontal and vertical convection velocities, respectively. They can be approximated by the horizontally averaged (height-dependent) rms values of the vertical and horizontal components of the fluid velocity, and respectively (for the definition and discussion of the horizontal average ⟨ · ⟩ _z, see Appendix A). In the left panel of Fig. 5 the ratio of υ_x,y,rms and υ_z,rms is plotted as a function of normalised gas pressure ⟨p⟩_z/p₀. In subsurface layers, we find υ_hor/υ_ver ≈ 1, hence Eq. (1) predicts that the horizontal scale of the flow pattern roughly follows the trend of the inwardly increasing density scale height.

The right panel of Fig. 5 shows the profile of D_hor as derived from Eq. (1) in units of the horizontal box size X_tot of the respective simulation. The density scale height obtained from the simulations was smoothed (convolution with a Gaussian kernel σ = 10Δz) to avoid a sharp maximum of D_hor at the optical surface of the two hottest models (cf. Fig. 9). The predicted horizontal scale of 12–20% of the horizontal box size at the optical surface and 25–40% of the horizontal box size at p = 100·p₀ matches the sizes of the patterns visible in Figs. 3 and 4.

	Fig. 3 Maps of the vertical velocity υz at constant geometrical depth z0 = ⟨z(τR = 1) ⟩ (average level of the optical surface) for four of the six models. Upward motions are blue, downward motions are red, colour scales saturate at ± 2·υz,rms(z0) (for values, see Table 1). Note that the horizontal scales are different (cf. Fig. 2).
Open with DEXTER

	Fig. 4 Same as Fig. 3, but for geometrical depth z2 with ⟨p(z2) ⟩ = 100 ⟨p(z0) ⟩, corresponding to 4.6 pressure scale heights below the optical surface.
Open with DEXTER

Fig. 5 Properties of the convective flows. Left: ratio of rms values of horizontal and vertical flow velocity on surfaces of constant geometrical depth. Right: estimated horizontal scale of the vertical velocity patterns as derived from Eq. (1) on surfaces of constant geometrical depth. The horizontal scale is given in units of the horizontal box size, Xtot, for an easier comparison with Figs. 3 and 4. The solid curves refer to the four simulations shown in these figures, the dashed curves to the remaining two simulations.

Open with DEXTER

	Fig. 6 Flow velocity rms. Left: rms of the vertical component of the flow velocity on surfaces of constant geometrical depth. Right: rms of the modulus of the flow velocity in units of the local sound speed, cs, (Mach number) on surfaces of constant geometrical depth.
Open with DEXTER

	Fig. 7 Properties of up- and downflows. Left: relative area covered by upflows (υz > 0) on surfaces of constant geometrical depth as functions of normalised averaged pressure. Right: average speed of the upflows, υu: = ⟨υz|υz > 0 ⟩ z, (solid) and of the downflows, υd: = ⟨υz|υz < 0 ⟩ z, (dashed) as functions of normalised average pressure for four of the six simulations.
Open with DEXTER

The left panel of Fig. 6 shows the depth dependence of υ_z,rms, which is a measure of the typical convective velocity. The profiles of υ_z,rms all peak near the optical surface, where radiative energy transport starts to become important (cf. Fig. 13). The peak rms velocity decreases with decreasing effective temperature. The position of the maximum of υ_z,rms shifts along the model sequence: in the F- and G-star simulations this maximum is almost directly at the optical surface while in the cooler models it is about one to two pressure scale heights below it. In the deeper layers, the convection velocity decreases monotonically with increasing depth in all simulations. In the optically thin upper layers, the overshooting large-scale convective motions slow down with increasing height above the optical surface where the stratification is stable against convection. However, the υ_z,rms drop only for about one to three scale heights, before they rise again, as shocks become more important.

The right panel of Fig. 6 shows profiles of the mean Mach number. Although the typical velocities in all simulations reach a substantial fraction of the local sound speed c_s, only in the atmosphere of the F3V simulation is an average Mach number of order unity reached. Surface convection is largely subsonic in our simulations of M, K, and G stars.

The left panel of Fig. 7 shows the depth dependence of the relative area of the upflows plotted as functions of normalised average pressure. Below the surface layers, the upflow area is very similar in all six simulations and almost constant at about 63 to 65% of the total area, which reflects the asymmetry between fast, dense downflows and slower upflows. The value for the upflow area of approximately 2/3 of the total area is in good agreement with the results of Trampedach & Stein (2011) who used another code and different stellar parameters. Near the optical surface, the area fraction of the upflows drops to about 50% as the strong correlation between vertical velocity and density weakens and the asymmetry between up- and downflows decreases. The low value of the relative upflow area in the upper layers of the F3V simulation of about 42 to 44% can be interpreted as an effect called reversed granulation in the subadiabatic atmospheric layers (Cheung et al. 2007, and references therein), which inverts the correlation between υ_z and ϱ. This effect is amplified by shock fronts in the F3V simulation: the material trailing the shocks which move upwards is over-dense compared to the average stratification.

The right panel of Fig. 7 shows the mean speed of the upflows and the downflows. Although both speed profiles peak slightly below the optical surface in all simulated stars, the asymmetry between up- and downflows in the convectively unstable layers leads to mean downflow speeds reaching 1.6 to 1.8 times the mean upflow speed in the lower part of the box. In the atmospheres, averaged up- and downflow speeds are almost equal for the four cooler stars. In the F3V simulation, the presence of shocks leads to a steep rise of the average upflow speed in the optically thin layers, whereas the gradient of the averaged downflow speed is flatter (indicating an asymmetry reversed to the one observed below the surface).

3.3. Temperature, pressure, and density

	Fig. 8 Temperature and pressure stratifications on surfaces of constant optical depth. Top panel: temperature averaged over surfaces of constant optical depth. Middle panel: same as top panel, but normalised by the effective temperature of the respective model. Bottom panel: pressure averaged on surfaces of constant optical depth.
Open with DEXTER

	Fig. 9 Horizontally averaged local scale heights as functions of normalised pressure. Left: pressure scale height. Right: density scale height.
Open with DEXTER

	Fig. 10 Profiles of mean temperature ⟨T⟩z (left panel) and mean density ⟨ϱ⟩z (right panel) as functions of normalised pressure.
Open with DEXTER

	Fig. 11 Profiles of the logarithmic temperature gradient (left panel) and of superadiabaticity (right panel) averaged on iso-τR surfaces as function of pressure.
Open with DEXTER

The upper two panels of Fig. 8 show the temperature (averaged over iso-τ_R surfaces; for a discussion of the different averages ⟨ · ⟩ _τ, ⟨ · ⟩ _z, and ⟨ · ⟩ _p, see Appendix A) as function of log τ_R, both in absolute units and normalised by T_eff. The simulations from F3V to K5V show a steep temperature gradient just beneath the optical surface, whereas, in the M dwarfs, the steepest temperature gradient occurs well below this layer (see also Figs. 10 and 11 and Sect. 3.4). In the normalised representation, all simulations have a similar profile in the atmosphere (log τ_R < 0), while their temperature curves diverge in the subphotospheric layers.

The bottom panel of Fig. 8 gives the pressure (averaged on iso-τ_R surfaces) plotted as function of log τ_R. In the atmosphere, where opacity and temperature are only mildly height-dependent, log p essentially depends linearly on log τ_R as the structure is governed by hydrostatic and radiative equilibrium and is also almost iso-thermal. In the layers just below the photosphere, the curves for the different simulations diverge. The diverging profiles of the subphotospheric temperature (middle panel) and pressure (bottom panel) reflect that the pressure and temperature structures are determined by convection below the photosphere and converge to different adiabat (depending on the stellar parameters) in the deep convective envelopes. Figure 9 shows the depth dependences of the pressure and density scale heights as functions of normalised pressure. We define horizontally averaged local scale heights as (2)As the gravitational acceleration increases and the photospheric temperature decreases monotonically from F3V to M2V the, local pressure scale height around τ_R = 1 decreases from ~500 km in F3V to ~35 km in M2V. In the atmosphere, where the temperature is mildly height-dependent, the local pressure scale height becomes roughly constant. In the convective layers, the strong temperature gradient entails also a strong increase of the pressure scale height towards deeper layers. At p = 100 p₀, near the bottom of the simulation boxes, the local pressure scale height of the F3V simulation is already ~2000 km, which poses a problem for the current implementation of the MURaM code with its fixed vertical cell size (high computational costs).

The peak of the density scale height near the optical surface in some simulations coincides with the strong photospheric temperature gradient of these simulations. Locally, the density scale height often becomes negative at the optical surface in the F3V and G2V simulations (density inversion). In the subsurface layers with high temperature gradient, the density scale heights are somewhat (M2V-G2V:15–30%, F3V: up to 45%) larger than the pressure scale heights whereas, in the almost isothermal atmospheres, the scale heights of pressure and density are almost equal.

Figure 10 shows the profiles of temperature and density averaged on iso-z surfaces as functions of the normalised pressure. As already seen in Fig. 8 on the τ_R scale, the coolest models lack the strong photospheric temperature gradient of the warmer models. This is shown quantitatively in the left panel of Fig. 11, which gives the mean profiles of the logarithmic temperature gradient, ⟨∇⟩_τ = d log ⟨T⟩_τ/d log ⟨p⟩_τ. Here, the iso-τ average was chosen because ∇ changes considerably near the optical surface. Since this represents a transition from a (highly) superadiabatic to a subadiabatic regime, averaging over iso-z planes would smear out the sharp photospheric feature in the temperature gradient and thus obscure the relevant physics in this layer.

In the right panel of Fig. 11, the profile of the superadiabaticity ⟨∇⟩_τ− ⟨∇_ad⟩_τ is given. The superadiabaticity in the lowest part of the simulation domain is small (~10^-3) for most of the models, with the exception of the F3V simulation, where the stratification remains substantially superadiabatic even 5 pressure scale heights below the optical surface. For the M dwarfs, the superadiabaticity is low compared to the hotter models, even in the layers directly beneath the optical surface. This is a consequence of the high densities (i.e. high heat capacity per volume) and low energy fluxes and consequently low horizontal temperature fluctuations.

	Fig. 12 Profiles of relative rms fluctuations. Top: temperature fluctuations on surfaces of constant optical depth. Middle: gas pressure fluctuations on surfaces of constant optical depth. Bottom: pressure fluctuations on surfaces of constant geometrical depth.
Open with DEXTER

The top panel of Fig. 12 shows the relative rms fluctuations of temperature on surfaces of constant optical depth. The relative rms fluctuations of temperature show a monotonic decrease from the hotter to the cooler stars at all depths. This is consistent with the trend in the bolometric intensity contrast (cf. Table 1). The subphotospheric peak in the depth profile of temperature fluctuations is at lower optical depth in the F- and G-star simulations compared to the cooler simulations. This has already been pointed out by Nordlund & Dravins (1990a), who coined the expression “hidden” or “veiled” granulation for stars cooler than the Sun, as the maximum temperature contrast occurs far below the optical surface. In the case of the K0V star, the relative temperature contrast at the optical surface is only about 34% of its peak value at log τ_R ≈ 1 (compared to 42% in the solar simulation and 64% in the F3V simulation). The reason for this effect is the lower temperature-sensitivity of opacity near the optical surfaces of the cooler K and M stars. This leads to the transition from convective to radiative energy transport occuring at somewhat larger optical depth or normalised pressure (particularly in the K-star simulations) and over a larger optical depth range or normalised pressure range (particularly in the M-star simulations; see Sect. 3.4).

The middle panel of Fig. 12 shows the rms fluctuations of gas pressure on surfaces of constant optical depth. They also diminish with decreasing effective temperature of the simulations, with the notable exception of the G-type star, where the rms fluctuations of pressure on surfaces of constant optical depth are very low in the upper atmosphere compared to the much cooler K stars. This can be explained as an opacity effect: while the temperature dependence of the Rosseland opacity κ_R(p,T) is usually much more important than the pressure dependence, in the temperature range between 4000 and 5000 K and at pressures between 10² and 10⁵ dyn cm^-2, κ_R is nearly independent of temperature. Moreover, the temperature fluctuations in the atmosphere are in general relatively small, so that density depends mainly on pressure. Therefore, the increment of the optical depth d τ_R = (κ_Rϱ) dz becomes (almost) independent of temperature, too, so that iso-τ_R and iso-p surfaces of an atmosphere in this temperature regime are almost identical, i.e. the fluctuations of pressure on surfaces of constant optical depth become small. In our simulation sequence, the G2V star is the only star where the temperature and pressure of the atmospheric layers at − 4 ≤ log τ_R ≤ −1 fall in the regime of nearly temperature-independent κ_R.

The bottom panel of Fig. 12 shows the rms fluctuations of pressure on planes of constant geometrical depth, which are mainly monotonically decreasing with increasing depth and with decreasing effective temperature. This illustrates that in the deeper layers, where velocities and horizontal temperature fluctuations are small, the deviations from hydrostatic equilibrium are small, too. Horizontal temperature fluctuations entail horizontally varying local pressure scale heights: the pressure fluctuations on a horizontal plane in the atmosphere can be regarded as the integrated effect of the temperature fluctuations below this plane.

	Fig. 13 Profiles of average energy flux, normalised to total flux. Left: convective energy flux split into net enthalpy flux and net kinetic energy flux. Right: radiative flux.
Open with DEXTER

Before the energy balance of the simulated stellar surface layers is analysed in Sect. 3.4, we want to point out that the trends observed in convective velocities and temperature fluctuations can be consistently explained as an effect of the stellar parameters: the lower temperature and higher gravitational acceleration in the atmospheres of cooler stars result in much higher densities. The resulting higher heat capacity per unit volume, c_p/ϱ, and the much lower net energy flux () then have the consequence that the convective motions in the coolest models are less vigorous. This is reflected in the rms velocities of about 0.5 km s^-1 for those simulations compared to more than 5 km s^-1 for the F3V-star simulation (see Fig. 6). A rough estimate of the convective heat flux F_conv (enthalpy flux) in the spirit of mixing-length theory gives (3)where ΔT is the temperature contrast between up and downflows, υ_conv is the convective velocity, and ϱ is the density, both of which are assumed to be equal in up- and downflows for this rough estimate. In the convective subsurface layers of cool stars radiative energy transport can be neglected, so that If we approximate the quantities υ_conv and ΔT by ⟨υ_z,rms ⟩ and ⟨δT⟩, respectively ( ⟨ · ⟩ denotes a horizontal and temporal average), and if we further ignore the variation of c_p with temperature due to ionisation, we obtain: (4)Directly beneath the optical surface, the product on the left-hand side decreases by a factor of about 500 from F3V to M2V, which entails a strong variation of ⟨δT⟩ · ⟨υ_z,rms(z₀) ⟩. In fact, the simulations show that both ⟨δT⟩ and ⟨υ_z,rms(z₀) ⟩ decrease monotonically through the model sequence by about an order of magnitude each (see Figs. 5 and 12).

3.4. Energy flux

The total energy flux leaving the stellar atmosphere in the form of radiation is supplied by energy flux from below. In the absence of nuclear energy sources, the (temporally averaged) total luminosity is constant throughout the upper layers of a star. In plane-parallel geometry, this means that the energy flux is independent of depth.

In the convection zone, the energy flux is almost entirely provided by convective energy transport and is mainly composed of two opposing fluxes: the net enthalpy flux which is directed towards the surface, and the net kinetic energy flux, which is directed inwards. This is a consequence of the asymmetry in temperature and velocity. The fast, cool, and dense downdrafts cary less enthalpy but more kinetic energy (per unit area) than the slow, hot upflows. The left panel of Fig. 13 shows the profiles of these two fluxes for the six simulated stars, all normalised to the respective total flux. Most models show profiles of the kinetic energy flux levelling off below the optical surface at values between ~20 to ~60% of the total flux, larger values corresponding to cooler stars. The F3V model, however, shows a monotonic increase of the kinetic energy flux for increasing depth. The right panel of Fig. 13 shows the radiative flux in the simulations. The transition from purely convective (F_rad/F_tot ≲ 0.1) to mainly radiative energy transport (F_rad/F_tot ≳ 0.9) is quite sharp (within one pressure scale height), except for the M-star simulations, for which this transition takes place over a more extended pressure range. For the K- and M-star simulations, the contribution of the radiative flux to the total flux is larger within the first few pressure scale heights below the surface. One pressure scale height below the surface (i.e. at p/p₀ = 2.72), in the M0V model, radiation carries already 17% of the total flux compared to 1.2% in the G2V simulation. Half a pressure scale height below the optical surface (p/p₀ = 1.65), radiation caries 60 to 70% of the flux in the K-star simulations, about 25 to 30% in the M-star models, but only 13.4% in the simulation of the G2V star. This means that the convective flux directly below the optical surface drops to a relatively small fraction of the total flux for the cooler simulations (especially the K-star simulations). As the convective flux is linked to vertical velocities and temperature contrast between up- and downflows, this is also reflected in the depth dependence of the rms of the vertical flow velocity (see Fig. 6) and the depth profile of the relative rms fluctuations of temperature (see top panel of Fig. 12). The peak of the vertical velocity rms is almost exactly at τ_R = 1 in the F- and G-star simulations but somewhat deeper in the simulations of cooler stars. Similarly, the peak of the relative temperature fluctuations shifts to higher optical depth from hotter to cooler models. Near the bottom of the simulation box, the contribution of the radiative flux eventually becomes negligible in all simulations except for F3V, where even at this depth radiation still carries about 0.5% of the energy flux. The much lower density and higher temperature in the subsurface layers of this star enable a somewhat more efficient radiative heat transport.

The inset in the right panel of Fig. 13 indicates the negative convective energy flux (which follows directly from ⟨F_rad ⟩ _z > F_tot) in the convective overshoot region, which is most prominent in the K-dwarf models.

Fig. 14 Comparison between the 3D MURaM simulations with 1D MLT models. The dashed curves represent results from a 1D MLT model atmosphere and the solid curves represent the averaged simulation results. Left: profile of temperature (3D: averaged on iso-p surfaces). Right: average vertical velocity, weighted by density (3D: averaged on iso-z surfaces; 1D: convective velocity). In the 1D models, the mixing-length parameter α was set to 1.5, 1.7, and 2.0 for F3V, G2V, and M0V, respectively.

Open with DEXTER

3.5. Comparison to 1D models

The horizontal averages of our 3D models can be compared to 1D mixing-length models. We used a model grid by Ludwig (priv. comm.) of 1D mixing-length calculations. Analogous to the definition of ⟨z⟩_τR = 1 ≡ 0 and ⟨p⟩_τR = 1 ≡ p₀ as reference points for the z- and p-scales in the averaged 3D stratifications, we chose z(τ_R = 1) ≡ 0 and p(τ_R = 1) ≡ p₀ as reference points for the 1D models.

Figure 14 shows depth profiles of temperature and rms of the vertical velocity as functions of pressure for three of our simulations (F3V, G2V, and M0V) and 1D models with the same g and T_eff. The mixing-length parameter was set to α = 1.5, 1.7, and 2.0 for F3V, G2V, and M0V, respectively. Although the temperature profiles in the almost adiabatic subsurface layers are influenced by the choice of α, they could not be brought into exact agreement with the profiles of the 3D results (which, moreover, depend on the averaging method). We therefore chose α close to the values by Trampedach & Stein (2011) for the mass mixing length, which are consistent with the literature values for the MLT-α cited there.

The left panel of Fig. 14 shows the run of temperature. The 3D results were averaged on surfaces of constant pressure (as a compromise between the iso-z and iso-τ_R averages, see Appendix A). The general shape of the curves does not differ strongly between averaged 3D and 1D results. However, the temperatures in the 1D models are in general somewhat lower than in the averaged 3D simulations, which can be partly explained by the way in which the 3D results were horizontally averaged but can also be an effect of the opacity in the optically thin layers, which differs between 1D and our 3D models (the 1D models were calculated with ATLAS6 opacities, see Kurucz 1979). The position and steepness of the strong photospheric gradient does not match the result of the more realistic 3D simulations. This has two reasons: First, this feature is very sensitive to the horizontal averaging method (see Appendix A) because of the strongly corrugated optical surfaces. Second, in this layer overshoot and the transition from convective to radiative energy transport play a major role. The physics behind these effects is essentially three-dimensional and has to be parameterised in a 1D model where only a very crude description of these effects is possible. In the lower part of the depth range considered, there is also a mismatch between the temperature gradients of the averaged 3D and the 1D results. A disparity in the superadiabaticity (particularly for F3V) and differences in the equation of state between the 3D and 1D models are responsible for this deviation.

In the right panel of Fig. 14, we show the run of the vertical velocity. For the 3D results, iso-z averages of the density-weighted vertical flow speed are shown. While the gradient of the subsurface velocity and the position of the velocity peak (at least for F3V and G2V) are similar between 1D and 3D models, the 1D models have lower velocities than the 3D simulations (by 10–30%) in the nearly adiabatic interior. The 1D models obviously lack any velocities in the convectively stable layers, where the simulations display overshooting flows.

In Fig. 15 the run of the superadiabaticity is shown. The 3D results are shown as iso-z and iso-τ_R horizontal averages. As discussed in Appendix A, near the photospheric transition, the iso-τ average is closer to the 1D description and therefore more useful for the comparison between 1D and averaged 3D models although the plain horizontal average is the physically more decisive quantity. Despite the necessary parameterisation of important physics in the 1D models, there is a qualitative similarity between 1D and 3D results. Below its peak, the superadiabaticity of the 1D models is higher by a factor of about 1.5 to 3 compared to the 3D models. For the F3V star, this deviation is the main cause for the difference of the temperature profiles, while for the simulations of the cooler stars, the differences in ∇_ad related to the equation of state between 1D and 3D models are larger than the small deviation in ∇−∇_ad. Around the superadiabatic peak, the profiles of 1D and 3D results differ more strongly. In this regime, which is more extended in terms of pressure scale heights for the M0V star than for the other two stars, the superadiabaticity in the 3D models is higher than predicted by the 1D models. The superadiabaticity in the atmospheric layers is qualitatively in agreement between 1D and 3D results.

Although there were small deviations in chemical abundances, equation of state, and opacities between 1D and 3D calculations, most of the differences in the upper part of the depth range shown can be attributed to the necessarily very crude treatment of convection and – most importantly – radiation in the 1D models versus the comprehensive 3D simulations.

Fig. 15 Comparison of the superadiabaticity of the 3D MURaM simulations with 1D MLT models. Top: 3D results averaged on planes of constant geometrical depth, z. Botom: 3D results averaged on iso-τR surfaces. The left subplots show the superadiabaticity of the atmospheric layers on a linear scale, the right subplots show the superadiabatic regime on a logarithmic scale. In the 1D models, the mixing-length parameter α was set to 1.5, 1.7, and 2.0 for F3V, G2V, and M0V, respectively.

Open with DEXTER

4. Concluding remarks

We analysed three-dimensional time-dependent hydrodynamical simulations of near-surface convection of six main-sequence stars of spectral types F3 to M2, including one with solar parameters. The coverage in depth (six to eight pressure scale heights below and above the optical surface) allows convective processes to be studied in the upper part of the convective envelope as well as the 3D structure of the atmospheres of the simulated stars.

In this first paper of the series, we present the averaged 3D stratification of the six simulations. Moreover, we analysed the energy transport through the upper part of the convective envelope.

Although the basic mechanisms driving and shaping the convective processes are similar in all models, the variation of stellar parameters (a factor of 5 in g and a factor of 2 in T_eff) entails differences in granule size and shape, intensity contrast, typical flow velocities, temperature fluctuations, scale heights etc. While some of these trends are more or less obvious effects of the variation of the stellar parameters, some are more subtle and are a consequence mainly of the highly non-linear temperature dependence of the opacity. Most notably, the low temperature sensitivity of opacity between 4000 and 5000 K, which is the main reason for “hidden” granulation in simulations of stars cooler than the Sun (especially of K stars), entails some differences in the atmospheric structure between the hotter two and the cooler four simulations of our sequence (see discussions of Figs. 6, 12, and 13).

The differences in the granulation necessarily have an impact on the formation of spectral lines, since line-of-sight velocity, pressure, and temperature (and thus also optical depth) are highly correlated with each other. The details of this impact will be analysed in the second paper of this series where we concentrate on the granulation, spectral line formation, and the 3D structure of the atmosphere.

Acknowledgments

The authors thank H.-G. Ludwig for providing 1D mixing-length models. They acknowledge research funding by the Deutsche Forschungsgemeinschaft (DFG) under the grant SFB 963/1, project A16. B.B. acknowledges financial support by the International Max Planck Research School (IMPRS) on Physical Processes in the Solar System and Beyond at the Universities of Braunschweig and Göttingen. A.R. acknowledges research funding from DFG grant RE 1664/9-1.

References

Online material

Appendix A: Horizontal averages

	Fig. A.1 Average profiles of the rms fluctuations of the geometrical height on iso-τR surfaces (left) and on iso-p surfaces (right) in units of the local pressure scale heights.
Open with DEXTER

Fig. A.2 Illustration of the different averages. Left panel: vertical cut through a part of the G2V simulation domain. Surfaces of constant optical depth (solid, red curves) and constant pressure (dashed, blue curves) are indicated; the underlying grey-scale image illustrate the density (normalised by the horizontal mean ⟨ϱ⟩z). Right panel: run of the horizontally averaged temperature vs. averaged pressure in the G2V simulation for the different averages ⟨ · ⟩ z (dotted, black), ⟨ · ⟩ τ (red, solid), and ⟨ · ⟩ p (blue, dashed); as a reference the position of some z and τ values is marked.

Open with DEXTER

We consider temporally and horizontally averaged quantities to study the mean stratification of the six simulated stars and compare them to 1D models. Depending on the context, the sensible “horizontal” average is an average over surfaces of constant geometrical depth, (Rosseland) optical depth, or (gas) pressure, denoted by ⟨ · ⟩ _z, ⟨ · ⟩ _τ, or ⟨ · ⟩ _p, respectively. The average ⟨ · ⟩ _τ is sensible for the layers around the optical surface and in the photosphere (especially for quantities involving the emergent intensity), while the average ⟨ · ⟩ _z is more relevant for the convective, deeper layers, where the stratifications are nearly adiabatic and largely independent from the radiation field. However, z is highly impractical as a coordinate scale for comparing simulations of stars of various spectral types, since the pressure scale height varies by more than a factor of ten between the six simulated stars presented in this paper (cf. Fig. 9). Therefore, we use the logarithm of the normalised pressure ⟨p⟩_z/p₀ as a more suitable depth coordinate (where p₀ = ⟨p⟩_τR = 1 is the average gas pressure at the optical surface) to illustrate depth dependences of quantities averaged on iso-z surfaces. This should not be confused with an average on surfaces of constant pressure, ⟨ · ⟩ _p.

The averages ⟨ · ⟩ _τ and ⟨ · ⟩ _p are defined as horizontal averages of a quantity on iso-τ and iso-p surfaces, respectively. For the average ⟨ · ⟩ _p, which is rarely used in this article, there is the problem that sometimes iso-p surfaces cannot be unambiguously defined since strong deviations from hydrostatic equilibrium in regions with Mach number of order unity can lead to a locally non-monotonic depth dependence of the gas pressure. Our iso-p surfaces are the deepest surfaces on which the pressure assumes the given value. This arbitrary choice does not significantly influence the results for the iso-p means in the simulations, except for the surface layers of the F3V simulation where the flows are mostly sonic and supersonic (cf. Fig. 6).

Figure A.1 shows the rms fluctuations of the geometrical depth on surfaces of constant optical depth (left panel) and pressure (right panel) as a measure of the corrugation of these surfaces. Due to this corrugation of the iso-p and iso-τ_R surfaces, profiles of quantities averaged in the different ways described above are not just distorted versions of each other, but can show a significantly different depth dependence of the same quantity. The differences between the three averaging methods are expected to be largest in the F- and G-type simulations, for which the corrugation of the iso-p and iso-τ_R surfaces is strongest.

The left panel of Fig. A.2 shows a vertical cut through some of the iso-p and iso-τ_R surfaces in the G2V simulation. Note that in the optically thin upper part of the simulation domain, the iso-p-surfaces follow the iso-τ_R surfaces. Although this is observed in all six simulations, this effect is most prominent in the G-type star (see discussion of Fig. 12 in Sect. 3.3). In the optically thick part, the temperature fluctuations determine the shape

of the iso-τ_R surfaces, since the opacity is highly temperature-sensitive in this regime. The iso-p surfaces, however, become almost flat planes in the deeper layers since the density contrast and deviations from hydrostatic equilibrium decrease with increasing depth.

The right panel of Fig. A.2 illustrates the different depth dependences of temperature, T, for the three different averages ⟨ · ⟩ _z, ⟨ · ⟩ _τ, and ⟨ · ⟩ _p. As for most of the figures in this paper, the gas pressure (here without normalisation) was used as depth coordinate. For ⟨T⟩_τ and ⟨T⟩_z, the pressure varies along the surfaces over which the average is performed. The depth coordinates are therefore averages themselves, namely ⟨p⟩_τ and ⟨p⟩_z, respectively, in these cases. As expected, the differences between the differently averaged temperature are largest at the optical surface and all three averages converge at large optical depth (log τ_R ≳ 4). As the deviations from hydrostatic equilibrium are small in the subsurface layers, ⟨ · ⟩ _p stays close to ⟨ · ⟩ _z. The average ⟨ · ⟩ _τ deviates more strongly near the photospheric transition since the big temperature fluctuations govern the opacity and thus lead to strongly corrugated iso-tau surfaces. In the atmosphere, the deviations between different averages of temperature become smaller with height. Especially, ⟨T⟩_p ≈ ⟨T⟩_τ, in these layers as expected, because the iso-τ surfaces roughly follow the iso-p surfaces.

If one aims at comparing 1D and averaged 3D results, one has to take into account that a 1D model does not have corrugated iso-τ surfaces. Profiles of quantities which change rapidly at the photospheric transition have a steep gradient in 1D models comparable to the local gradient in a 3D simulation. As the depth of the photospheric transition varies across the surface in a 3D simulation, these strong local gradients are smeared out and the similarity between 1D and averaged 3D results is obscured, if a plain horizontal average, ⟨ · ⟩ _z is used. The average ⟨ · ⟩ _τ is more appropriate in these cases for a comparison between 1D and averaged 3D profiles. This average however has no relevance below the photospheric transition, where ⟨ · ⟩ _z is more useful. For the stellar parameters used in our simulations, the ⟨ · ⟩ _p average seems a good compromise between the two other averaging methods for comparison with 1D models as it is converging towards ⟨ · ⟩ _τ in the atmosphere and towards ⟨ · ⟩ _z in the convection zone.

In order to obtain the temporally averaged profiles presented in this paper, first one of the horizontal averaging methods described above was applied to several snapshots with a time separation of 5–7 min, depending on the star. Then, for each quantity and at each depth point, the values of the different snapshots were averaged. The time dependence of the horizontal averages of most quantities under consideration (such as T, p, ϱ, etc.) was found to be very small, so that we found a small number of snapshots to be sufficient for a sensible temporal average. The mean profiles presented in this article are averages over six snapshots each.

All Tables

All Figures

	Fig. 1 Stellar parameters of the six models along with three isochrones by Bressan et al. (2012), solid line: zero-age main sequence (ZAMS), dashed line: age of 1 Ga, dotted line: age of 4.5 Ga (approximate solar age) on the log g–log Teff plane.
Open with DEXTER
In the text

	Fig. 2 Maps of bolometric intensity emerging vertically from single snapshots of the six simulations. The grey scale of each image is saturated at ±2 σI, where σI is the rms contrast (cf. Table 1). The significant difference in the length scales of the images is illustrated by the inset in the upper left panel, which shows all other images on the same scale as the one from the F3V simulation.
Open with DEXTER
In the text

	Fig. 3 Maps of the vertical velocity υz at constant geometrical depth z0 = ⟨z(τR = 1) ⟩ (average level of the optical surface) for four of the six models. Upward motions are blue, downward motions are red, colour scales saturate at ± 2·υz,rms(z0) (for values, see Table 1). Note that the horizontal scales are different (cf. Fig. 2).
Open with DEXTER
In the text

	Fig. 4 Same as Fig. 3, but for geometrical depth z2 with ⟨p(z2) ⟩ = 100 ⟨p(z0) ⟩, corresponding to 4.6 pressure scale heights below the optical surface.
Open with DEXTER
In the text

	Fig. 5 Properties of the convective flows. Left: ratio of rms values of horizontal and vertical flow velocity on surfaces of constant geometrical depth. Right: estimated horizontal scale of the vertical velocity patterns as derived from Eq. (1) on surfaces of constant geometrical depth. The horizontal scale is given in units of the horizontal box size, Xtot, for an easier comparison with Figs. 3 and 4. The solid curves refer to the four simulations shown in these figures, the dashed curves to the remaining two simulations.
Open with DEXTER
In the text

	Fig. 6 Flow velocity rms. Left: rms of the vertical component of the flow velocity on surfaces of constant geometrical depth. Right: rms of the modulus of the flow velocity in units of the local sound speed, cs, (Mach number) on surfaces of constant geometrical depth.
Open with DEXTER
In the text

	Fig. 7 Properties of up- and downflows. Left: relative area covered by upflows (υz > 0) on surfaces of constant geometrical depth as functions of normalised averaged pressure. Right: average speed of the upflows, υu: = ⟨υz|υz > 0 ⟩ z, (solid) and of the downflows, υd: = ⟨υz|υz < 0 ⟩ z, (dashed) as functions of normalised average pressure for four of the six simulations.
Open with DEXTER
In the text

	Fig. 8 Temperature and pressure stratifications on surfaces of constant optical depth. Top panel: temperature averaged over surfaces of constant optical depth. Middle panel: same as top panel, but normalised by the effective temperature of the respective model. Bottom panel: pressure averaged on surfaces of constant optical depth.
Open with DEXTER
In the text

	Fig. 9 Horizontally averaged local scale heights as functions of normalised pressure. Left: pressure scale height. Right: density scale height.
Open with DEXTER
In the text

	Fig. 10 Profiles of mean temperature ⟨T⟩z (left panel) and mean density ⟨ϱ⟩z (right panel) as functions of normalised pressure.
Open with DEXTER
In the text

	Fig. 11 Profiles of the logarithmic temperature gradient (left panel) and of superadiabaticity (right panel) averaged on iso-τR surfaces as function of pressure.
Open with DEXTER
In the text

	Fig. 12 Profiles of relative rms fluctuations. Top: temperature fluctuations on surfaces of constant optical depth. Middle: gas pressure fluctuations on surfaces of constant optical depth. Bottom: pressure fluctuations on surfaces of constant geometrical depth.
Open with DEXTER
In the text

	Fig. 13 Profiles of average energy flux, normalised to total flux. Left: convective energy flux split into net enthalpy flux and net kinetic energy flux. Right: radiative flux.
Open with DEXTER
In the text

	Fig. 14 Comparison between the 3D MURaM simulations with 1D MLT models. The dashed curves represent results from a 1D MLT model atmosphere and the solid curves represent the averaged simulation results. Left: profile of temperature (3D: averaged on iso-p surfaces). Right: average vertical velocity, weighted by density (3D: averaged on iso-z surfaces; 1D: convective velocity). In the 1D models, the mixing-length parameter α was set to 1.5, 1.7, and 2.0 for F3V, G2V, and M0V, respectively.
Open with DEXTER
In the text

	Fig. 15 Comparison of the superadiabaticity of the 3D MURaM simulations with 1D MLT models. Top: 3D results averaged on planes of constant geometrical depth, z. Botom: 3D results averaged on iso-τR surfaces. The left subplots show the superadiabaticity of the atmospheric layers on a linear scale, the right subplots show the superadiabatic regime on a logarithmic scale. In the 1D models, the mixing-length parameter α was set to 1.5, 1.7, and 2.0 for F3V, G2V, and M0V, respectively.
Open with DEXTER
In the text

	Fig. A.1 Average profiles of the rms fluctuations of the geometrical height on iso-τR surfaces (left) and on iso-p surfaces (right) in units of the local pressure scale heights.
Open with DEXTER
In the text

	Fig. A.2 Illustration of the different averages. Left panel: vertical cut through a part of the G2V simulation domain. Surfaces of constant optical depth (solid, red curves) and constant pressure (dashed, blue curves) are indicated; the underlying grey-scale image illustrate the density (normalised by the horizontal mean ⟨ϱ⟩z). Right panel: run of the horizontally averaged temperature vs. averaged pressure in the G2V simulation for the different averages ⟨ · ⟩ z (dotted, black), ⟨ · ⟩ τ (red, solid), and ⟨ · ⟩ p (blue, dashed); as a reference the position of some z and τ values is marked.
Open with DEXTER
In the text