Issue 
A&A
Volume 557, September 2013



Article Number  A26  
Number of page(s)  30  
Section  Stellar atmospheres  
DOI  https://doi.org/10.1051/00046361/201321274  
Published online  15 August 2013 
The Staggergrid: A grid of 3D stellar atmosphere models
I. Methods and general properties^{⋆,}^{⋆⋆}
^{1} MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85741 Garching, Germany
email: magic@mpagarching.mpg.de
^{2} Research School of Astronomy & Astrophysics, Cotter Road, Weston ACT 2611, Australia
^{3} StarPlan, Natural History Museum of Denmark/Niels Bohr Institute, Øster Voldgade 5–7, 1350 Copenhagen, Denmark
^{4} JILA, University of Colorado and National Institute of Standards and Technology, 440 UCB, Boulder, CO 80309, USA
^{5} Laboratoire Lagrange, UMR 7293, CNRS, Observatoire de la Côte d’Azur, Université de Nice SophiaAntipolis, Nice, France
^{6} Department of Physics & Astronomy, Michigan State University, East Lansing, MI 48824, USA
Received: 11 February 2013
Accepted: 10 July 2013
Aims. We present the Staggergrid, a comprehensive grid of timedependent, threedimensional (3D), hydrodynamic model atmospheres for latetype stars with realistic treatment of radiative transfer, covering a wide range in stellar parameters. This grid of 3D models is intended for various applications besides studies of stellar convection and atmospheres per se, including stellar parameter determination, stellar spectroscopy and abundance analysis, asteroseismology, calibration of stellar evolution models, interferometry, and extrasolar planet search. In this introductory paper, we describe the methods we applied for the computation of the grid and discuss the general properties of the 3D models as well as of their temporal and spatial averages (here denoted ⟨3D⟩ models).
Methods. All our models were generated with the Staggercode, using realistic input physics for the equation of state (EOS) and for continuous and line opacities. Our ~ 220 grid models range in effective temperature, T_{eff}, from 4000 to 7000 K in steps of 500 K, in surface gravity, log g, from 1.5 to 5.0 in steps of 0.5 dex, and metallicity, [Fe/H], from − 4.0 to + 0.5 in steps of 0.5 and 1.0 dex.
Results. We find a tight scaling relation between the vertical velocity and the surface entropy jump, which itself correlates with the constant entropy value of the adiabatic convection zone. The range in intensity contrast is enhanced at lower metallicity. The granule size correlates closely with the pressure scale height sampled at the depth of maximum velocity. We compare the ⟨3D⟩ models with currently widely applied onedimensional (1D) atmosphere models, as well as with theoretical 1D hydrostatic models generated with the same EOS and opacity tables as the 3D models, in order to isolate the effects of using selfconsistent and hydrodynamic modeling of convection, rather than the classical mixing length theory approach. For the first time, we are able to quantify systematically over a broad range of stellar parameters the uncertainties of 1D models arising from the simplified treatment of physics, in particular convective energy transport. In agreement with previous findings, we find that the differences can be rather significant, especially for metalpoor stars.
Key words: convection / hydrodynamics / radiative transfer / stars: abundances / stars: atmospheres / stars: fundamental parameters
Appendices A–C are available in electronic form at http://www.aanda.org
Full Table C.1 is available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/557/A26
© ESO, 2013
1. Introduction
The primary source of information for stellar objects is the light they emit, which carries information about the physical conditions at its origin. However, in order to interpret the information correctly, one first needs either theoretical or semiempirical models of the atmospheric layers at the surface of stars from where the stellar radiation escapes. Therefore, models of stellar atmospheres are essential for much of contemporary astronomy.
In the case of latetype stars, the theoretical modeling of stellar atmospheres is complicated by the presence of convective motions and turbulent flows as well as of magnetic fields in their envelopes (see review by Nordlund et al. 2009, and references therein). In particular, convection can significantly affect both the atmospheric stratification and emergent spectral energy distribution in these stars. Hence, in order to correctly represent the temperature stratifications in the outer layers of stars, from where the stellar light escapes, it is vital to accurately account for the interaction between radiative and convective energy transport at the optical surface.
The first realistic grids of lineblanketed atmosphere models for latetype stars appeared with the publication of MARCS (Gustafsson et al. 1975, 2008) and ATLAS models (Kurucz 1979; Castelli & Kurucz 2004). Subsequently, other onedimensional (1D) atmosphere codes, e.g. PHOENIX (Hauschildt et al. 1999) and MAFAGS (Grupp 2004), were developed to model the atmospheres of stars. In general, these theoretical 1D atmosphere models assume hydrostatic equilibrium, flux constancy, and local thermodynamic equilibrium (LTE). For the modeling of convective energy transport, they commonly employ the mixinglength theory (MLT, see BöhmVitense 1958), which is characterized by several free parameters, the most commonly known being the mixinglength l_{m}, or equivalently, the parameter α_{MLT} = l_{m}/H_{P}. Alternatively, some relatives thereof are available, such as the full turbulence spectrum (FTS) theory by Canuto & Mazzitelli (1991), which itself also has a free parameter. The values of these free parameters are not known from first principles and need to be calibrated based on observations or simulations. The mixinglength theory has in total four free parameters (see BöhmVitense 1958; Henyey et al. 1965; Mihalas 1970). These free parameters can be calibrated based on their effect on synthetic spectra, but usually only α_{MLT} is calibrated based on the reproduction of selected lines (Fuhrmann et al. 1993; Barklem et al. 2002; Smalley et al. 2002). Moreover, the free mixing length is calibrated in stellar evolutionary calculations by matching the observed luminosity and radius of the Sun at its current age (e.g. Magic et al. 2010). To construct simple yet realistic 1D models of convection is rather difficult, in particular convective overshooting beyond the classical Schwarzschild instability criterion is normally not considered in 1D atmospheric modeling. Attempts have been made at including its effects in 1D model atmospheres albeit with only limited success (Castelli et al. 1997).
The first numerical 1D model stellar atmosphere codes usually assumed a planeparallel geometry for the atmospheric stratification. This was later improved upon by changing to a spherical symmetry, leading to lower temperatures in the upper layers, in particular for giant stars, due to the dilution of the radiation field with increasing radial distance, which can cover a significant fraction of the stellar radius at low log g (see Gustafsson et al. 2008). Initially line blanketing was included by means of opacity distribution functions (ODFs, Gustafsson et al. 1975) with a few hundred ODFs covering the entire spectrum, eventually replaced by opacity sampling (OS) including thousands of wavelength points (Johnson & Krupp 1976). Nowadays, thousands of ODFs or hundreds of thousands of OS wavelengths are used. Despite such high resolution in wavelength, the computational costs for 1D atmosphere models are currently quite small, at least for LTE models. Large, homogeneous grids of atmospheres with up to ~ 10^{5} models exist (Gustafsson et al. 2008; Cassisi et al. 2004; Hauschildt et al. 1999), covering a wide range of stellar atmosphere parameters (T_{eff}, log g, and [Fe/H]).
Even though the 1D atmosphere models are based on numerous simplifications, they have demonstrated high predictive capabilities owing to major improvements in the atomic and molecular data (e.g. line lists by Kurucz 1993, or VALD by Piskunov et al. 1995). Also, the continuum opacity sources and the EOS have undergone similar developments. Thanks to these, 1D atmosphere models are in many respects very successful in comparisons with observations and are widely applied in astronomy today.
Another approach, almost exclusively used for solar atmosphere modeling, is the use of semiempirical models. In these models, the temperature stratification is inferred from observations (e.g. from lines forming at different heights or continuum centertolimb variations). Oftenused semiempirical 1D solar atmosphere models are the Holweger & Mueller (1974), VAL3C (Vernazza et al. 1976), Maltby et al. (1986) and MISS (Allende Prieto et al. 2001) models. A similar approach can be used to integrate spatially resolved observations and thus infer the threedimensional (3D) atmosphere structures using inversion techniques (Ruiz Cobo & del Toro Iniesta 1992; SocasNavarro 2011). Semiempirical modeling is rarely attempted for other stars, although exceptions exist (e.g. Allende Prieto et al. 2000).
Constructing more realistic models requires one to go beyond the 1D framework and model convection without relying on MLT. Stellar convection is an inherently 3D, timedependent, nonlocal, and turbulent phenomenon. Therefore, one cannot expect 1D models to reproduce all observed properties accurately, even with access to free parameters to tweak. The next natural step is to abandon some of these crude simplifications by constructing realistic 3D atmosphere models of solar convection. Early hydrodynamic simulations (Nordlund 1982; Nordlund & Dravins 1990; Steffen et al. 1989) revealed that stellar surface convection operates in a distinctly different fashion from the MLT picture. Instead of the homogeneous convective elements, they displayed highly asymmetrical motions with slow broad steady upflows interspersed with fast narrow turbulent downdrafts, sometimes even supersonic (e.g. Stein & Nordlund 1998, hereafter SN98; Asplund et al. 2000; Nordlund et al. 2009; Carlsson et al. 2004; Ludwig et al. 1999). The advent of 3D simulations, which are constructed from first principles, has enabled astronomers to predict various observables such as solar granulation properties and spectral line profiles astonishingly well. More recent solar 3D simulations are remarkably good at reproducing the observed centertolimb variation (e.g. Pereira et al. 2009a; Asplund et al. 2009; Ludwig 2006).
3D atmosphere models are by design free from the adjustable parameters of MLT and other parameters such as micro and macroturbulence that have hampered stellar spectroscopy for many decades. Instead, in 3D simulations, convection emerges naturally, by solving the timedependent hydrodynamic equations for mass, momentum and energyconservation, coupled with the 3D radiative transfer equation in order to account correctly for the interaction between the radiation field and the plasma. Also, the nonthermal macroscopic velocity fields associated with convective motions are rendered realistically, and various natural kinetic consequences such as overshooting and excitation of waves emerge from the simulations, without the need for further ad hoc modeling or additional free parameters. The inhomogeneities in the convective motions arise spontaneously and selforganize naturally to form a distinct flow pattern that exhibits the characteristic granulation at the surface. Furthermore, additional spectral observables such as limbdarkening and detailed spectral line shapes, including asymmetries and shifts, are also modeled unprecedentedly accurately with 3D models for the Sun (Nordlund et al. 2009; Pereira et al. 2009b).
For metalpoor latetype stars it has been shown (Asplund et al. 1999b; Collet et al. 2006, 2007) that the assumption of pure radiative equilibrium in the convectively stable photospheric layers of classical hydrostatic models is generally insufficient. In particular in the upper photosphere, the thermal balance is instead primarily regulated by radiative heating due to spectral line reabsorption of the continuumradiation from below and adiabatic cooling due to the expansion of upflowing gas. In metalpoor stars, the balance between heating by radiation and cooling by mechanical expansion of the gas occurs at lower temperatures because of the weakness and scarcity of spectral lines at low metallicities. By contrast, 1D MLT models have no velocity fields outside their convection zones, and are therefore in pure radiative equilibrium. The temperature stratification there is therefore regulated solely by radiative heating and cooling, thus neglecting altogether the adiabatic cooling component. This results in an overestimation of the temperatures by up to ~ 1000 K in 1D models at very low metallicities, which can potentially lead to severe systematic errors in abundance determinations based on 1D models (see Asplund et al. 1999b; Asplund & García Pérez 2001; Ludwig et al. 2010; Collet et al. 2009; González Hernández et al. 2010). These shortcomings of 1D models are manifested as inconsistencies in the analysis of observed spectra, such as abundance trends with excitation potential of the lines (e.g. analysis of NH lines in the very metalpoor star HE13272326, by Frebel et al. 2008) and discrepant abundances between atomic and molecular lines involving the same elements (e.g. Nissen et al. 2002). For further discussion, we refer to a review of possible impacts of 3D models on stellar abundance analysis by Asplund et al. (2005).
Additionally, there are discrepancies between observations and predictions from 1D models of the solar structure in the context of helioseismology, which point to mistakes in the outer layers of theoretical 1D stellarstructure models, and which are usually referred to as surface effects (Rosenthal et al. 1999). With classical 1D stellar structures, higher frequency pmodes of the Sun are systematically shifted due to discrepancies at the upper turning points of the modes, which occur in the superadiabatic peak at the top of the convection envelope. Rosenthal et al. (1999) found better agreement of stellar structures with helioseismic observations, when including the mean stratification of solar 3D models at the top, since the turbulent pressure, usually neglected in 1D models, extends the resonant cavity. Also, it was found that with 3D solar models the predicted pmode excitation rates are much closer to helioseismic observations (Nordlund & Stein 2001; Stein & Nordlund 2001). Ludwig et al. (2009b) compared the power spectra of the photometric microvariability induced by granulation and found good agreement between the theoretical predictions of 3D solar models and observations with SOHO.
With the aid of 3D simulations, stellar radii have been derived for a number of red giants from interferometric observations (Chiavassa et al. 2010, 2012). The determined stellar radii are slightly larger than estimated with the use of 1D models, which has an impact on the zero point of the effective temperature scale derived by interferometry. Furthermore, Chiavassa et al. (2012) showed that for interferometric techniques a detailed knowledge of the granulation pattern of planethosting stars is crucial for the detection and characterization of exoplanets.
Several 3D magnetohydrodynamics codes with realistic treatment of radiative transfer have been developed and applied to the modeling of stellar surface convection. Here, we make use of the Staggercode, which is developed specifically to run efficiently on the massively parallel machines available today (Nordlund & Galsgaard 1995 ^{1}; Kritsuk et al. 2011). The Bifrostcode is an Oslo derivative of the Staggercode (Gudiksen et al. 2011), tailored for simulations of the solar photosphere and chromosphere, and therefore including true scattering (Hayek et al. 2010). Other widely used codes are CO^{5}BOLD (Freytag et al. 2012), MURaM (Vögler et al. 2005) and ANTARES (Muthsam et al. 2010), which have been independently developed in the last decades. Beeck et al. (2012) compared solar models from three of the above 3D stellar atmosphere codes (Stagger, CO^{5}BOLD and MURaM), and showed that the models are overall very similar, despite the distinct numerical approaches. Most of the available 3D stellar convection codes are now highly parallelized, which when coupled with the computational power available today makes it feasible to construct grids of 3D convection simulations within a reasonable timescale. Grids of 2D and 3D atmosphere models already exist (Ludwig et al. 1999, 2009a; Trampedach 2007; Trampedach et al. 2013; Tanner et al. 2013). Clearly, the age of 3D atmosphere modeling has arrived, partly driven by the rising demand created by improved highresolution spectroscopic and asteroseismic observations.
In this paper, we present a new large grid of 3D model atmospheres for latetype stars, covering an extensive range of effective temperatures, surface gravities, and metallicities. In Sect. 2 we describe the methods that we followed in order to compute the 3D model atmospheres, with emphasis on the convection code we used (Sect. 2.1) and on the tools we developed for scaling the grid models and postprocess the results of the numerical calculations (Sect. 2.3). In Sect. 3, we present an overview of general properties (Sect. 3.1) of our simulations, and discuss the temporally and spatially averaged atmospheres (Sect. 3.2) from the 3D model atmospheres. We also compare our results with theoretical 1D models in Sect. 3.3 corresponding to the same stellar parameters ^{2}. These have been computed with a specifically newly developed 1D code that employs exactly the same EOS and opacities as the 3D simulations. Also, we compare our 3D model atmospheres with 1D models from grids widely adopted by the astronomical community (MARCS and ATLAS). Finally, in Sect. 4, we summarize our findings and outline a roadmap for our future ambitions on the many possible applications of the Staggergrid.
2. Methods
2.1. The STAGGERcode
The 3D model atmospheres presented here were constructed with a custom version of the Staggercode, a stateoftheart, multipurpose, radiativemagnetohydrodynamics (RMHD) code originally developed by Nordlund & Galsgaard (1995), and continuously improved over the years by its user community. In pure radiationhydrodynamics mode, the Staggercode solves the timedependent hydrodynamic equations for the conservation of mass (Eq. (1)), momentum (Eq. (2)), and energy (Eq. (3)) in a compressible flow coupled to the radiation field via the heating and cooling (per unit volume) term (4)which is computed from the solution of the radiative transfer equation (5)to account properly for the energy exchange between matter and radiation (here ρ denotes the density, v the velocity field, p the thermodynamic pressure, e the internal energy per unit volume ^{3}, g the gravity, the viscous stress tensor, q_{visc} = ∑ _{ij}τ_{ij}∂v_{i}/∂r_{j} the viscous dissipation rate, κ_{λ} the monochromatic opacity, I_{λ} the monochromatic intensity, J_{λ} = 1/4π∫_{Ω}I_{λ}dΩ the monochromatic mean intensity averaged over the entire solid angle, and S_{λ} the source function). We have ignored magnetic fields in the present grid of 3D convection simulations, however, we will study their effects in a future work.
2.1.1. Details on the numerics
The Staggercode uses a sixthorder explicit finitedifference scheme for numerical derivatives and the corresponding fifthorder interpolation scheme. The solution of the hydrodynamic equations is advanced in time using an explicit thirdorder RungeKutta integration method (Williamson 1980). The code operates on a staggered, Eulerian, rectangular mesh: the thermodynamic variables, density and internal energy per volume, are cellcentered, while momentum components are defined at cell faces. Also, in the MHD mode, the components of the magnetic field B (electric field E) are defined at the cell faces (edges). This configuration allows for a fluxconservative formulation of the magnetohydrodynamic equations, at the same time ensuring that the magnetic field remains divergencefree. The solution of the discretized equations is stabilized by hyperviscosity which aims at minimizing the impact of numerical diffusion on the simulated flow, while providing the necessary diffusion for largeeddy simulations with finitedifference schemes (see also Nordlund & Galsgaard 1995, for further details). The values of the numerical viscosity parameters ^{4} are empirically tuned for the solar surfaceconvection simulation: they are set large enough to stabilize the numerical solution of the hydrodynamic equations and, at the same time, kept small enough to reduce their smoothing of the flow’s structures. The same optimized values of the parameters are then applied to all other simulations in the grid.
The version of the Staggercode we used for this work is fully MPIparallel. The parallelization scales well with the number of cores. For this project, the simulations were typically run on 64 cores.
2.1.2. Geometrical properties
The setup of the simulations is of the socalled boxinastar type: the domain of the simulations is limited to a small representative volume located around the stellar photosphere and including the top portion of the stellar convective envelope. The boundary conditions of the simulation box are periodic in the horizontal directions and open vertically. Gravity is assumed to be constant over the whole extent of the box, neglecting sphericity effects. However, since the size of the simulation domains correspond to only a fraction of the total radii of the stars (0.4% and ~ 10% of the stellar radius for the solar simulation and for a typical log g = 1.5 red giant simulation, respectively) such effects can be regarded as small for the purposes of the current grid of models. Also, for simplicity, the effects of stellar rotation and associated Coriolis forces are neglected in the present simulation setup, as it would add two more dimensions to the grid.
At the bottom, the inflowing material has a constant value of specific entropy per unit mass, which ultimately determines the emerging effective temperature. While the domains of our simulations cover only a small fraction of the convective zone, the boxinastar setup is still valid because the bulk upflows at the bottom boundary of the simulations carry essentially the same entropy value as in deeper layers and are mostly unaffected by entrainment with cooler downflows. At the beginning of each simulation, the entropy of the inflowing gas at the bottom is adjusted in order to yield the desired T_{eff} and, after that, is kept unchanged during the entire run (see Sect. 2.3). Furthermore, pressure is assumed to be constant over the whole bottom layer.
The physical dimensions in the horizontal directions are chosen to be large enough to cover an area corresponding to about ten granular cells. The vertical dimensions are extended enough for the simulations to cover at least the range of − 5.0 < log τ_{Ross} < + 6.0 in terms of Rosseland optical depth (in fact they range on average from − 7.3 < log τ_{Ross} < + 7.5), which typically corresponds to approximately six orders of magnitude in pressure (about 14 pressure scale heights). All of the simulations have a mesh resolution of 240^{3}, since a resolution of about 200^{3} − 250^{3} was found to be adequate (see Asplund et al. 2000). Five layers at the bottom and the top are reserved for the socalled ghostzones: these extra layers serve to enforce boundary conditions for the highorder derivatives in the vertical direction. The spacing between cells in the horizontal direction (Δx,Δy) is constant, ranging from about 6 km in dwarfs to about 25 Mm in giants, while it varies smoothly with depth in the vertical direction, in order to resolve the steep temperature gradients near the optical surface. These are the layers from where the continuum radiation escapes; they are characterized by a sharp transition between stellar interior and outer layers in terms of thermodynamic quantities such as temperature, internal energy, and entropy that marks the beginning of the photosphere. Also, the steepest temperature gradients are found in the superadiabatic region just below the optical surface (0.0 < log τ_{Ross} < 2.0). Therefore, it is very important that the thin transition layer around the optical surface is wellresolved in order to ensure an accurate modeling of the radiative transfer and to avoid spurious numerical artifacts from insufficient spatial resolution.
2.1.3. Equation of state
We use the realistic equation of state (EOS) by Mihalas et al. (1988), which explicitly treats excitation to all bound states of all ionization stages, of all included elements. We have custom computed tables for a mix of the 17 most abundant elements (H,He,C,N,O,Ne,Na,Mg,Al,Si,S,Ar,K,Ca,Cr,Fe and Ni). The only molecules that are included in the EOS are H_{2} and , and they are treated on equal footing with the atoms and ions. For the solar abundances we employed the latest chemical composition by Asplund et al. (2009), which is based on a solar simulation performed with the same code and atomic physics as presented here. Our choice for the EOS, is supported by Di Mauro et al. (2002) who showed that solar models based on the EOS by Mihalas et al. (1988) show better agreement with helioseismology in the outer 20 Mm (≥0.97 R_{⊙}), compared to models based on the OPALEOS. We inverted the Mihalas et al. (1988) EOS tables, hence the temperatures and the thermodynamic pressures are tabulated as a function of density and internal energy. This inversion exploits the analytical derivatives provided in the EOS tables to minimize losses in accuracy. These analytical derivatives are also used in the bicubic spline interpolation in the inverted tables.
2.1.4. Opacity
We use the continuum absorption and scattering coefficients listed in detail and with references by Hayek et al. (2010). These include the sophisticated calculations by Nahar (2004)^{5} for the first three ions of all metals we include, except for K and Cr. These calculations are improvements over those forming the basis for the OP opacities ^{6} (Badnell et al. 2005). The line opacity is supplied by the opacity sampling (OS) data that was also used for the newest MARCS grid of stellar atmospheres Gustafsson et al. (2008), which are in turn based on the VALD2 database ^{7} (Stempels et al. 2001) of atomic and molecular lines.
2.1.5. Radiative transfer
The radiative heating and cooling rate (Eq. (4)) is evaluated by solving the radiative transfer equation (6)where τ_{λ} = ∫ρκ_{λ}ds denotes the monochromatic optical depth along a given direction s, with a method similar to that by Feautrier (1964). The equation is solved at each time step and grid point on long characteristics, along the vertical direction and along eight additional inclined angles (two μ = cosθ and four ϕangles) by tilting the (domaindecomposed) 3D cube. Given the opacity κ_{λ} and the source function S_{λ}, the monochromatic intensity I_{λ} can be obtained by solving Eq. (6) and the radiative heating and cooling rate computed by integrating ρκ_{λ}(I_{λ} − S_{λ}) over solid angle and wavelength. We use the Radau quadrature to determine the optimal ray directions to approximate the angular integral in the calculation of the radiative heating and cooling rate as a weighted sum. For the radiative transfer calculations, we employ opacities as described above (Sect. 2.1.4).
Computing the full monochromatic solution to the radiative transfer equation in 3D at each time step is extremely expensive. The cost of the radiative transfer calculations however can be reduced enormously by opting instead for an approximated solution based on the opacity binning or multigroup method (Nordlund 1982; Skartlien 2000). Following this method, we sort all sampled wavelength points into different bins based on the spectral range they belong to and on their associated opacity strength or, better, their formation depth, i.e. the Rosseland optical depth , where the monochromatic optical depth equals unity. In this way, wavelength points characterized by similar formation heights and belonging to the same spectral interval are grouped together (see Fig. 3). For each simulation, we use the 1D temporal and spatial mean stratification to estimate the formation heights of the various wavelengths and sort the wavelengths into the different opacity bins. The bin selection and wavelength sorting process is performed twice during the simulation’s relaxation phase after updating the individual mean stratifications, but is kept unchanged during the production runs that make up the timeseries presented in this work.
To each bin, we assign a mean opacity κ_{i} which accounts for the contribution from both continuum and line opacities. To compute the mean opacities, we differentiate between diffusion and freestreaming limits, i.e. between the optical thick and optical thin regimes, below and above the photospheric transition zone, respectively. The mean binopacity κ_{i} is calculated as a Rosselandlike average (7)in the optical thick regime, and as a meanintensityweighted mean opacity (8)in the optical thin regime, where is the set of wavelength points assigned to bin i. For bin i, the transition from one regime to the other around that bins optical surface is achieved by means of an exponential bridging of the two averages: (9)All simulations presented here have been run with the radiative transfer in the strict LTE approximation, i.e. under the assumption that he monochromatic source function S_{λ} (in Eq. (4)) is the Planck function at the local gas temperature, i.e. . For each bin i, we compute an integrated source function by summing up the contributions from all wavelength points in the bin; (10)Collet et al. (2011) showed that, with this opacity binning implementation, the approximation of strict LTE results in a temperature stratification is very similar to the case, where scattering is properly treated, as long as the contribution of scattering to the extinction is excluded when averaging the mean opacities κ_{J,i} (Eq. (8)) in the optically thin layers (“free streamingregime”), but include it as true absorption when averaging the mean opacities κ_{Ross,i} (Eq. (7)) in the optically thick layers (“diffusion approximation regime”). They also showed that including scattering as true absorption leads to erroneous atmosphere structures due to overestimated radiative heating in the optically thin layers. However, these findings have so far being verified for only a small sample of stellar parameters, therefore we cannot rule out that scattering needs to be accounted for properly in certain cases. Nonetheless, evaluating the radiative transfer in strict LTE greatly eases the computational burden compared to the case, where the contribution of scattering is included to the total extinction (Hayek et al. 2010).
The radiative transfer equation is solved for the individual opacity bins (Eq. (5)) for all layers that have , while in the deeper layers, we use instead the diffusion approximation, which is fulfilled to a high degree at such depths. With the opacity binning approximation, the radiative heating rate term (Eq. (4)) takes then the form (11)where J_{i} is the mean intensity computed from the solution of the radiative transfer equation for bin i.
For the relaxation phase of the simulation runs we considered six bins, while for the final models we used twelve opacity bins. We have developed an algorithm for the bin selection, which will be explained further below (see Sect. 2.3.2). Towards lower surface gravities (log g ≲ 2.0) and higher effective temperatures, numerical artifacts in the radiative transfer can occasionally develop and manifest as a Moiré pattern in the integrated outgoing intensities due to very steep temperature gradients in the photosphere. For those situations, we solve the radiative transfer equation on an adaptive mesh with finer vertical resolution, which is dynamically optimized to resolve regions where temperature gradients are steeper. The radiative heating and cooling rates computed on the adaptive mesh are then interpolated back to the coarser hydrodynamic depth scale under the consideration of energy conservation.
2.2. The STAGGERgrid models
The Staggergrid covers a broad range in stellar parameters with 217 models in total. The range in effective temperature is from T_{eff} = 4000 K to 7000 K in steps of 500 K, while the gravity ranges from log g = 1.5 to 5.0 in steps of 0.5. The grid also covers a broad range in metallicity starting from [Fe/H] = −4.0 to + 0.5 in steps of 1.0 below − 1.0, and steps of 0.5 above that ^{8}. We decided to apply the same parameters T_{eff} and log g for all metallicities, in order to facilitate the interpolation of (averaged) models within a regular grid in stellar parameters. In addition, the grid also includes the Sun with its nonsolar metallicity analogs, and four additional standard stars, namely HD 84937, HD 140283, HD 122563 and G 6412 that are presented in Bergemann et al. (2012). For metalpoor chemical compositions with [Fe/H] ≤ − 1.0 we applied an αenhancement of [α/Fe] = + 0.4 dex, in order to account for the enrichment by corecollapse supernovae (Ruchti et al. 2010).
In Fig. 1, we present an overview of our simulations in stellar parameter space. Therein, we also show evolutionary tracks (Weiss & Schlattl 2008) for stars with masses from 0.7 to 1.5 M_{⊙} and solar metallicity, in order to justify our choice of targeted stellar parameters. Hence, the grid covers the evolutionary phases from the mainsequence (MS) over the turnoff (TO) up to the redgiant branch (RGB) for lowmass stars. In addition, the RGB part of the diagram in practice also covers stars with higher masses, since these are characterized by similar stellar atmospheric parameters.
2.3. Scaling and relaxing 3D models
Generating large numbers of 1D atmosphere models is relatively cheap in terms of computational costs, but the same is not true for 3D models. Based on our experiences from previous simulations of individual stars, we designed a standard workflow of procedures for generating our grid. More specifically, we developed a large set of IDLtools incorporating the various necessary steps for generating new 3D models, which we then applied equally to all simulations. The steps are:

Scale the starting model from an existing, relaxed 3Dsimulation, and perform an initial run with six opacity bins, sothat the model can adjust to the new stellar parameters.

Check the temporal variation of T _{eff} and estimate the number of convective cells. If necessary, adjust the horizontal sizes, in order to ensure that the simulation box is large enough to enclose at least ten granules.

If the optical surface has shifted upwards during the relaxation, add new layers at the top of it to ensure that ⟨log τ_{Ross}⟩_{top} < −6.0.

Determine the period π _{0} of the radial pmode with the largest amplitude, then damp these modes with an artificial exponentialfriction term with period π _{0} in the momentum equation (Eq. (2)).

Let the natural oscillation mode of the simulation emerge again by decreasing the damping stepwise before switching it off completely.

Recompute the opacity tables with 12 bins for the relaxed simulation.

Evolve the simulations for at least ~7 periods of the fundamental pmode, roughly corresponding to ~ 2 convective turnover times, typically, a few thousand timesteps, of which 100–150 snapshots equally spaced were stored and used for analysis.
During these steps the main quantities of interest are the time evolution of effective temperature, pmode oscillations, and drifts in the values of the mean energy per unit mass and of the mean density at the bottom boundary, which indicate the level of relaxation. When the drifts in these above properties stop, we regard the simulation as relaxed. If these conditions were not fulfilled, we continued running the model, to give the simulation more time to properly adjust towards its new quasistationary equilibrium state. Also, when the resulting effective temperature of an otherwise relaxed simulation deviated more than 100 K from the targeted T_{eff}, we rescaled the simulation to the targeted value of T_{eff} and started over from the top of our list of relaxation steps.
Fig. 1 Kiel diagram (T_{eff} − log g diagram) showing the targeted Staggergrid parameters for the 217 models, comprising seven different metallicities (colored circles). Four additional standard stars (see text) are also indicated (squares). In the background, the evolutionary tracks for stellar masses from 0.7 to 1.5 M_{⊙} and for solar metallicity are shown (thin grey lines). 
The interplay between EOS, opacities, radiative transfer and convection can shift the new location of the photosphere, when the initial guess made by our scaling procedure slightly misses it. This is the case for a few red giant models leading to upwardsshifts of the optical surface and of the entire upper atmosphere during the adjustment phase after the scaling, with the average Rosseland optical depth ending up to be larger than required, i.e. ⟨log τ_{Ross}⟩_{top} ≥ − 6.0. In order to rectify this, we extended those simulations at the top by adding extra layers on the top, until the top layers fulfilled our requirements of ⟨log τ_{Ross}⟩_{top} < − 6.0.
2.3.1. Scaling the initial models
To start a new simulation, we scale an existing one with parameters close to the targeted ones, preferably proceeding along lines of constant entropy of the inflowing gas at the bottom in stellar parameter space (see Fig. 6). In this way, we find that the relaxation process is much faster. In order to generate an initial model for a set of targeted parameters, we scale temperature, density, and pressure with depthdependent scaling ratios derived from two 1D models, with parameters corresponding to the current and intended 3D model (Ludwig et al. 2009a). For this, we used specifically computed 1D envelope models (MARCS or our own 1D models, see Sect. 3.3.1), which extend to log τ_{Ross} > 4.0. The reference depthscale for all models in the scaling process is the Rosseland optical depth above the photosphere and gas pressure normalized to the gas pressure at the optical surface below it (log τ_{Ross} > 0.0).
Fig. 2 Top panel: nonequidistant vertical spacing Δz of the depth scale as a function of geometrical depth in our solar model (solid line). The zscale is optimized to resolve the flows and thermal structure, which naturally results in the highest spatial resolution around the photosphere. Furthermore, we also show the maximum of the vertical gradient of the absorption coefficient max⟨dlnα_{Ross}/dz⟩ as a function of depth (dotted line). Bottom panel: aspect ratio Δz/Δx (solid line) and we also indicated its lower allowed limit with 1:4 (dashed line). The actual verticaltohorizontal aspect ratio ranges from 0.26 at the photosphere to 1.18 at the bottom of the simulation domain. 
After the initial scaling, we construct the geometrical depth scale z for the new simulation by enforcing the same (quasi)hydrostaticequilibrium condition as in the starting simulation, but with the newly scaled pressure and density. The resulting new zscale is usually not smooth, therefore we generate a new zscale, which is optimized to resolve the region with the steepest (temperature) gradients, as shown in Fig. 2. The density, energy, and velocity cubes are then interpolated to this new geometrical depth scale. The new zscale is constructed using the variation with depth of the (smoothed) maximum of the derivative of the Rosseland absorption coefficient, max⟨dlnα_{Ross}/dz⟩, as a guide. The basic idea behind this approach is to vertically distribute the mesh points as evenly as possible on the opticaldepth scale. With such an optimized zscale we can efficiently resolve the same features with fewer gridpoints, compared to an equidistant vertical mesh. Furthermore, a limiting verticaltohorizontal aspect ratio (Δz/Δx and Δz/Δy) of 1:4 over the whole vertical extent is enforced. We find that this value represents in practice an optimal lower limit to the aspect ratio, with respect to numerical stability and accuracy of the solution of the radiative transfer equation along inclined rays. Finally, the position of the zeropoint in the depth scale is adjusted to coincide with the position of the mean optical surface, i.e. ⟨τ_{Ross}⟩(z = 0) = 1.
At fixed surface gravity and metallicity, the mean diameter of the granules, which is used for determining the horizontal extent of the simulation, increases with higher effective temperature (see Figs. 11 and Sect. 3.1.5). The number of granules present in the simulation box is retrieved with the aid of the contour routine in IDL. Based on the map of the temperature below the surface (the vertical velocity would serve equally well), a contour chart of the significantly hotter granules is extracted, from which the number of granules is counted. Concerning the temporal resolution of the simulation sequences of the final production runs, the frequency, at which snapshots are stored, is based on the soundcrossing time of one pressure scale height, H_{P}, in the photosphere, i.e. (12)(see Cols. 14 and 15 in Table C.1). With the help of functional fits of the dependence of granule sizes and soundcrossing time scales on stellar parameters, the horizontal sizes of the simulation boxes and the snapshot sampling times can be estimated rather accurately in advance (see Appendix B).
2.3.2. Selection of the opacity bins
As we mentioned earlier, in Sect. 2.1.5, the purpose of the opacitybinning approximation is to reproduce the radiative heating and cooling rates as accurately as possible with a small number of opacitybins, in order to reduce the computational burden. For the assignment of wavelength points to bins, we first compute the opacity strengths for all of the ≳10^{5} wavelength points in the opacitysampling (OS) data from the MARCS package (Gustafsson et al. 2008). The histogram of their distribution as a function of wavelength (see Fig. 3) exhibits a characteristic “L”shape. Shorter wavelengths (UV) require more bins to resolve the wide range in opacity strength, while the lower part of the Lshaped distribution at longer wavelengths (optical and IR) calls for a better resolution in terms of wavelength. Therefore, we initially make a division in wavelength at λ_{X}, between the UV and the optical/IR (see boundaries of bin 1, 11 and 12 in Fig. 3) and comprising approximately an equal number of wavelength points. These two regions are then in turn subdivided evenly into opacity bins according to the number of λpoints. By trial and error, we found that a binning scheme with three bins in the λ < λ_{X} region and eight bins for λ > λ_{X}, one of which being a large one and comprising the stronger lines in the optical and IR (bin number 10 in Fig. 3) gives a good representation of the monochromatic radiative heating and cooling. We iterate the bin selection with slight differences (e.g., one additional division in opacity strength for the 8 bins in the lower part of the optical and IR) and by small adjustments, and choose the bin selection with the smallest relative difference between the total heating rates computed with opacity binning q_{bin} and the full monochromatic solution q_{λ} for the average stratification of the 3D simulation, i.e. (13)We found that the individual selection of some of the bins displays a highly nonlinear response to small changes. In most cases an even distribution was favored by the minimization. Naturally, our method will typically find only a local minimum due to the small sample of iterations instead of a true global minimum. However, our method is a fast, repeatable, and automatic selection of the opacity bins, which minimizes the human effort significantly, while at the same time yielding very satisfactory results. Moreover, the possible deviation from the global minimum due to our automated bin selection and its resulting uncertainties are anyways smaller than the overall uncertainties associated with the opacity binning method.
Fig. 3 Twelve opacity bins selected for the solar simulation by plotting the opacity strength (or, more precisely, the formation height) against wavelength for all sampled wavelength points. The individual bin elements are indicated by colored symbols. For clarity, we plotted only a subset of the wavelength points considered for the opacity binning procedure. In the background, we included the smoothed histogram of the opacity strength distribution (blue contour). This shows how the majority of λpoints are mostly concentrated close to the continuumforming layers and only a smaller fraction contributes to lines. 
In Fig. 4, we compare the resulting radiative heating and cooling rates from the monochromatic calculation against those from the opacity binning solution for the mean stratification of our solar model. The radiative heating and cooling rates from the simplified opacity binning appear rather similar to those from the monochromatic solution, thereby supporting our approach. For the solar model, our algorithm finds a bin selection that is just slightly less accurate (max[δq_{bin}] = 2.78%) than an optimized manual bin selection (1.86%). Incidentally, with six bins, we get max[δq_{bin}] = 3.54%. We obtain an average max[δq_{bin}] for all the grid models of , while with six bins we get . We find that max[δq_{bin}] increases slightly with T_{eff} and [Fe/H]. We note that the opacity binning method with its small number of bins states an approximation for the radiative transfer, therefore, despite the small values for max[δq_{bin}] further improvement is necessary.
3. Results
The spatially and temporally averaged mean 3D stratifications (hereafter ⟨3D⟩) from all of our 3D models will be available online. The methods we applied to average our models are explicitly described in a separate paper. We provide the models in our own format, but also in various commonly used formats suited for standard 1D spectrum synthesis codes such as MOOG (Sneden 1973), SYNTHE (Kurucz 1993) and Turbospectrum (Plez 2008; de Laverny et al. 2012), together with routines to interpolate the ⟨3D⟩ models to arbitrary stellar parameters. In this paper the discussion will be confined to global properties and mean stratifications only. More extensive discussions and presentations of the wealth of details present in the data of our 3D RHD models will be performed systematically in subsequent papers.
Fig. 4 Comparison of the radiative heating and cooling resulting from monochromatic computations q_{λ} (filled dots) and the opacity binning method q_{bin} (solid line) for the solar model mean stratification. In the top panel we show both q_{rad} vs. optical depth, while in the bottom panel, we compare the two against each other. 
3.1. Global properties
In Table C.1, we have listed the stellar parameters together with the thermodynamic values fixed for the inflows at the bottom, i.e. the internal energy ε_{bot}, density ρ_{bot} and entropy s_{bot}, as well as important global properties for our 3D simulations. Before we consider the ⟨3D⟩ stratifications in Sect. 3.2, we briefly discuss some (temporally averaged) global properties.
3.1.1. Stellar parameters
Surface gravity and metallicity are input parameters for a simulation, while the effective temperature is a property ensuing from the fixed entropy of the inflowing material at the bottom s_{bot}. We calculate the effective temperature from the spatially averaged emergent radiative energy flux F_{rad} and the StefanBoltzmann law T_{eff} = [F_{rad}/σ]^{1/4}, with σ being the StefanBoltzmann constant. In Column 1 of Table C.1 we have listed the resulting temporally averaged T_{eff} of our final, relaxed simulations. These differ somewhat from the targeted T_{eff}s, since we do not know a priori, the relation between s_{bot} and T_{eff}. However, the majority of our models (72%) deviate less than 50K, and the mean deviation for the whole grid is .
3.1.2. Constant entropy of the adiabatic convection zone
Fig. 5 Overview of the constant entropy value of the adiabatic convection zone, which is given by the fixed entropy of the inflowing plasma at the bottom, s_{bot}, (circles) as well as the atmospheric entropy minimum, s_{min}, (squares) for two metallicities ([Fe/H] = −0.5 and 0.0, blue and black respectively) against T_{eff}. The jump in entropy, Δs, is indicated with vertical dotted lines. Simulations with the same gravity are connected with functional fits for s_{bot} and s_{min} (solid and dashed lines respectively; see Appendix B), while similar simulations with different [Fe/H] are connected with short solid black lines. 
The main input parameter that has to be adjusted is s_{bot}, which has the same value as the entropy in the deep convection zone due to the adiabaticity of convection, i.e. s_{bot} = s_{ad} (Steffen 1993). This is also the reason, why the results from our rather shallow boxes are valid. We set s_{bot} by specifying a fixed value for the density and energy per unit mass for the inflowing material at the bottom, ρ_{bot} and ε_{bot}. The actual values of ε_{bot}, ρ_{bot} and s_{bot} applied in our simulations are given in Table C.1. Furthermore, we provide functional fits for s_{bot} (see Appendix B). We compute the entropy by integrating ^{9} the first law of thermodynamics in the form (14)adding an arbitrary integration constant in order to shift the zeropoint of the entropy to a similar value as in Ludwig et al. (1999). In Fig. 5, we show s_{bot} against T_{eff} for [Fe/H] = 0.0 and − 0.5, as an example. The value for s_{bot} increases exponentially with higher T_{eff} and with lower log g, and linearly with metallicity with a moderate slope.
In order to increase the effective temperature solely, i.e. the emergent radiative flux F_{rad} at the top boundary, the total energy flux ascending from the convection zone has to increase by that same amount due to conservation of energy (see Sect. 3.2.8). On the other hand, when we keep T_{eff} fixed and decrease the surface gravity, this in turn will cause the density to decrease correspondingly (see Sect. 3.2.4). Therefore, to maintain the same energy flux, either the transported heat content (Δs, ε) or the mass flux (ρ or v_{z}) needs to be enhanced. When the energy flux is carried by a larger mass flux, then we speak of an enhancement in convective efficiency ^{10} (see Sects. 3.2.4, 3.2.2 and 3.2.8). When one considers ε_{bot} with stellar parameters, then it clearly depicts qualitatively the same characteristic changes as s_{bot}. By inserting the perfect gas law ^{11} in Eq. (14) one obtains (15)from which one can immediately see that the entropy increases with internal energy, s ∝ ε, and also increases with lower density s ∝ − lnρ.
On the other hand, an increase in metallicity leads to a higher entropy of the adiabat and also a larger atmospheric entropyjump (see Fig. 5). Furthermore, we find increased velocities (and Δs) and decreased densities at higher [Fe/H] (see Sects. 3.2.2 and 3.2.4), which in turn affects the convective efficiency. The dependence on metallicity can be unveiled with the following approximation. The opacity (and absorption coefficient) increases with higher [Fe/H], since the opacity depends directly on the metallicity. The hydrostatic equilibrium can be written in terms of optical depth as dp_{th}/dτ_{Ross} = g/κ_{Ross}. From the EOS (also ideal gas law), one can see that the pressure scales with the density, p_{th} ∝ ρ. Therefore, when one would fix T_{eff} and log g, but increase metallicity (and opacity), then the hydrostatic balance will be realized at a lower density stratification (see bottom and middle panel in 18), which is also given by 1D MLT models. The lower density stratification will result in higher s_{bot} and Δs (Eq. (15) and top panel in18).
We emphasize that the dependence of both s_{bot} and Δs with stellar parameters is quite nontrivial, since not only is it coupled to the changes in the total energy fluxes, but it is also affected by the differences in the transition from convective to radiative transport of energy with stellar parameters. In particular, the nonlocal radiative transfer depends nonlinearly on the conditions present in stellar atmospheres, especially changes in the opacity and the EOS will strongly influence the radiative transfer. Additionally, s_{bot} will be influenced by changes in the efficiency of the convective energy transport, that is in the convective fluxes arising from the hydrodynamics (see Sect. 3.2.8 for more details).
Fig. 6 Lines of constant entropy of the adiabat s_{bot} from 1.3 to 4.0 × 10^{13}erg/g/K in steps of 0.3 × 10^{13}erg/g/K and for [Fe/H] = −3.0 and 0.0 (blue solid and black dashed lines respectively) in the Kiel diagram. Additionally, we show evolutionary tracks for 0.7 to 1.5 M_{⊙} with solar metallicity (thin grey lines). 
An analytical derivation of s_{bot} as a function of stellar parameters is rather difficult, as explained above. Nonetheless, we can fit s_{bot} with stellar parameter in a functional form based on the results from our simulations (see Appendix B). This has been done previously, based on 2D RHD models with solar metallicity by Ludwig et al. (1999). In Fig. 6, we show how s_{bot} varies across the Kiel diagram (T_{eff} − log g diagram) for [Fe/H] = 0.0 and −3.0. In the case of s_{bot}, our results with solar metallicity are qualitatively in good agreement with those of Ludwig et al. (1999), despite the inherent differences between 2D and 3D convection simulations, the adopted EOS, opacities, and radiative transfer treatment. We find that s_{bot} (which depicts s_{ad}) lies on lines of constant entropy in the Kiel diagram, in fact for the solar metallicity these lines with the same entropy of the adiabat run almost diagonally. Moreover, towards lower metallicity, we find two significant differences for the lines of constant s_{bot}, the first one being that the slopes of the lines steepen, and the second being that the distances in the T_{eff} − log g plane between the lines decrease. The latter implies metalpoor stars feature a broader range in s_{bot} compared to metalrich ones. As we will see in Sect. 3.2, this has important consequences for the stratifications.
3.1.3. Entropy jump
The upflows enter the simulation box at the bottom with the constant entropy value of the adiabatic convection zone, s_{bot}, and ascend until they reach the superadiabatic region (SAR) just below the visible surface, where the convective energy is converted to radiative energy. In the photosphere, the mean free path for the continuum radiation grows large enough for the gas to become transparent, and the overturning upflow at the surface loses its internal energy as photons escape and carry away entropy. Further above in the nearly isothermal atmosphere (with constant ε, see Fig. 16) with an exponentially decreasing density the entropy increases again due to the EOS (see Eq. (15)). This leads to a minimum, s_{min} = min[⟨s⟩], just above the surface (log τ_{Ross} < 0.0) in the temporal and horizontal averaged entropy (see Fig. 24). We determined the entropy jump from the difference of the entropy minimum and the fixed entropy at the bottom, i.e. (16)In order to calculate s_{min}, we used averages of the entropy on constant Rosseland optical depth ^{12}, since it is the radiation losses that cause the sharp changes in the thermodynamic state around the optical surface and, therefore, the optical depth scale offers a better reference frame for comparisons. The averages on constant geometrical height ⟨3D⟩_{z} smear out and thereby overestimate s_{min} increasingly towards higher T_{eff} due to the increasing level of corrugation of isos surfaces on the geometrical scale (see Fig. 24). The constant entropy at the bottom s_{bot}, on the other hand, is a fixed input value for each simulation. It is worthwhile to mention that the main contribution to the variation in Δs as a function of stellar parameters is due to s_{bot}, since s_{min} varies just slightly with stellar parameter compared to s_{bot} (see Fig. 5).
In Fig. 5 we show s_{min} (dashed lines) as well as Δs (dotted lines), and it is obvious that the minimum in entropy increases just slightly with increasing T_{eff}, while the jump increases with the constant entropy value of the adiabatic convection zone quasiexponentially at higher T_{eff} and lower log g. This can be concluded more easily from Fig. 18 (top panel), where we display Δs vs. T_{eff} (see also Col. 8 of Table C.1 and for Δs and s_{min} in Appendix B). We note that the location of the entropy jump essentially represents the boundary of stars and the jump is to be regarded as physically realistic, which is a result of 3D RHD simulations. The entropy minimum coincides with the position of the upper end of the SAR. A similar sharp drop occurs for most of the thermodynamic quantities of interest (ε, T and n_{el}), whereas ρ and p_{tot} display a marked change of gradient. Moreover, the jump in entropy is an important value, since it is a direct measure of the efficiency of convective energy transport (see Trampedach et al. 2013). The latter is in 1D modeling set by the four MLT parameters, especially α_{MLT}, in the framework of MLT (see BöhmVitense 1958; Henyey et al. 1965). Towards cool dwarfs Δs becomes smaller, indicating a higher convective efficiency, while towards hotter stars the large entropy jumps reflect a low convective efficiency. We will present a calibration of α_{MLT} based on the entropy jump in a subsequent study, as previously done by using multidimensional convection simulations (see Ludwig et al. 1999; Trampedach et al. 1999). While we found s_{bot} to be qualitatively similar to those of Ludwig et al. (1999), in the case of Δs our findings are also similar, but the differences are here distinctively larger.
Fig. 7 Comparison of the entropy jump Δs against the constant entropy value of the adiabatic convection zone s_{bot} for all grid models (filled circles; colorcoding explained in the legend). Furthermore, we show also hyperbolic tangent functional fits (see Eq. (B.4) and Appendix B) between T_{eff} = 4000 and 6000 K (red lines; top panel) and between log g = 1.5 and 4.5 (grey lines; bottom panel). 
It is quite remarkable how closely the variation of Δs resembles the change of s_{bot} with stellar parameters (see Fig. 5). Motivated by this, we compare directly Δs against s_{bot} in Fig. 7. We find a nice correlation between Δs vs. s_{bot}. At lower s_{bot} values, Δs seems to converge towards 0.0 (a negative jump is not expected, since the atmosphere is losing energy in form of radiation from the photosphere), while for s_{bot} ≳ 1.7, Δs grows linearly with s_{bot} and only a modest level of scatter. In Fig. 7, we colorcoded the T_{eff} and log gvalues respectively, to show how the residuals depend systematically on atmospheric parameters. Models with higher T_{eff} (bright orange dots) and higher log g (dark grey dots) settle along higher Δs and vice versa. In order to illustrate this better, we have fitted a set of hyperbolic tangent functions (see Eq. (B.4)), which we show also in Fig. 7. We included functional fits between T_{eff} = 4000 and 6000 K (red/orange lines in top panel) and between log g = 1.5 and 4.5 (grey lines in bottom panel). Hence, we find hotter dwarfs along lines at larger Δs, while cooler giants settle along lines at smaller entropy jumps.
Interestingly, in the linear part (s_{bot} ≳ 1.7) of displays a rather universal slope of Δs/s_{bot} ~ 0.85, while the actual offset in the ordinate depends mainly on T_{eff} and log g. Another interesting aspect is that T_{eff} shows a similar strong influence as log g. The latter, however, is obviously expressed in logarithmic scale, therefore the influence of T_{eff} is much stronger. On the other hand, when one performs a similar hyperbolic tangent functional fit for a fixed value of [Fe/H], then Δs is dispersed around the functional fit with such a large scatter that a fit is rather meaningless. Therefore, in contrast to T_{eff} and log g we find no systematic trends with metallicity.
Based on the strong correlation between the entropy jump Δs and s_{bot}, it is of interest to investigate what other scaling relations may be manifested for other stellar properties. With Δs as an inverse measure of convective efficiency, we expect that in light of such scaling relations, important quantities depending on the entropy jump will also similarly scale systematically with s_{bot}, in particular density and velocity (see Sects. 3.2.2 and 3.2.4), and therefore also the calibrated mixinglength of a particular MLT implementation. We note briefly that qualitatively similar relations can be achieved with 1D MLT models with a fixed mixing length.
3.1.4. Emergent intensity
Fig. 8 Overview of the emergent (bolometric) intensity for a selection of stars, namely mainsequence (MS), turnoff (TO), Kgiant and Kdwarf (from left to right, respectively) at a given time instant. For each star, we show four metallicities [Fe/H] = 0.0, − 1.0, − 2.0 and − 3.0 (from top to bottom, respectively). To facilitate comparisons between the different metallicity of each star, the intensity scale and the horizontal geometrical size of the metalpoor simulations are identical to [Fe/H] = 0.0, and the individual intensity contrasts [%] are indicated in each box. 
While classic 1D models are inherently horizontally symmetric, therefore lacking a visible granulation pattern, the emergent intensity of 3D models features inhomogeneities exhibiting rich details, which arise due to the presence of turbulent convective motions. We give an overview the emergent intensity of our simulations in Fig. 8. Therein we display a mainsequence (MS) simulation (the Sun), a turnoff (TO) simulation, a Kgiant, and a Kdwarf model, each with four different metallicities. To facilitate direct comparisons among the four metallicities, we kept the horizontal sizes and the color scales for the continuum intensities fixed from [Fe/H] = 0.0 for the individual stellar categories (we extended the metalpoor simulations by exploiting the periodic horizontal boundary conditions). The dark regions depict the cold intergranular lanes, while the brighter areas are the hot granules. The radiation above the granules originate at higher geometrical heights, while for downdrafts it comes from much lower heights. This is because the opacity is highly nonlinear due to the strong temperature sensitivity of the H^{−}opacity (κ_{H−} ~ T^{10}, see SN98), which is the by far the dominant continuum opacity source in the visible for latetype stellar photospheres. Since the temperature difference between the granules and the intergranular lanes is very large (>10^{3} K), layers of constant optical depth will be increasingly more corrugated and become largest around the peak of the SAR. Therefore, the radiation above granules is emerging from higher geometrical depths, z_{up}, while above downdrafts it originates from deeper geometrical heights, z_{dn} (for the Sun the largest difference between the averaged geometrical heights can amount up to ⟨z⟩_{up} − ⟨z⟩_{dn} ≃ − 140 km at log τ_{Ross} = 2.0).
An immediate, interesting aspect that leaps to the eye from the overview presented in Fig. 8 is the qualitative selfsimilarity of the granulation patterns despite the large variations in sizescales. The emergent intensity increases towards higher T_{eff} and decreases for lower surface gravities, as expected. From Fig. 8, it is also clear that the granule sizes decrease with metallicity (due to smaller H_{P}, see Sect. 3.1.5; see also Collet et al. 2007). Also apparent is the change of intensity contrast with stellar parameters, as we will discuss below.
In order to discuss the changes in the intensity, we show in Fig. 9 the temporally averaged histograms of the intensity I normalized to their individual mean intensity , thereby enabling a direct comparison between different stellar parameters. The histograms of the intensity show two components: a peak at lower (darker, ) intensities, resulting from the cool downdrafts, and an often broader component at higher (brighter, ) intensities, arising from the upflowing hot granules. We note that these findings are qualitatively to be expected (see SN98; Trampedach et al. 2013).
As clearly depicted in Fig. 9, the shapes of the two components change with stellar parameters, in particular the amplitudes and widths, thereby changing the overall shape. The two components can be clearly extracted from histograms at higher T_{eff}, where the intensity contrast is increasingly enhanced and eventually produces a distinctly bimodal distribution, which is a manifestation of the hidden or naked granulation (see Fig. 10 and Nordlund & Dravins 1990). In order to better illustrate this, we also included the full width at half maxima (FWHM) of the individual intensity histograms I_{FWHM} in Fig. 10. On the other hand, at lower T_{eff} the intensity contrast decreases in general, so that the two components overlap, leading to a single narrower higher peak in the histogram, thereby becoming indistinguishable from each other in the histogram. Ludwig & Kučinskas (2012) found also an unimodal intensity distribution in the context of a 3D giant model with solar metallicity. Furthermore, we find that the individual contribution to the intensity from upflows and downflows is often asymmetric, meaning that the amplitudes of the two peaks in the bimodal distribution are unequal (see Fig. 9). In general, for dwarfs, we find that the relative importance of downflows with respect to upflows in terms of the peak contribution to the intensity distribution increases with increasing T_{eff}. However, we also find exceptions, e.g. at lower metallicity where the behavior at T_{eff} = 4500 K is actually the opposite. Also, the balance between upflows and downflows varies with surface gravity. The intensity histograms for giants are in general broader (higher contrast) compared to dwarfs of the same T_{eff} (see I_{FWHM}), hence exhibiting a larger intensity contrast. For dwarfs at lower metallicity (right bottom panel) the bimodality is more pronounced and the I_{FWHM} (contrast) is broader (higher) towards higher T_{eff}, while at lower T_{eff} the I_{FWHM} (contrast) becomes narrower (lower) compared to solar metallicity (left bottom panel). The latter hints at an enhancement of the effect of hidden or naked granulation.
Fig. 9 Temporally averaged histograms of the (bolometric) intensity (solid lines) for various stellar parameters (the binsize is 100). To ease comparisons between the different stellar parameters, we normalized the individual intensity scales with its mean value. Furthermore, we have indicated the respective FWHM of the intensity histograms, I_{FWHM}, (dotted lines), a measure of the intensity contrast. Note the different ordinate scale in the top panel. The bimodal distributions are given due to the asymmetry and the large contrast in the up and downflows. 
Fig. 10 Overview of the (bolometric) intensity contrast ΔI_{rms} against T_{eff} for [Fe/H] = −2.0 and 0.0 (blue and black respectively). Models with the same gravity, but different T_{eff} are connected. The enhanced naked or hidden granulation at lower metallicity can be extracted in the larger range of intensity contrasts (see text for details). 
To illustrate the latter in more detail, we show in Fig. 10 the rms of the bolometric disk center intensity fluctuations for [Fe/H] = 0.0 and − 2.0, which is commonly referred as the intensity contrast (17)with being the (spatial) mean intensity and N the number of data points (see Roudier & Muller 1986). We remark that the shown ΔI_{rms} are is temporal averages. It is essentially defined as the relative standard deviation, hence it reflects the width of the intensity distribution (see Fig. 9). This often measured value is very suitable for quantifying the range of brightness fluctuations due to granulation. The intensity contrast increases with higher T_{eff} and lower log g. For our solar simulation we get an intensity contrast of 15%, which is close to the one found by SN98 with 16% (see Col. 10 in Table C.1).
Towards higher T_{eff}, we find that s_{bot}, Δs, and the vertical velocity increase, as shown in Sects. 3.1.2 and 3.2.2. For increasingly hotter stars, the top of the convective zone, z_{top,cz}, penetrates higher and higher above the optical surface due to larger vertical velocities (see Fig. 17). Additionally, at higher T_{eff} (higher s_{bot} and Δs), the overall temperatures and their fluctuations also increase, implying that one observes increasingly higher layers, since the dominant H^{−}opacity, hence the optical depth, depends sensitively on the temperature. Therefore, the granulation pattern is enhanced at higher T_{eff}, while on the contrary for lower T_{eff} the granulation becomes less visible, since z_{top,cz} recedes below the optical surface in the latter case (see overview in Fig. 8). This phenomenon has been already described by Nordlund & Dravins (1990) as naked granulation.
Interestingly, in our simulations, we find that at lower metallicity the effect of naked and hidden granulation is more pronounced, in the sense that the range of contrast from cool, lowcontrast dwarfs, to hotter, highcontrast dwarfs, is 61% larger (from 10.9 to 17.6) for our [Fe/H] = −2.0 simulations than for solar metallicity (see Fig. 10). At lower metallicity the major electron donors (metals) are depleted, therefore the formation of the dominant opacity source H^{−} depends primarily on the ionization of hydrogen, which is the reason for the steep increase of intensity contrast with higher T_{eff} (Nordlund & Dravins 1990).
The variations in the intensity and in the intensity contrast with stellar parameters have important ramifications for observations. At the one hand, the enhanced naked or hidden granulation at lower metallicity affects the formation of spectral lines and the limb darkening. On the other hand, it should also lead to distinct signatures in the granulation background of asteroseismological observations and spectrointerferometric imaging. Here, we limit ourselves to the discussion of the global properties of the emergent intensity patterns, and a detailed analysis will be performed in subsequent papers.
3.1.5. Granule size
Fig. 11 Overview of the granule diameter d_{gran} derived from the maximum of the mean 2D spatial power spectrum of the bolometric intensity against T_{eff} for [Fe/H] = −2.0 and 0.0 (blue and black respectively). Models with the same gravity are connected by their respective functional fits (solid lines; see Appendix B). 
The physical dimensions of the simulations boxes (s_{x},s_{y} and s_{z} in Cols. 11 and 12 of Table C.1) are selected based on the mean diameter of granules (see Sect. 2.3.1 and Table C.1) and a target of about 10 granules in the box. Additionally, we measured the granule sizes by calculating the 2D spatial power spectrum of the bolometric intensity for the time series, and determining its maximum from the smoothed time average (see Fig. 9 in Trampedach et al. 2013). This method is quite robust despite the large variations in gravity. In Fig. 11, we present the measured granule sizes d_{gran} (given in Col. 13 in Table C.1; see also Appendix B), showing that they become larger with smaller surface gravity. Also, the granules of the simulations with fixed log g, the lowest T_{eff} are typically ~ 50% smaller compared to the simulations with the hottest T_{eff}, while for the models with the lowest metallicity they are typically ~ 30% smaller than for the metalrich ones.
Fig. 12 Top panel: comparison of the granule size approximated with the pressure scale height d_{NSR} (see Eq. (18)) against the mean granule diameter d_{gran} (same as Fig. 11) for all models. The individual stellar parameters are indicated (T_{eff}, log g and [Fe/H] with red, gray, and blue respectively). We indicated the position of d_{NSR} = d_{gran} the solid diagonal line. Bottom panel: relative residuals. We indicated also the mean residual (dotted line), which amounts to ~ 30%. Here, models with the same gravity are connected (solid grey lines). 
We find a remarkable validation for the approximation of the maximal horizontal extent of a granule based on massconservation considerations made by SN98 (see also Nordlund et al. 2009). Hereafter, we denote the following relation as the Nordlund scaling relation (NSR). The ascending buoyant plasma inside a cylindrical granule with radius r gives rise to a vertical mass flux with j_{z} = [πr^{2}] ρv_{z}. This mass flux has to deflect and overturn increasingly towards the top. Due to conservation of mass, the upflow has to drain off sideways through the edge of the granule within approximately one pressure scale height H_{P}, hence resulting in a horizontal mass flux j_{h} = [2πrH_{P}] ρv_{h}. The pressure is a quantity that preserves its characteristic shape with stellar parameters, i.e. the pressure of two different simulations look rather similar on a uniform depth scale, therefore the pressure scale height is preferred over the density scale height. Equating j_{z} and j_{h} we can solve for the (maximal) granular diameter, d = 2r: (18)We show in Fig. 12 a comparison of the granular diameters estimated with d_{NSR} and from the maximum of 2D spatial power spectra d_{gran}, which is shown in Fig. 11. The astonishing tight correlation can solely be interpreted as clear evidence for the validity of the NSR. We find that the mean pressure scale height taken at the height of the maximum vertical rmsvelocity below the optical surface (log τ_{Ross} ~ 2.0, see Fig. 17) gives the best match between d_{NSR} and the granule sizes. Furthermore, we also confirm that the relevant scaleheight is that of the total pressure scale height, H_{P} = p_{tot}/ρg, since we find a better agreement with the latter. The granular diameters found from the peak of 2D spatial power spectra are about 30% larger than the estimate from Eq. (18), i.e., d_{gran} ~ 1.3d_{NSR} (see lower panel of Fig. 12). The variation of the velocity ratio v_{h}/v_{z} in the convection zone is rather small (v_{h}/v_{z} ~ 1.0) as both are of the order of the sound speed, therefore the variation in Eq. (18) stems predominantly from H_{P}. In hydrostatic equilibrium the pressure scale height is inversely proportional to the surface gravity (H_{P} ∝ 1/g), which explains the strong correlation between the granular sizes and log g. On the other hand, with increasing T_{eff} and [Fe/H], the pressure scale height increases slightly because of the increase in the ratio of pressure and density (H_{P} ∝ p_{tot}/ρ). The ratio actually increases even though both values decrease, since the density drops with height slightly more rapidly than the pressure.
Finally, we want to mention our finding on the filling factor for upflows and downflows, f_{up} and f_{dn} respectively. We derived the filling factor from the sign of the velocity field in the unaltered simulations on layers of constant geometrical height. Then we computed the mean filling factor in the convection zone, which yields on average for all simulations f_{up} ≃ 0.65 with a minute deviation of σ = 0.014. Therefore, we find that the mean filling factor is rather universal, and close to previous findings by SN98 with f_{up} ~ 2/3 and f_{dn} ~ 1/3. In deeper solar simulations, which reach down to 20 Mm (Stein et al. 2011), we find very similar values for the filling factor.
3.2. The mean atmosphere
In the following, we want to discuss the properties of the mean stratifications and the temporal and spatial averages of various important quantities. Unless specified otherwise, the ⟨3D⟩ stratifications presented here are averages on surfaces of constant Rosseland optical depth, i.e. ⟨3D⟩ = ⟨3D⟩_{Ross}. Whenever we employ alternative averages in the text, e.g., on constant geometrical height ⟨3D⟩_{z}, we indicate that explicitly. We remark briefly that only the averages on constant geometrical height ⟨3D⟩_{z} strictly fulfill the equations of conservation (Eqs. (1)–(3)), therefore also the hydrostatic equilibrium, while all other averages exhibit slight deviations. This and the actual methods we used to compute the mean stratifications are discussed in a separate paper (see Magic et al., in prep.). For the sake of clarity, we display here only a subsample of our grid models including MS and RGB stars (log g = 4.5 and 2.0, which we refer to as dwarfs and giants, respectively) with solar and subsolar metallicity ([Fe/H] = 0.0 and − 2.0, solar and metalpoor, respectively). Whenever possible, we compare with corresponding 1D models that are obtained with our 1D code (see Appendix A).
Before continuing our discussion, we would like to point out the importance of the superadiabatic region (SAR), as it will be referred repeatedly in the following. It is the region, where the transport of energy changes character, from convective to radiative. The top of the SAR, where ∇_{sad} = 0, marks the top the convection zone, since it is the uppermost point, where the Schwarzschild criterion is fulfilled. At the location of the peak of the superadiabatic gradient, one also finds the largest fluctuations and inhomogeneities in the thermodynamic variables due to the nonadiabatic transition to the photosphere. Furthermore, it is here in the SAR, where the entropy jump and the peak in the vertical velocity occur. In fact, the SAR effectively represents the physical outer boundary of the convective envelope. It is the most dynamic part in the interior of latetype stars, where the largest fluctuations are found. This is the reason why hydrostatic 1D modeling has the greatest challenges in this rather small region.
3.2.1. Temperature stratification
Fig. 13 ⟨3D⟩ stratifications of the temperature vs. optical depth for various stellar parameters (solid lines) and 1D models with α_{MLT} = 1.5 (dashed lines). Furthermore, we have marked the positions of entropy minimum (plus), vertical peak velocity (diamond), and maximum in ∇_{sad} = ∇ − ∇_{ad} (triangle). 
We first consider the temperature stratifications in optical depth, which we show in Fig. 13. We also show the corresponding stratifications of 1D theoretical model atmosphere with α_{MLT} = 1.5 based on our 1D code (dashed lines) with identical EOS and opacity tables as the 3D models. In the continuum forming layers around the optical surface (− 1.0 < log τ_{Ross} < 0.5), the differences between ⟨3D⟩ models with different T_{eff}s, but same log g and [Fe/H], are rather small besides the shift in the temperature stratification corresponding to the difference in effective temperature ΔT_{eff}, which is to be expected since T_{eff} ≈ T(τ = 2/3). Well above and below the optical surface, on the other hand, we find significant differences between the ⟨3D⟩ models depending on the stellar parameters.
In the upper layers (log τ_{Ross} < − 2.0) of atmospheres with solar metallicity, we find that the behavior in mean temperature is similar between 3D and 1D models. On the other hand, the metalpoor ⟨3D⟩ models exhibit significantly cooler temperature stratifications compared to ⟨3D⟩ models with solar metallicity (ΔT/T(log τ_{Ross} = −0.5) ~ − 1% and ~ − 14% for [Fe/H] = −1.0 and − 3.0 respectively), in particular for dwarfs (log g = 4.5). The temperature stratification in the upper photospheres of solarmetallicity models is largely controlled by radiative equilibrium, while for lowmetallicity models this is not generally the case: for metalpoor models, the absorption features become considerably weaker, therefore, the radiative heating by spectral line reabsorption (q_{rad}) is dominated by the adiabatic cooling due to expansion of the ascending gas (− p_{th}∇·v) in the energy balance (Eq. (3)), leading to an equilibrium structure at cooler temperatures (Asplund et al. 1999b). For cool, metalpoor giants (e.g., T_{eff} = 4000 K, log g = 2.0), we recognize the effects of molecule formation on the structure of the high atmosphere. At sufficiently low temperatures, molecules start to form, which contribute with a large line opacity, shifting the balance from adiabatic to radiative heating and cooling, resulting in a stratification closer to the radiative equilibrium one (see Gustafsson et al. 2008). On the other hand, for giants with solar metallicity the radiative equilibrium is even more dominating, since these exhibit hotter stratifications than 1D models in the upper layers. Ludwig & Kučinskas (2012) find the same but on a much milder level. These effects are rather nonlinear, and we find no simple systematic trends within our grid models.
Fig. 14 ⟨3D⟩ temperature stratifications with the variation of one stellar parameter at a time, while the two others are fixed (T_{eff}, log g, and [Fe/H], from top to bottom, respectively). We show our standard averages on constant optical depth ⟨3D⟩ (solid line) and on constant geometrical depth ⟨3D⟩_{z} (dotted line). 
We would now like to examine the influence of individual stellar parameters on the temperature stratification. Therefore, we show in Fig. 14 the temperature stratifications of models where we separately vary one at the time (T_{eff}, log g, and [Fe/H]), while keeping the other two parameters constant. Figure 14 (top panel) shows, as expected, that with increasing T_{eff}, the temperature stratification becomes overall hotter above, but also below the optical surface, in order to provide the required total energy flux (higher enthalpy, Eq. (22)). We find in our simulations (both 1D and 3D) that the increased T_{eff}s with hotter stratifications are accompanied by lower densities and higher vertical velocities below the surface (see p_{tot}(log τ_{Ross} = 2.0) in Fig. 22 representatively for the ρ). The net effect on the convective flows are lower mass fluxes for higher T_{eff}s, since the decrease in density is predominating the increase in velocity, therefore resulting in a more inefficient convection. This is compensated with higher entropy jumps (see Fig. 18 with Δs as an inverse measure for convective efficiency), hence higher temperatures and steeper temperature gradients. On the other hand, the temperatures in the upper, radiative layers increase less with increasing T_{eff} than in the deeper, convective ones. We find with decreasing surface gravity (middle panel in Fig. 14) the same correlations as with increasing T_{eff}s before, the temperature stratifications become hotter below the photosphere, and due to lower densities we find a more inefficient convection, while the upper atmosphere is less affected. For lower metallicities (bottom panel), the temperature stratifications are significantly cooler, both above and below the optical surface (ΔT/T(log τ_{Ross} = 2.0) ~ − 5% and − 15% for [Fe/H] = −1.0 and − 3.0 respectively). At the top the stratifications are cooler at lower [Fe/H] due to the dominance of adiabatic cooling over radiative heating. Below the optical surface, we find higher densities with lower velocities and entropy jumps (while the mass flux is increasing), therefore, leading to an efficient convection with shallow temperature gradients at lower metallicities. We find cooler models that fall below the opacity edge, which we describe below (compare T_{eff} = 5000 K in Fig. 16), follow an adiabatic temperature stratification even in the atmosphere, which coincides with the rather sudden change between [Fe/H] = −1.0 and − 2.0 in Fig. 14 (bottom panel). Besides our standard averages on constant Rosseland optical depth, we show also the averages on constant geometrical depth scale ⟨3D⟩_{z} (here is z fixed and τ_{Ross} = ⟨τ_{Ross}⟩_{z}), which are systematically different, in particular below the optical surface, but behave qualitatively in a similar way with stellar parameters.
Fig. 15 Overview of the maximum temperature gradient ∇_{peak} (top panel) and Rosseland opacity κ_{Ross} taken at the height τ_{Ross} ≈ 3.0 (bottom panel) against T_{eff} for [Fe/H] = −2.0 and 0.0 (blue and black respectively). Models with the same surface gravity are connected by their respective functional fits in the top panel (solid lines; see Appendix B). 
Fig. 16 Mean internal energy against mean density for dwarf models (log g = 4.5) with [Fe/H] = 0.0 and − 2.0 (left and right respectively). The specific isocontours for the entropy s_{bot} (green solid) and s_{min} (green dotted), Rosseland opacity per volume, ρκ_{Ross}, (blue) and temperature T (orange) are underlayed. Moreover, the positions of entropy minimum s_{min} (plus), optical surface (large square), vertical peak velocity (diamonds) and maximum in ∇_{sad} (triangle) are marked respectively. The amplitude of is indicated with horizontal bars with markings in 1 km s^{1}. We also included the 1D models with α_{MLT} = 1.5 (blue dashed lines). The range in optical depth is shown from log τ_{Ross} = −5.0 to + 5.0 for each dex (small squares). However, we note that our simulations boxes are much deeper (⟨log τ_{Ross}⟩ ≈ + 7.5). 
In the subphotospheric region (log τ_{Ross}> 0.5), where convection dominates, the temperature gradients ∇ = dlnT/dlnp_{tot}^{13} become increasingly steeper with higher T_{eff}, reflecting the hotter interior stratifications. This can be illustrated with the maximum in temperature gradient, ∇_{peak} = max∇, which we show in Fig. 15 (functional fits are also given in Appendix B). The increase of ∇_{peak} with T_{eff} is close to linear, but it seems to saturate at higher T_{eff} (see T_{eff} ≥ 6500 K). We find that the maximum in temperature gradient ∇_{peak} reproduces qualitatively a similar behavior as the intensity contrast ΔI_{rms} with stellar parameter (compare Figs. 10 and 15), which is consistent with the strong temperature sensitivity of the H^{−} opacity. Furthermore, our metalpoor simulations exhibit a larger range of ∇_{peak}values than their solar metallicity counterparts, and ∇_{peak} is similarly enhanced at lower metallicity (compare also T(log τ_{Ross} = 2.0) for dwarfs in Fig. 13), as the intensity contrast (see Sect. 3.1.4). Our cool metalpoor simulations have flatter and hot metalpoor simulations have steeper temperature stratifications than the metalrich part of our grid (see Fig. 14). Curiously, ∇_{peak} is close to constant with metallicity for the solar T_{eff} and log g.
We identify three main reasons for the given variations in temperature gradients with stellar parameters in the SAR, which are rooted in the hydrodynamics and the radiative transfer: velocity field, convective efficiency, and radiative backwarming.

1.
As we discussed above (see Sect. 3.1.3), theentropy jump Δs increases with effective temperature according to a power law (see Fig. 18). This behavior arises due to the variations in the radiative losses (see Sect. 3.2.8), which is accompanied by changes in internal energy and density (see Fig. 16). The velocities rise rapidly, as exhibited by the growth of and also with T_{eff} (see Figs. 17 and 21 respectively). Similar to ∇_{peak}, both and occur in the SAR, and both increase towards higher T_{eff} and lower log g.

2.
The mass mixing length α_{m} changes with stellar parameters (Trampedach & Stein 2011). The latter is evaluated as the inverse gradient of the vertical mass flux, separately in the up or downflows, hence , with j_{z,up} being the vertical mass flux in the upflows. Therefore, the mass mixing length is composed of the gradients of the density and the vertical velocity, i.e. α_{m} ∝ 1/dlnρ + 1/dlnv_{z}. We find that α_{m} increases for lower Δs, and ∇_{peak}. We will publish our findings on the mass mixing length in a separate paper.

3.
In the lower photospheric layers, where the continuum forms, radiation is absorbed (blocked) by spectral lines; this implies that less radiative flux can be transported at the wavelengths corresponding to spectral lines and, conversely, that more flux has to be pushed through continuum windows, an effect commonly referred to as lineblanketing. This in turn leads to a steepening of the temperature gradient and to additional heating of the subsurface layers, also known as backwarming (see Mihalas 1970; Nordlund et al. 2009). In 1D models it is straightforward to quantify the effect of backwarming, as done for example by Gustafsson et al. (2008), who found it to contribute a slight increase in temperatures below the surface ( for solar metallicity stars with T_{eff} ≈ 5000 K and log g = 3.0). In our 3D RHD atmosphere models, lineblanketing and backwarming effects are also naturally included through our opacitybinning method. Isolating the radiative backwarming effect in our 3D simulations is, however, a little more involved than in 1D and we defer the analysis of this mechanism to a future paper in this series.
The three mentioned effects are nonlinearly coupled and compete with each other, making it difficult to disentangle the individual contributions. A quantitative analysis will be presented in a later paper.
We would like now to examine more closely a sample of ⟨3D⟩ models in the ε − ρplane, as shown in Fig. 16, in order to better illustrate the variations with stellar parameters. One can clearly distinguish three different regimes: the adiabatic convection zone, the photospheric transition, and the almost isothermal upper atmosphere.
At the bottom boundary, sufficiently deep in the convection zone where entropy fluctuations become small, ⟨δs⟩ = 0.3%, the models follow closely the associated adiabats with s = s_{bot} (green lines). They deviate increasingly from their adiabats, as the top of the convection zone is approached. This is due to the entropy deficient downdrafts (cooled in the photosphere) becoming less diluted by the entropic upflows, as the optical surface is approached. For the 1D models (blue dashed lines), the value of the entropy at the bottom of the stratifications is evidently overestimated, particularly at higher T_{eff}, but this is because we haven’t calibrated the α_{MLT} parameter here, and we have used a value of 1.5 for all models. The transition of energy transport from fully convective to fully radiative is clearly visible, since, at the optical surface, one finds a sharp isochoric (Δρ ~ 0.0) drop in internal energy (this is basically the enthalpyjump Δh in Eq. (22)). The εjump coincides with the subphotospheric region (0.0 < log τ_{Ross} < 2.0), where the atmosphere starts to become transparent. The transition zone ends eventually at the optical surface (log τ_{Ross} ≃ 0.0, marked with big squares). Above the optical surface (log τ_{Ross} < 0.0) the atmosphere is almost isothermal (compare with the orange isotherms in Fig. 16), with exponentially decreasing density and almost constant internal energy (Δε ~ 0.0).
The entropy s_{bot} at the bottom grows exponentially with increasing T_{eff} and decreasing log g. We showed above that the entropy jump increases in a similar way (see Fig. 7). Here we find a similar behavior for the jump in internal energy (Δε, hence Δh) in the photosphere. Moreover, we show in Fig. 16 the positions of and located in the εjump, and, again, we find both and to scale exponentially with T_{eff}. For we have indicated the amplitudes as well, which also increases exponentially with T_{eff}. All of the aforementioned positions are distributed rather regularly in the ε − ρplane, while they are less so on the log τ_{Ross}scale (see Fig. 13). The position of s_{min} is close to the optical surface and shows little variation in optical depth.
At lower energies and densities in Fig. 16 (log ε ~ 0.3 and log ρ ~ 3.0 to 0.0) we notice the effect of H i and He i opacity in form of an edge in the opacity contours (log κ_{Ross} ~ − 5.0 to − 1.0), since the boundfree absorption increases (more excited states) towards higher energy below the ionization energy, and they fade away again above it. Models that fall below this edge exhibit a rather different stratification. In particular, towards cool metalpoor dwarfs, i.e. lower T_{eff}, higher log g, and lower [Fe/H], the models more closely follow adiabats than isotherms, in the atmosphere. This effect of the competition between radiative and dynamic heating (see beginning of this Sect.) above the convection zone becomes particularly evident at lower metallicity (for [Fe/H] ≤ − 2.0). However, for the 1D models (blue dashed lines), this is obviously not the case, since these always follow isotherms due to the enforcement of radiative equilibrium. Furthermore, the cool metalpoor ⟨3D⟩ models also display higher densities at the optical surface, thereby spanning a larger ρrange for different T_{eff}s. The stratifications of simulations of hotter dwarfs, on the other hand, depend little on metallicity. For the simulations, we have only plotted the range log τ_{Ross} ∈ [− 5.0,5.0], and the top of this is reached at much higher densities for the metalpoor dwarfs than for the solar metallicity dwarfs. Therefore, the density ranges covered above the optical surface by the individual atmospheres is small for metalpoor models (min[log ρ] ~ − 1.0 and − 2.0, respectively; see Fig. 16).
3.2.2. Velocity field
Fig. 17 Vertical rmsvelocity, v_{z,rms}, from the ⟨3D⟩ stratifications (solid lines) and convective velocity v_{MLT} from the corresponding 1D MLT models (dashed lines) as function of optical depth log τ_{Ross} for various stellar parameters. 
Next, we consider the velocity field in our simulations, which arise selfconsistently from the solution of the hydrodynamic equations. In Fig. 17 we show the rmsvelocity of the vertical component v_{z,rms}, being the flux carrying component of the convective flows, and being the broadening component of spectral lines at disk center.
The buoyant uprising plasma will experience increasingly a decrease in density towards the photosphere, hence a strong density gradient, and due to mass conservation, the convective motions will eventually overturn. Therefore, v_{z,rms} peaks in the SAR around log τ_{Ross} ~ 1.5 for dwarfs and ~ 2.3 for giants. Furthermore, since in the SAR, the transition region from convective to radiative transport of energy takes place due to decrease in opacity and the subsequent radiative losses, here we find the strongest turbulent motions concomitant with the greatest fluctuations in all thermodynamical quantities (in particular entropy, temperature and density, see Figs. 24, 13 and 18). Further towards the interior, v_{z,rms} drops as the density increases. From the slightly subphotospheric maximum, velocities fall off to a minimum above the optical surface, then increases again in the higher atmosphere (see Fig. 17). Towards upper layers, v_{z,rms} increases again due to pmodes, excited in the SAR but leaking out of the acoustic cavity as they have frequencies above the acoustic cutoff. The metal poor simulations show a slightly smaller increase in v_{z,rms}, since their density gradients are shallower due to steeper Tgradients. The declining velocity above the surface is due to the fact that the convective motions overshoot well above the top of the convection zone. We find the velocity minimum to occur between log τ_{Ross} ~ − 2.3 for dwarfs and − 1.5 for giants.
As expected, the magnitude of the velocity field is enhanced towards higher T_{eff}, lower log g, and higher [Fe/H], similar as Δs. The symmetry of the velocity profile changes with log g and metallicity, while it is little affected by T_{eff}. For lower log g, the peak in the velocity field is increasingly shifted to optically deeper layers (e.g. at solar metallicity the average peak position for dwarfs is ⟨log τ_{Ross}⟩ ~ 1.5, while for giants it is ~ 2.0). The coolest metalpoor simulations display a flatter profile, and the position of the minimum is increasingly shifted towards higher layers, especially for extreme metalpoor dwarfs ([Fe/H] < − 3.0), and the profile is therefore stretched and skewed.
For comparison, we also show in Fig. 17 the convective velocity v_{MLT} of our 1D models determined by MLT. It is apparent that the general trends of increasing velocities with increasing T_{eff} and [Fe/H] and decreasing log g, are common between the simulations and the 1D MLT models, although much less pronounced in 1D. Furthermore, v_{MLT} drops rather sharply at the top of the convection zone (as given by the Schwarzschild criterion), as no overshooting is allowed for in our implementation of MLT. Several nonlocal variants of MLT exists, and they allow for overshooting, but none of them produce velocity profiles close to that of our simulations. We also want to mention the large asymmetry in velocities of the up and downflows (SN98): in 3D simulations, the latter are much faster than the former (up to ~ 2 faster, in particular for cool dwarfs), contrary to what is normally assumed in 1D descriptions of convection.
Fig. 18 Overview of the entropy jump (top panel), maximal vertical rmsvelocity below the surface (middle panel) and the density at the same height ρ_{peak} (bottom panel) vs. T_{eff} for [Fe/H] = −2.0 and 0.0 (blue and black respectively). Models with the same gravity are connected with their respective functional fits (solid lines; see Appendix B). Note the inverse correlation between density and velocity. 
The peak vertical rmsvelocity, (see bottom panel of Fig. 18), is a good measure of the global magnitude of velocities in the simulations. It also serves as a measure of the amount of turbulence present in the simulations. The actual values are also given in Col. 9 in Table C.1, together with a functional fit in Appendix B. The variation of with stellar parameters resembles that of the entropy jump Δs (compare top and bottom panel in Fig. 18), namely it also increases exponentially with higher T_{eff} and lower log g and linearly with [Fe/H]. An interesting aspect is the increase of and Δs with T_{eff}, which are close to exponential, indicating a correlation between the two. The characteristic variations of correspond to the inverse variations of the density taken at the same heights as the peaks in v_{z,rms} (see middle panel in Fig. 18). This behavior arises due to conservation of mass (Eq. (1)), which can be expressed as (19)for a stationary flow (∂_{t}ρ = 0). Of course, under this assumption, this equation is strictly speaking valid only locally, while we compare here averaged values. Despite that, we find that this relation is, in general, qualitatively fulfilled across the whole depth of our simulations. Towards the optical surface, the density decreases, which has to be compensated by faster velocities, in order to fulfill conservation of mass as well as sustain the energy flux. The velocity field profile results ultimately from the interplay between the vertical and horizontal acceleration due to buoyancy and overturning respectively. The latter in turn is set by the radiative losses that arises from the prevailing opacity conditions according to the thermodynamic state of the plasma (see Sect. 3.2.8). Furthermore, one can also reason that at a higher effective temperature, hence hotter temperature stratification, the density will be lower (ideal gas gives T ∝ 1/ρ; see also middle panel in Fig. 18), however, at the same time, more energy (enthalpy) has to be carried to the surface, which necessitates a faster flow (as is given in Eq. (23)). The entropy jump, density, and velocity are coupled intimately with each other (the vertical mass flux is j_{z} = ρv_{z}). Therefore, changes in one quantity imply corresponding variations in the values of the other quantities, and vice versa. The radiative energy losses at the photospheric transition generate the entropy fluctuations according to the prevailing opacity and the irradiationduration, hence it sets the amplitude of the entropy jump Δs. On the other hand, the entropy deficient plasma with its density excess determines the buoyancy force, f_{B} ~ Δρ, and therefore the vertical velocities v_{z} of the downdrafts. The downdrafts in turn will settle the upflows in order to deliver the required convective energy flux. The subtle details in the chain of causalities are nontrivial and beyond the scope of the present paper.
Fig. 19 Comparison of the entropy jump Δs (top panel) and the density at the optical depth of the maximal vertical velocity below the optical surface ρ_{peak} (bottom panel) vs. maximal vertical rmsvelocity for all models. A linear fit is also indicated in both panels (red lines). 
Similar to our finding in Sect. 3.1.3 of a scaling relation between the entropy jump and the constant entropy value of the adiabatic convection zone , we find here again another interesting, tight scaling relations between ρ_{peak}, and Δs, which we show in Fig. 19. The values are plotted on a double logarithmic scale, to more clearly illustrate the powerlaw character of the relations. From the above discussion, it follows that the vertical velocity is also correlated with the constant entropy value s_{bot} of the adiabatic convection zone and the density. We also show linear fits of the density ρ_{peak} and entropy jump Δs as a function of in loglog scale (red lines in Fig. 19), exhibiting the slopes of and , hence a scaling with the respective slopes.
In 3D RHD simulations, the nonthermal, macroscopic velocity fields arising from convective instabilities are computed selfconsistently from first principles and therefore have an immediate physical meaning. They represent the buoyant motions associated with convection and its turbulent features, and their statistical properties carry equally important physical information as the mean temperature or density stratifications. By contrast, in 1D atmosphere modeling, a freeparameterdependent velocity field v_{MLT} is derived for the convective flux in MLT. Also, for radiative transfer and spectral line formation calculations, two adhoc Gaussian velocity distributions – the socalled micro and macroturbulence (ξ_{turb} and χ_{turb}, respectively) – are usually introduced to model Doppler broadening of spectral lines associated with nonthermal (e.g., convective, turbulent, oscillatory, etc.) gas motions in stellar atmospheres. The values of the micro and macroturbulence parameters are determined by comparing synthetic and observed spectral line profiles and line strengths. Usually, a depthindependent value of the microturbulence ξ_{turb} and one global value of the macroturbulence χ_{turb} are applied in theoretical spectrum syntheses with 1D model atmospheres. Full3D line formation calculations using 3D models similar to those described here, have demonstrated that in latetype stars the required nonthermal Doppler line broadening is indeed primarily the result of Doppler shifts from the convective motions and to a lesser extent oscillations in the atmosphere (Asplund et al. 2000). As such this nonthermal velocity field is clearly depthdependent, while micro and macroturbulence are almost always assumed to be nonvarying with depth. Furthermore, v_{MLT} is solely assigned to satisfy the necessary amount of convective flux by the individual prescription of MLT. While, interestingly, v_{MLT} mimics to a certain extent the run of v_{z,rms} in the interior for cooler dwarfs. We remark that this interpretation is however not physically consistently motivated. Moreover, the convective velocity varies depending on its actual implementation (e.g. BöhmVitense 1958; Henyey et al. 1965) and, as such, v_{MLT} should be interpreted and used with caution. We point out that one important motivation for conducting 3D RHD atmosphere models is the fact that the before mentioned spurious inconsistent velocities become redundant. The hydrodynamical simulations account consistently for only one unique velocity field.
Fig. 20 ⟨3D⟩ stratification of the absolute value of the vorticity ω vs. optical depth is shown for various stellar parameters. 
Another good measure for the turbulence of a velocity field is the absolute value of the vorticity ω = ∇ × v, which is shown in Fig. 20. The vorticity arises below the surface in SAR due the overturning of the upflows and the turbulent downdrafts experiencing the density gradient. The peak in ω is associated with pronounced shear flows, which arise due deflection of the horizontal flows into downdrafts of the overturning plasma (see SN98). The vorticity is concentrated in tubelike structures in the intergranular lanes around the edges of granules. The run of the vorticity follows closely that of v_{z,rms} (see Fig. 17).
3.2.3. Turbulent pressure
Fig. 21 Fraction of turbulent pressure to total pressure p_{turb}/p_{tot} vs. optical depth log τ_{Ross} is displayed for various stellar parameters. 
The turbulent pressure, , is an additional (dynamical) pressure that arises from the (macroscopic) vertical bulk flows due to the convective motions. It appears when considering the horizontal averages of the momentum equation (Eq. (2)), more specifically of the advection term therein. The ratio of turbulent to total pressure, p_{turb}/p_{tot}, shown in Fig. 21, follows qualitatively very closely the run of v_{z,rms} (compare with Fig. 17), namely, it peaks in the SAR (log τ_{Ross} ~ 1.5), reaches a minimum around log τ_{Ross} ~ − 2.0, and increases in the upper layers (a functional fit for [p_{turb}/p_{tot}]_{peak} is given in Appendix B). In the SAR, the shape of the p_{turb}/p_{tot} profile with optical depth looks similar to a Gaussian function, however, towards lower T_{eff} and metallicity, it becomes increasingly skewed. Averages on constant geometrical depth ⟨3D⟩_{z} are similar, only the peak and the upper layers are slightly lower at higher T_{eff}s.
For hotter stars, in particular metalrich giants, the turbulent pressure becomes comparable to the gas pressure (p_{turb}/p_{tot} ~ 0.4) in the SAR, and the atmosphere is increasingly supported by p_{turb}. This means that neglecting the turbulent pressure, as is usually done in 1D models, would significantly overestimate the gas pressure. The consequence of this is a faulty, inconsistent stratification, since the overestimation in gas pressure comes at the cost of altering other physical quantities like the density, even when the temperature stratification looks similar compared to a ⟨3D⟩ model.
We find that v_{z,rms} is very close to v_{turb} = [⟨p_{turb}⟩/⟨ρ⟩]^{1/2}, since the latter is basically the densityweighted analog of the former. In 1D stellar structure models that include turbulent pressure, the convective velocity from MLT is considered, i.e. . However, as shown in Fig. 17, the convective velocities, v_{MLT}, are underestimating v_{turb} systematically towards higher T_{eff} and lower log g. Therefore, the 1D models can be improved by using v_{z,rms} resulting from 3D simulations, or one can fit the scaling factor β to match v_{z,rms}, which would clearly reduce the error.
As shown by Wende et al. (2009), we also expect v_{z,rms} to correlate with the microturbulence ^{14}ξ_{turb}, since v_{z,rms} is a horizontal average of the velocity field. In 1D line formation calculations, a depthindependent ξ_{turb} is introduced in order to compensate for missing Doppler broadening in the line extinction profile. The actual correlation of the velocity field with ξ_{turb} is nontrivial due to the nonlocality of the radiation field, which is additionally impeded by nonlinear atomic physics. Therefore, this correlation has to be determined empirically by comparing the results of 3D lineformation calculations with their ⟨3D⟩ counterparts (see also Steffen et al. 2009). We intend to perform such calculations in an upcoming work.
Finally, Chiavassa et al. (2011) showed that using a realistic turbulent pressure contribution to the hydrostatic equilibrium in 1D red supergiant atmospheres, greatly improves the derived surface gravity in these stars. This extra pressure component also leads to an expansion of the atmosphere compared to a 1D model stratification without turbulent pressure. This is referred to as atmospheric levitation (see Trampedach 2001). This will affect pmodes by affording them a larger cavity, and hence lowering their frequencies. This is part of the seismic nearsurface effect which has plagued helio and asteroseismology.
3.2.4. Total pressure and density
Fig. 22 ⟨3D⟩ total pressure p_{tot} against the optical depth log τ_{Ross} for various stellar parameters (solid lines). For comparison we also plot the 1D models with dashed lines. Note the different ordinate scale in the top panel. 
Fig. 23 ⟨3D⟩ stratifications of the local electron number density n_{el} against optical depth log τ_{Ross} for various stellar parameters (solid lines). For comparison, we also show the corresponding stratifications from 1D models (dashed lines). 
The total pressure is defined as the sum of thermodynamic and turbulent pressure, p_{tot} = p_{th} + p_{turb}, and the former consists of gas and radiation pressure, p_{th} = p_{gas} + p_{rad}. In Fig. 22, we show the total pressure for various stellar parameters. In contrast to the previous quantities, p_{tot} decreases with higher T_{eff}, lower log g, and higher metallicity. From the three stellar parameters, the influence of the metallicity is the strongest. We find the highest pressures (and densities) in the coolest metalpoor dwarfs and the lowest pressures in the hottest metalrich giants. In the upper layers of hot metalpoor dwarfs, we find pressures systematically increased with respect to their 1D counterparts, which is accompanied by similar behavior in ρ, p_{gas} and p_{turb}. As we showed above, a significant fraction of the total pressure is contributed by turbulent pressure in the SAR and in the upper layers (see Fig. 21), in particular towards higher T_{eff} and lower log g. Moreover, we note that the temporal and horizontal ⟨3D⟩_{z}averages from our relaxed simulations are very close to hydrostatic equilibrium, and the turbulent pressure contributes significantly to this equilibrium.
Since the mean density stratifications look qualitatively similar to the total pressure ones, we refrain from showing them. Instead we prefer to show, in Fig. 18, the peak density ^{15}ρ_{peak}, which is the density at the height of the maximum rmsvertical velocity (see Fig. 17). The density ρ_{peak} increases with lower T_{eff}, higher log g, and lower [Fe/H]. These variations with stellar parameters arise due to the radiative transfer, since the cooling and heating rates (see Eq. (4)) depend on density ρ and opacity κ. We showed in Sect. 3.1.2 that with higher metallicity and opacity, the hydrostatic stratification is set at lower ρ.
3.2.5. Electron number density
Next, we discuss the properties of the electron number density n_{el} (Fig. 23), which is the temporal and spatial average of the local electron density on layers of constant Rosseland optical depth. The electron number density drops by about ~ 2 dex at the transition from the interior to the photosphere. This is due to the fact that the density itself decreases here, and due to the recombination of hydrogen at the photospheric transition. The convective flux consists to ~ 1/3 of F_{ion} (see Sect. 3.2.8), therefore, as the hot ionized plasma reaches the surface, it radiates away energy, recombines, and overturns into downdrafts, thereby reducing the number of free electrons. The electron density increases with higher T_{eff}, lower log g, and higher [Fe/H]. The electron pressure p_{el} = n_{el}k_{B}T follows similar trends as the electron density in terms of variations with stellar parameters and depth.
3.2.6. Entropy
Fig. 24 ⟨3D⟩ and ⟨3D⟩_{z} mean stratifications (solid and dashed respectively) of the entropy s vs. the total pressure normalized to the pressure at the optical surface log p_{tot}/p_{surf} for various stellar parameters. We show also the sstratifications of the 1D models (dotted lines). Note the different ordinate scale in the top panel. 
Local, boxinastar, 3D RHD atmosphere models have well defined boundary conditions at the bottom boundary because of the adiabaticity of the convection zone, even though they are relatively shallow and comprise only a small fraction of the convection zone. Indeed, the specific entropy per unit mass of the plasma stays constant across most of the convective zone, in particular for the upflows. In Fig. 24, we show the average entropy. Below the optical surface, the entropy converges asymptotically against s_{bot} into deeper layers, especially the averages on constant geometrical height ⟨3D⟩_{z} (dashed lines). As the hot plasma in the granules reaches the optical surface, it becomes transparent, thereby a large fraction of the energy is radiated away. This results in a decrease in entropy, until it reaches a minimum at the top of the convection zone (log τ_{Ross} ~ 0.0). Further up, the entropy then increases again due to the decoupling of the radiation and matter above the photosphere, which results in an almost isothermal atmosphere. The 1D models (dotted lines) exhibit larger entropy stratifications in the deeper convection zone, in particular for higher T_{eff}, thereby overestimating the entropy jump increasingly due to the fixed mixinglength parameter α_{MLT} with 1.5 for all stellar parameters.
3.2.7. Superadiabatic temperature gradient
Fig. 25 ⟨3D⟩_{Ross} (solid lines) and ⟨3D⟩_{z} (dashed lines) superadiabatic gradient ∇_{sad} vs. optical depth log τ_{Ross} for various stellar parameters. Furthermore, for comparison, we also show the corresponding values from 1D models (dotted lines). Note the different ordinate scale in the top panels. 
We limit ourselves to show only the superadiabatic gradient ∇_{sad} = ∇ − ∇_{ad}, since it combines the important properties of both the total and the adiabatic temperature gradient (∇ and ∇_{ad} respectively). In Fig. 25, we show ∇_{sad} averaged over constant geometrical height or Rosseland optical depth. The peak in ∇_{sad}, which looks like a skewed Gaussian function, arises solely from the temperature gradient ∇. The superadiabatic gradient peaks around log τ_{Ross} ~ 1.0 − 2.0, and becomes subadiabatic, i.e. ∇_{sad} < 0.0, above the optical surface at log τ_{Ross} < 0.0 − 0.5. The entropy jump correlates directly with the superadiabatic gradient, since ∇_{sad} = 1/c_{p}[∂s/∂lnp_{tot}] and one can show that ∂s/dz = c_{P}/H_{P}(∇ − ∇_{ad}). Hence, it is no surprise that they exhibit similarity in the peak amplitude and position. In particular, the peak amplitude increases with increasing T_{eff} and log g (see Δs in Fig. 18; a functional fit for is given in Appendix B). The position of on the optical depth scale ⟨3D⟩_{Ross} (triangles), hence the position of the steepest temperature gradient, changes slightly with stellar parameters. However, similar to the position of , in the ε − ρplane, the distribution of is regular (see Fig. 16), namely it shifts systematically towards higher ε and lower ρ with increasing T_{eff}.
As it is clear in Fig. 25, one finds substantial differences in ∇_{sad} when comparing the two ⟨3D⟩ stratifications with their 1D counterparts, namely, the 1D gradients exhibit distinctively larger amplitudes. These differences arise partly due the missing turbulent pressure in the 1D case, but do not resolve the discrepancies. Furthermore, we find an asymmetrically skewed shape towards the optical surface in the 1D gradients, which is shared by the geometrical averages ⟨3D⟩_{z}, but is not the case for the averages on constant Rosseland optical depth ⟨3D⟩_{Ross}. A main reason for the shown differences between ⟨3D⟩_{Ross} and 1D comes from the averaging over layers of constant τ_{Ross}. The underlying ∇_{ad}s are rather insensitive to the deviations between the ⟨3D⟩_{Ross} and 1D stratifications, so the differences arise mainly due to ∇. Between ⟨3D⟩_{z} and 1D the adiabatic gradients differ in the subphotospheric gradient.
3.2.8. Transport of energy
The individual energy fluxes are quantities worthy of further consideration. The energy flux is conserved only on averages of constant geometrical height ⟨3D⟩_{z}, therefore, we show and discuss the latter here. The total energy flux F_{tot} = F_{rad} + F_{conv} + F_{visc} emerges from the photosphere solely in the form of radiative energy flux. The total energy flux is supplied from the convection zone by the convective energy flux, which is the sum of the enthalpy flux (20)(δj_{z} being the horizontal fluctuations of the vertical mass flux; the average vertical mass flux vanishes) carried in the upflows and the kinetic energy flux (21)arising from the downdrafts (see SN98 and Nordlund et al. 2009). Since the mean kinetic energy flux F_{kin} is negative, the positive enthalpy flux F_{enth} is the major component of the convective energy flux F_{conv}. The enthalpy flux in turn consists of the energy fluxes due to ionization , thermal heat and acoustic (sound) waves F_{acous} = ⟨p_{th}v_{z}⟩ − ⟨p_{th}⟩⟨v_{z}⟩. In Fig. 26, we show the energy fluxes F_{rad}, F_{enth}, F_{kin}, F_{ion}, and F_{th} normalized to the total emergent energy flux (for clarity, we refrain from showing F_{visc} and F_{acous}, since their contribution to F_{tot} is very small). We vary one stellar parameter at a time, while the other two are fixed (T_{eff}, log g, and [Fe/H], from top to bottom in Fig. 26, respectively). Just below the optical surface (0.5 < log p_{tot}/p_{surf} < 1.0), both F_{kin} (solid lines) and F_{enth} (dashed lines) increase towards cool metalpoor dwarfs, i.e. lower T_{eff}, [Fe/H] and higher log g, due to higher densities and velocities. The increased reduction of the total flux by F_{kin} (<0.0) is compensated by a simultaneously higher F_{enth} (>1.0). On the other hand, in deeper layers (log p_{tot}/p_{surf} > 1.5), both converge to similar fractions for all stellar parameters (− 0.17 and 1.14 for F_{kin} and F_{enth} respectively). This convergence to very similar values is rather remarkable. The convective motions seem to follow an exact guideline, which might be correlated to the universal filling factor, however, we will study this more carefully in a future work. We remark that in deeper solar simulations ^{16} (20 Mm) that F_{enth} and F_{kin} increase with depth, while their sum remains constant (see Stein et al. 2009).
Fig. 26 Behavior of the normalized energy fluxes F/F_{tot} against the total pressure normalized to the pressure at the optical surface log p_{tot}/p_{surf} as a function of variations in the individual stellar parameters (T_{eff}, log g, and [Fe/H], from top to bottom, respectively). In each panel, the various curves are shown varying one of the parameters while keeping the other two fixed. The individual normalized components of F_{tot} are F_{enth} (dashed), F_{kin} (dotted), F_{ion} (dashtripledotted), F_{th} (long dashed) and F_{rad} (solid). Averages are evaluated at constant geometrical height. 
The majority of the total energy flux F_{tot} in the convection zone is carried in form of ionized hydrogen^{17} with F_{ion} ≃ 0.67, while thermal heat is the second most important component with F_{th} ≃ 0.29. The acoustic energy constitutes only a small fraction with F_{acous} ≃ 0.04. SN98 found similar fractions with F_{kin} ~ − 0.10 to − 0.15, F_{ion} ~ 2/3 and F_{th} ~ 1/3 for the Sun. The F_{ion} and F_{th} fractions, which are the major constituents of the enthalpy flux, undergo a significant change below the surface, as we show in Fig. 26 for models with different stellar parameters. In particular, the fraction of thermal heat F_{th} becomes more significant at the cost of F_{ion} towards cool metalpoor dwarfs. The thermal flux F_{th} reaches a maximum (up to F_{th,max} ≃ 0.5) just below the surface, but eventually converges close to the above mentioned fractions in deeper layers (long dashed lines).
In 1D MLT models, the convective flux is assumed to consist of the enthalpy flux only, F_{conv,1D} (see Appendix A.1). This is a result of the MLT assumption of symmetric flows which means that the kinetic energy fluxes in the up and downflows cancel exactly. As remarked by Henyey et al. (1965), the details on F_{kin}, F_{ion}, F_{th} and F_{acous} are not at hand due to the lack of a selfconsistent velocity field. The energy fluxes from 3D RHD simulations, on the other hand, arise selfconsistently from solving the coupled equations of radiative hydrodynamics, without further assumptions.
As mentioned above, the emergent total energy flux F_{tot} is carried in the convection zone mainly by the positive enthalpy flux F_{enth} (Eq. (20)). Therefore, one can approximate the convective energy flux with the mean jump in enthalpy ^{18} Δh times the mean vertical mass flux of the upflows below the optical surface, hence (22)At the transition region, the enthalpy jump Δh is primarily caused by the strong drop in internal energy ε, hence entropy s, and the thermodynamic pressure work is rather small (note the change of p_{tot} below the surface log τ_{Ross} > 0.0 in Fig. 22), i.e. Δh ≈ TΔs, where T is the temperature at the surface. By approximating T ≃ T_{eff}, one can expect the total energy flux to depend to first order on the mean entropy jump ^{19}, density, and vertical velocity: (23)This approximation can already be retrieved on dimensional grounds, however, we derived the latter in order to explain the systematic variations of Δs, and ρ_{peak} with stellar parameters, which we have observed above (see Figs. 5 and 18, respectively). The emergent radiative energy flux is correlated with Δs, ρ and v_{z}, and the respective composition resulting from the individual contributions varies with stellar parameters.
The interplay between the radiative heating and cooling rates q_{rad} (Eq. (4)) and hydrostatic equilibrium, require a different density stratification for different stellar parameters due to the dependence of opacity on thermodynamic variables, as we showed in Sect. 3.1.2. The resulting density variations will induce adjustments in the vertical velocity and entropy jump. Furthermore, we find with increasing log g or lower [Fe/H] at a fixed T_{eff}, the density increases, which is compensated by lower Δs and v_{z} (Δs ∝ Δρ^{1}, see Eq. (15)). We would like also to emphasize the remarkably important (nonlocal) influence of the rather thin photospheric transition region on basically the whole convection zone, since the entropy deficiency of the turbulent downdrafts are generated mainly here. The latter sets the entropy jump and the convective driving (see Nordlund et al. 2009).
Fig. 27 Normalized cooling and heating rates vs. optical depth log τ_{Ross} for various stellar parameters. The shown averages are retrieved on constant geometrical height ⟨3D⟩_{z}. Note the different ordinate scale in the top panel. 
The radiative heating and cooling rates q_{rad} (Eq. (4)) due to radiative losses enter the hydrodynamic equations as a source and sink term in the energy equation (Eq. (3)). It is the divergence of the radiative flux q_{rad} = −∇·F_{rad}, and a large, negative q_{rad}, the cooling peak, marks the transition of energy transport from fully convective below the optical surface to fully radiative close to the photosphere. To better illustrate the depth dependence of q_{rad} and the comparison among different models, in Fig. 27, we show the normalized cooling and heating rates . One can see that the amplitude of increases with higher T_{eff}, accompanied by an increase in the width of the cooling peak. The position of the maximum absolute amplitude coincides with the position of , since the cooling rate (radiative loss) is setting the entropy fluctuations, hence the superadiabatic gradient (see Sect. 3.2.7). Furthermore, this location moves into upper layers for higher T_{eff} (from log τ_{Ross} ≃ 1.0 up to 0.2 for T_{eff} = 4000 to 7000 K respectively). On the other hand, the width of the photospheric transition region Δ_{ph} = Δlog τ_{Ross}(q_{rad} < 0) clearly widens for hotter T_{eff}, but also, in particular, for metalpoor giants (see top right panel in Fig. 27). While for cool dwarfs the width is typically Δ_{ph} ≈ 3.0 dex, for hot metalpoor giants, it reaches Δ_{ph} ≈ 5.0 dex (see, e.g., model with 5000 K in right top panel of Fig. 27).
3.3. Comparison with 1D models
3.3.1. 1D models
A differential comparison between 1D and 3D in terms of approaches in the modeling of stellar atmospheres is of obvious relevance here. Therefore, we developed a planeparallel, hydrostatic, 1D atmosphere code (hereafter simply referred to as the 1D code) that is based on a similar physical treatment as the MARCS code with a few simplifications (see Appendix A and Gustafsson et al. 2008 for more details). We employ exactly the same EOS and opacities as in the individual 3D models, thereby excluding differences due to dissimilar input physics. Also, we applied the ⟨3D⟩ models as initial stratifications for the 1D models. These mean ⟨3D⟩ stratifications are defined on an equidistant optical depth scale from log τ_{Ross} = −5.0 to + 5.0 in steps of 0.1. The wellresolved optical depth scale reduces discretization errors in the 1D atmosphere calculations, thereby making the 1D⟨3D⟩ comparison more reliable.
In Fig. 28, we show a comparison of the 1D and ⟨3D⟩ temperature stratifications. One can immediately extract that the upper layers of the atmospheres are systematically overestimated in the 1D models by up to ~ 1000 K, in particular for metalpoor stars [Fe/H] ≤ − 2.0 (for solar models the maximal difference is ~ 500 K). In the optically thin layers of 1D models, stable against convection, radiative equilibrium is enforced. However, in the upper layers of the metalpoor ⟨3D⟩ models, the effect of the nonvanishing adiabatic cooling rate is to shift the balance with radiative heating to lower temperatures due to a scarcity and weakness of spectral lines at lower metallicities (Asplund et al. 1999a; Collet et al. 2007). Interesting are also the hotter temperature stratifications for a few giants (T_{eff}/log g = 4500 K/1.5 and 5000 K/2.0) towards higher metallicity ([Fe/H] > − 2.0), which results from the radiative equilibrium at higher temperatures. On the other hand, with the 1D models, we find systematically cooler temperatures below the photosphere log τ_{Ross} ≃ 2.0 with up to ~ 1000 K (here there is no difference with different metallicities). Therefore, one has to keep in mind that, with 1D atmosphere models, and for metalpoor stars in particular, these severe effects on the stratifications can lead to large systematic errors in spectroscopic abundance determinations, up to 0.5 dex or more in logarithmic abundance, depending on the formation height of the individual spectral lines used in the analysis (e.g. Asplund et al. 1999a; Asplund & García Pérez 2001; Collet et al. 2006, 2007; Caffau et al. 2008, 2011; González Hernández et al. 2010; Kučinskas et al. 2013). We will return to this issue using our new grid of 3D stellar models in subsequent investigations.
In the 1D model calculations the mixinglength parameter is kept constant with α_{MLT} = 1.5, which is the commonly applied value (see Gustafsson et al. 2008). However, it is wellknown that α_{MLT} varies with stellar parameters (see Ludwig et al. 1999; Bonaca et al. 2012, Magic et al., in prep.). Therefore, we caution that a single fixed value will lead to severe differences in atmospheric stratification. The systematic deviations beneath the optical surface in the temperature stratification between 1D and ⟨3D⟩ towards cool dwarfs can be interpreted as the manifestation of the wrong α_{MLT} (see Fig. 28). Furthermore, p_{turb} is neglected in the 1D code, which affects the stratification by reducing the gas pressure (see Sect. 3.2.3).
3.3.2. MARCS and ATLAS models
Fig. 28 Comparison of the temperature stratification for the ⟨3D⟩ and the 1D models by showing the difference ⟨1D⟩ − ⟨3D⟩ (colored bars) between log τ_{Ross} = −5.0 and 2.0 in the Kieldiagram. We present four different metallicity ([Fe/H] = −3.0, − 2.0, − 1.0 and 0.0). 
Last, we would also like to briefly compare our ⟨3D⟩ stratifications with the currently widely applied MARCS and ATLAS models (see Fig. 29, we show only the comparison with MARCS modes, since the ATLAS models look qualitatively rather similar). We find qualitatively similar deviations as with the 1D models above. At the same time, here we also have additional differences due to the different input physics (EOS and opacities). The largest differences between the ⟨3D⟩ and 1D MARCS stratifications of metalpoor stellar atmospheres are slightly higher, with ~ 1300 K at [Fe/H] = −3.0, while for solar metallicity the temperatures are underestimated in 1D by ~ 500 K at mosts. Below the surface, the differences amount to ~ 1000 K. The ATLAS models are up to ~ 850 K hotter at the top and ~ 1000 K cooler below the surface. In both cases, the deviations at the top increase towards lower [Fe/H].
4. Conclusions
We presented here a comprehensive grid of realistic, stateoftheart, threedimensional (3D), timedependent, radiativehydrodynamic (RHD) stellar atmosphere models for latetype stars, covering a substantial portion of stellar parameter space, and provided a detailed description of the approach we followed for the construction of models. With the aid of our realistic 3D RHD simulations, we are able to access and render details of stellar atmospheres and subsurface convection, that are out of reach for 1D models and also inaccessible by observations. We presented and discussed a number of important global physical properties of the simulations as well as the mean stratifications resulting from the relatively large amount of data.
The constant entropy value of the adiabatic convection zone has a profound influence on several aspects and properties of the 3D RHD simulations. In particular, we find systematic correlations among the constant entropy value of the adiabatic convection zone, the entropy jump, and the vertical velocity, which we interpreted as scaling relations. In addition, we find that the variation in intensity contrast is enhanced at lower metallicity. Also, we determined that the granule size scales basically with the pressure scale height close to the surface, which can be explained in the picture of what we refer to as Nordlund scaling relation (NSR).
Fig. 29 Similar as Fig. 28, however here we compare the ⟨3D⟩ with the MARCS models for [Fe/H] = −3.0, − 2.0, − 1.0 and 0.0. For a better comparison, we applied the same temperature range as given Fig. 28. 
We discussed in great detail the depthdependent temporal and spatial averages of various important physical quantities. In particular, we determined and examined various systematic trends in the variations of the entropy jump, the density, and the vertical velocity with stellar parameters. The latter can be discussed by regarding the changes in the transition of energy transport from convective to radiative at the photosphere. Namely, for different stellar parameters, the coupling between radiation and matter through the radiative transfer necessitates specific physical conditions due to changes in the opacity, which in turn alters the density. These variations in the density on the other hand require adjustments in the entropy jump and the vertical velocity. This can be illustrated under consideration of the total energy flux and conservation of energy. The named important values are coupled with each other, and these also set basically the general physical framework of the stellar atmosphere. The actual particular connections of these correlations have to be studied carefully in more detail, thus possibly leading to an improved understanding of the physical mechanisms operating in subsurface convection, hence stellar atmospheres.
We compared our 3D models and their mean stratifications to 1D models employing the same input physics, thereby revealing important systematic differences between the two kinds of models due to the incomplete treatment of convection by the 1D mixinglength theory (MLT) and the assumption of radiative equilibrium. The latter leads to an overestimation of the temperature stratification in metalpoor stars. While below the optical surface, we find that the temperatures are typically underestimated due to a fixed mixing length (α_{MLT} = 1.5), in particular for higher T_{eff} and lower log g. Also, we find that MLT fails to render a realistic vertical velocity field. The often neglected turbulent pressure has towards giants a nonnegligible contribution on the total pressure, thereby, indicating that the thermal gas pressure is also overestimated significantly. We also quantified the differences with widely used 1D atmosphere models, in particular ATLAS and MARCS. For a number of important values we provide functional fits with stellar parameters, so that these can be accessed immediately. Thereby, one can easily scale new 3D models based on these informations.
The present work is meant to be an introduction to a series of papers on the Staggergrid. The material discussed here is in fact just a small fraction of the actual information contained in the complete simulation data set. On the other hand, the list of potential applications for 3D models is also long. Because of space constraints, we will discuss in more detail the different methods we have applied for computing temporal and spatial averages of our 3D RHD models in a separate paper. Therein, we will also discuss our routines for the interpolation of ⟨3D⟩ atmosphere models, and explore the differences between ⟨3D⟩ and 3D line formation. We have also compiled details on statistical properties, which we will present individually. These can be later on utilized for a socalled 1.5D line formation (see Ayres et al. 2006). An obvious plan on the agenda is also to analyze carefully the detailed properties of granulation and the intergranular lane for different stellar parameters across the HRdiagram in the future.
As the major purpose of theoretical atmosphere models, we are computing full 3D synthetic spectra with OPTIM3D (see Chiavassa et al. 2009) for all of our models (Chiavassa et al., in prep.). These will be made publicly accessible and can be used for various applications, e.g. spectroscopic parameter determination. Furthermore, we will derive new limbdarkening coefficients based on the full 3D synthetic spectra (Magic et al., in prep.), which are vital for applications such as the characterization of transiting exoplanets (Hayek et al. 2012). We intend to calibrate photometric colors and radial velocities from the spectra. For abundance determinations, we will construct an extensive library of synthetic highresolution spectral lines. These will serve for deriving 3D effects on lineshifts, lineasymmetries and bisectors.
Also an important application is the calibration of free parameters in 1D models by employing our 3D simulations, in particular MLT and micro and macroturbulence. We will incorporate our ⟨3D⟩ models into stellar evolutionary models, which will result in improved stellar structures useful for asteroseismology. And with 3D line formation calculations we will calibrate the microturbulence.
Despite the enormous success and the abinitio nature of 3D atmosphere modeling, as last we want to mention the weaknesses. In order to keep the computation costs reasonable, the radiative transfer is simplified with the opacity binning method, which may influence the outcome. Also, the numerical resolutions of these so called largeeddy simulations are not resolving the microscopic viscous dissipation length scales, hence, the need to introduce numerical diffusion. However, these do not affect the main properties of the macroscopic flows and of the physical stratification. Also, we minimized the diffusion coefficients once under the constraint of numerical stability, and then applied for all the simulations. These issues can be solved, and the easiest rectification will be enhancement of the numerical resolution, in particular for giants. Furthermore, we will conduct improvements in the radiative transfer by employing more bins and possibly more more angles in the future. By inserting magnetic fields with different field strength, we will explore the influence of the latter on convection for different types of stars. Accounting for nonLTE effects is extremely expensive for 3D simulations, therefore, these are usually neglected, however, even these will be included eventually in the future.
Online material
Appendix A: The Staggergrid 1D atmospheres
The following discussion concerns solely 1D atmosphere models and MLT, therefore, similar quantities as discussed above may deviate (e.g. F_{conv}). The numerical code that we used for computing 1D atmospheres for the Staggergrid models solves the coupled equations of hydrostatic equilibrium and energy flux conservation in 1D planeparallel geometry. The 1D models use the same EOS and opacity package in order to allow consistent 3D1D comparisons. The set of equations and numerical methods employed for their solution are similar to those of the MARCS code (Gustafsson et al. 2008) with a few changes and simplifications that will be outlined in the following. The resulting model atmospheres yet maintain very good agreement with MARCS models (see Sect. 3.3.2).
Appendix A.1: Basic equations
Assuming 1D planeparallel geometry with horizontal homogeneity and dominance of hydrostatic equilibrium over all vertical flow simplifies the equation of motion (Eq. (2)) to the hydrostatic equilibrium equation (A.1)where κ_{std} and τ_{std} are a standard opacity and corresponding optical depth (e.g. the Rosseland mean), p_{gas} and p_{turb} denote gas pressure and turbulent pressure, ρ is the gas density, and g is the surface gravity. Radiation pressure is neglected, as in the 3D simulations. Turbulent pressure is estimated using the expression (A.2)with the scaling parameter β that corrects for asymmetries in the velocity distribution and the mean turbulent velocity v_{turb} that is used as a free, independent parameter.
The depthintegral of the energy equation (Eq. (3)) reduces to the flux conservation equation, (A.3)where F_{rad} is the radiative energy flux, F_{conv} is the convective energy flux, σ is the StefanBoltzmann constant and T_{eff} is the stellar effective temperature. Contrary to the 3D case, effective temperature now appears as a boundary value and is thus a free parameter. Owing to numerical instabilities of the formulation, Eq. (A.3) is replaced in the higher atmosphere (τ_{Ross} ≲ 10^{2}) with the radiative equilibrium condition (A.4)where J_{λ} and S_{λ} are the mean intensity and the source function, similar to Eq. (6). In the 3D case, q_{rad} is explicitly calculated and is nonzero in general. Enforcing the condition of radiative equilibrium q_{rad} ≡ 0 in 1D leads to an atmospheric stratification where an exact balance of radiative heating and cooling in each layer is achieved, ignoring the effects of gas motion.
The mean intensity and the radiative energy flux at each depth are obtained by solving the radiative transfer equation, (A.5)where μ = cosθ with the polar angle θ off the vertical, I_{λ} is the specific intensity at wavelength λ, and τ_{λ} is the vertical monochromatic optical depth (with τ_{λ} = 0 above the top of the atmosphere). A Planck source function S_{λ} = B_{λ} is assumed. The monochromatic mean intensity and radiative flux are then delivered by the integrals In the absence of an explicit convection treatment, convective energy transfer is estimated using a variant of the mixing length recipe described in Henyey et al. (1965). The convective flux is given by the expression (A.8)where ρ is the gas density, c_{p} is the specific heat capacity, T is the temperature, and v_{MLT} is the convective velocity. The wellknown free mixing length parameter α_{MLT} = l_{m}/H_{p} sets the distance l_{m} in units of the local pressure scale height H_{p} over which energy is transported convectively. See Gustafsson et al. (2008) for details of the expressions used to obtain the convective velocity v_{MLT} and the factor δΔ = Γ/(1 + Γ)∇_{sad}, which scales superadiabaticity ∇_{sad} = ∇ − ∇_{ad} of the atmospheric stratification (see also Sect. 3.2.7), by a convective efficiency factor with the optical thickness τ_{e} = κ_{Ross}l_{m}. We adopt the same parameters y = 0.076 and ν = 8 as Gustafsson et al. (2008) for the radiative heat loss term and turbulent viscosity that enter the above quantities.
Appendix A.2: Numerical methods
The system of equations is solved using a modified NewtonRaphson method with an initial stratification of temperature T and gas pressure p_{gas} on a fixed Rosseland optical depth grid. Discretized and linearized versions of the hydrostatic equation and the energy flux equation (or radiative equilibrium condition, respectively) provide the inhomogeneous term and the elements of the Jacobian matrix for the system of 2N linear equations, where N is the number of depth layers. The radiation field is computed for each NewtonRaphson iteration using the integral method, based on a secondorder discretization of the fundamental solution of the radiative transfer equation (Eq. (A.5)).
The corrections ΔT and Δp_{gas} derived from the system of linear equations are multiplied by a variable factor < 1 that is automatically regulated by the code to aid convergence. Convergence is assumed when the (relative) residuals of the 2N equations decrease beneath a preset threshold. Note that, contrary to the 3D simulations, the effective temperature is now an adjustable parameter; the requirement of minimal residuals automatically leads to an atmospheric stratification with correct T_{eff} through the energy flux equation.
In order to obtain a 1D model, a given ⟨3D⟩ stratification provides the initial input for the NewtonRaphson iterations, along with the targeted effective temperature and surface gravity. The same EOS tables that were used for the 3D simulation provide gas density, specific heat capacity, and adiabatic gradient as a function of T and p_{gas}. Likewise, the tables containing group mean opacities and the Rosseland mean opacity provide the required microphysics for solving the radiative transfer equation, ensuring maximal consistency with the 3D simulations.
Once convergence has been achieved for the 1D stratification, the mixing length parameter α_{MLT} can be calibrated to obtain a better approximation to the ⟨3D⟩ stratification in the convection zone beneath the stellar surface.
Appendix B: Functional fits
The resulting amount of data from our numerical simulations is enormous. A convenient way to provide important key values is in form of functional fits, which can be easily utilized elsewhere (e.g. for analytical considerations). In the present paper we have frequently discussed various important global properties that are reduced to scalars. Some of them are global scalar values and
some are determined at a specific height from the ⟨3D⟩ stratifications, i.e. temporal and spatial averages on layers of constant Rosseland optical depth. We fitted these scalars with stellar parameters for individual suitable functions, thereby enforcing a smooth rendering. However, we would like to warn against extrapolating these fits outside their range of validity, i.e. outside the confines of our grid. Also, one should consider that possible small irregularities between the grid steps might be neglected, which arise due to nonlinear response of the EOS and the opacity. On the other hand, we provide also most of the actual shown values in Table C.1.
We use three different functional bases for our fits and we perform the leastsquares minimization with an automated LevenbergMarquardt method. Instead of the actual stellar parameters, we employ the following transformed coordinates: x = (T_{eff} − 5777)/1000, y = log g − 4.44 and z = [Fe/H]. Furthermore, we find that in order to accomplish an optimal fit with three independent variables, f_{i}(x,y,z), simultaneously, the metallicity should be best included implicitly as nested functions in the form of second degree polynomial , each resulting in three independent coefficients a_{i}. The linear function (B.1)is applied for the following quantities: s_{min} (Fig. 5), log ρ_{peak} (Fig. 18), (Fig. 18), log d_{gran} (Fig. 10), log Δt_{gran} and . The resulting coefficients are given in Table B.1. On the other hand, we considered the exponential function (B.2)for s_{bot}, Δs (Figs. 5 and 18) and [p_{turb}/p_{tot}]_{peak} (Fig. 21). For ∇_{peak} and 15 we applied the following function (B.3)with coefficients for f_{2} and f_{3} listed in Table B.2. Finally, we showed in Fig. 7 the entropy jump Δs as a function of s_{bot}, which we fitted with (B.4)The resulting coefficients are listed in Table B.3.
Appendix C: Tables
In Table C.1 we have listed important global properties with the stellar parameters. Due to the lack of space, we show only an excerpt with solar metallicity ([Fe/H] = 0.0). The complete table is available at CDS http://cds.ustrasbg.fr.
Stellar parameters: effective temperature, surface gravity (Cols. 1 and 2 in [K] and [dex]).
The actual values we used are n_{1} = 0.005 and n_{2} = 0.8 (see Eq. (9) in Kritsuk et al. 2011).
In 1D MLT modeling the term convective efficiency is commonly referred to the mixinglength. The latter is in 3D RHD simulations referred to the mass mixinglength, which is the inverse gradient of the mass flux (see Trampedach & Stein 2011).
For example, Δh can be determined at the top and bottom of the photospheric transition region (see Fig. 16).
Acknowledgments
We acknowledge access to computing facilities at Rechen Zentrum Garching (RZG) through MaxPlanck Institute for Astrophysics (MPA) and the National Computational Infrastructure (NCI) of Australia, through Australian National University (ANU), where the simulations were carried out. We are grateful to W. Däppen for access to the code and data tables for the EOS. And we thank B. Plez and B. Edvardsson for providing the MARCS lineopacities. Also, we acknowledge the Action de Recherche Concertée (ARC) grant provided by the Direction générale de l’Enseignement non obligatoire et de la Recherche scientifique – Direction de la Recherche scientifique – Communauté francaise de Belgique, and the F.R.S. FNRS. Remo Collet is the recipient of an Australian Research Council Discovery Early Career Researcher Award (project number DE120102940). We thank the referee for the helpful comments.
References
 Allende Prieto, C., García López, R. J., Lambert, D. L., & Ruiz Cobo, B. 2000, ApJ, 528, 885 [NASA ADS] [CrossRef] [Google Scholar]
 Allen de Prieto, C., Lambert, D. L., & Asplund, M. 2001, ApJ, 556, L63 [NASA ADS] [CrossRef] [Google Scholar]
 Asplund, M., & García Pérez, A. E. 2001, A&A, 372, 601 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Asplund, M., Nordlund, Å., & Trampedach, R. 1999a, 173, 221 [Google Scholar]
 Asplund, M., Nordlund, Å., Trampedach, R., & Stein, R. F. 1999b, A&A, 346, L17 [NASA ADS] [Google Scholar]
 Asplund, M., Nordlund, Å., Trampedach, R., Allen de Prieto, C., & Stein, R. F. 2000, A&A, 359, 729 [NASA ADS] [Google Scholar]
 Asplund, M., Grevesse, N., & Sauval, A. J. 2005, 336, 25 [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Ayres, T. R., Plymate, C., & Keller, C. U. 2006, ApJS, 165, 618 [NASA ADS] [CrossRef] [Google Scholar]
 Badnell, N. R., Bautista, M. A., Butler, K., et al. 2005, MNRAS, 360, 458 [NASA ADS] [CrossRef] [Google Scholar]
 Barklem, P. S., Stempels, H. C., Allen de Prieto, C., et al. 2002, A&A, 385, 951 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beeck, B., Collet, R., Steffen, M., et al. 2012, A&A, 539, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bergemann, M., Lind, K., Collet, R., Magic, Z., & Asplund, M. 2012, MNRAS, 427, 27 [NASA ADS] [CrossRef] [Google Scholar]
 BöhmVitense, E. 1958, ZAp, 46, 108 [NASA ADS] [Google Scholar]
 Bonaca, A., Tanner, J. D., Basu, S., et al. 2012, ApJ, 755, L12 [NASA ADS] [CrossRef] [Google Scholar]
 Caffau, E., Ludwig, H.G., Steffen, M., et al. 2008, A&A, 488, 1031 [CrossRef] [EDP Sciences] [Google Scholar]
 Caffau, E., Ludwig, H.G., Steffen, M., Freytag, B., & Bonifacio, P. 2011, Sol. Phys., 268, 255 [NASA ADS] [CrossRef] [Google Scholar]
 Canuto, V. M., & Mazzitelli, I. 1991, ApJ, 370, 295 [NASA ADS] [CrossRef] [Google Scholar]
 Carlsson, M., Stein, R. F., Nordlund, Å., & Scharmer, G. B. 2004, ApJ, 610, L137 [NASA ADS] [CrossRef] [Google Scholar]
 Cassisi, S., Salaris, M., Castelli, F., & Pietrinferni, A. 2004, ApJ, 616, 498 [NASA ADS] [CrossRef] [Google Scholar]
 Castelli, F., & Kurucz, R. L. 2004 [arXiv:astroph/0405087] [Google Scholar]
 Castelli, F., Gratton, R. G., & Kurucz, R. L. 1997, A&A, 318, 841 [NASA ADS] [Google Scholar]
 Chiavassa, A., Plez, B., Josselin, E., & Freytag, B. 2009, A&A, 506, 1351 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Chiavassa, A., Collet, R., Casagrande, L., & Asplund, M. 2010, A&A, 524, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chiavassa, A., Freytag, B., Masseron, T., & Plez, B. 2011, A&A, 535, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chiavassa, A., Bigot, L., Kervella, P., et al. 2012, A&A, 540, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Collet, R., Asplund, M., & Trampedach, R. 2006, 306 [Google Scholar]
 Collet, R., Asplund, M., & Trampedach, R. 2007, A&A, 469, 687 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Collet, R., Nordlund, Å., Asplund, M., Hayek, W., & Trampedach, R. 2009, Mem. Soc. Astron. It., 80, 719 [NASA ADS] [Google Scholar]
 Collet, R., Hayek, W., Asplund, M., et al. 2011, A&A, 528, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de Laverny, P., RecioBlanco, A., Worley, C. C., & Plez, B. 2012, A&A, 544, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Di Mauro, M. P., ChristensenDalsgaard, J., RabelloSoares, M. C., & Basu, S. 2002, A&A, 384, 666 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Feautrier, P. 1964, Comptes Rendus Académie des Sciences (série non spécifiée), 258, 3189 [NASA ADS] [Google Scholar]
 Frebel, A., Collet, R., Eriksson, K., Christlieb, N., & Aoki, W. 2008, ApJ, 684, 588 [NASA ADS] [CrossRef] [Google Scholar]
 Freytag, B., Steffen, M., Ludwig, H.G., et al. 2012, J. Comput. Phys., 231, 919 [Google Scholar]
 Fuhrmann, K., Axer, M., & Gehren, T. 1993, A&A, 271, 451 [NASA ADS] [Google Scholar]
 González Hernández, J. I., Bonifacio, P., Ludwig, H.G., et al. 2010, A&A, 519, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Grupp, F. 2004, A&A, 420, 289 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gudiksen, B. V., Carlsson, M., Hansteen, V. H., et al. 2011, A&A, 531, A154 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gustafsson, B., Bell, R. A., Eriksson, K., & Nordlund, A. 1975, A&A, 42, 407 [NASA ADS] [Google Scholar]
 Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hauschildt, P. H., Allard, F., & Baron, E. 1999, ApJ, 512, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Hayek, W., Asplund, M., Carlsson, M., et al. 2010, A&A, 517, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hayek, W., Sing, D., Pont, F., & Asplund, M. 2012, A&A, 539, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Henyey, L., Vardya, M. S., & Bodenheimer, P. 1965, ApJ, 142, 841 [NASA ADS] [CrossRef] [Google Scholar]
 Holweger, H., & Mueller, E. A. 1974, Sol. Phys., 39, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Johnson, H. R., & Krupp, B. M. 1976, ApJ, 206, 201 [NASA ADS] [CrossRef] [Google Scholar]
 Kritsuk, A. G., Nordlund, Å., Collins, D., et al. 2011, ApJ, 737, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Kurucz, R. L. 1979, ApJS, 40, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Kurucz, R. L. 1993, SYNTHE spectrum synthesis programs and line data (Kurucz CDROM, Cambridge, MA: Smithsonian Astrophysical Observatory) [Google Scholar]
 Kučinskas, A., Steffen, M., Ludwig, H.G., et al. 2013, A&A, 549, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ludwig, H.G. 2006, A&A, 445, 661 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ludwig, H.G., & Kučinskas, A. 2012, A&A, 547, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ludwig, H.G., Freytag, B., & Steffen, M. 1999, A&A, 346, 111 [NASA ADS] [Google Scholar]
 Ludwig, H.G., Caffau, E., Steffen, M., et al. 2009a, Mem. Soc. Astron. Italiana, 80, 711 [Google Scholar]
 Ludwig, H.G., Samadi, R., Steffen, M., et al. 2009b, A&A, 506, 167 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ludwig, H.G., Caffau, E., Steffen, M., et al. 2010, 265, 201 [Google Scholar]
 Magic, Z., Serenelli, A., Weiss, A., & Chaboyer, B. 2010, ApJ, 718, 1378 [NASA ADS] [CrossRef] [Google Scholar]
 Maltby, P., Avrett, E. H., Carlsson, M., et al. 1986, ApJ, 306, 284 [Google Scholar]
 Mihalas, D. 1970 [Google Scholar]
 Mihalas, D., Dappen, W., & Hummer, D. G. 1988, ApJ, 331, 815 [NASA ADS] [CrossRef] [Google Scholar]
 Muthsam, H. J., Kupka, F., LöwBaselli, B., et al. 2010, New A, 15, 460 [NASA ADS] [CrossRef] [Google Scholar]
 Nahar, S. N. 2004, Phys. Rev. A, 69, 042714 [NASA ADS] [CrossRef] [Google Scholar]
 Nissen, P. E., Primas, F., Asplund, M., & Lambert, D. L. 2002, A&A, 390, 235 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nordlund, A. 1982, A&A, 107, 1 [NASA ADS] [Google Scholar]
 Nordlund, A., & Dravins, D. 1990, A&A, 228, 155 [NASA ADS] [Google Scholar]
 Nordlund, Å., & Stein, R. F. 2001, ApJ, 546, 576 [Google Scholar]
 Nordlund, Å., Stein, R. F., & Asplund, M. 2009, Liv. Rev. Sol. Phys., 6, 2 [Google Scholar]
 Pereira, T. M. D., Asplund, M., & Kiselman, D. 2009a, Mem. Soc. Astron. It., 80, 650 [NASA ADS] [Google Scholar]
 Pereira, T. M. D., Kiselman, D., & Asplund, M. 2009b, A&A, 507, 417 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Piskunov, N. E., Kupka, F., Ryabchikova, T. A., Weiss, W. W., & Jeffery, C. S. 1995, A&AS, 112, 525 [NASA ADS] [Google Scholar]
 Plez, B. 2008, Phys. Scr. T, 133, 014003 [Google Scholar]
 Rosenthal, C. S., ChristensenDalsgaard, J., Nordlund, Å., Stein, R. F., & Trampedach, R. 1999, A&A, 351, 689 [NASA ADS] [Google Scholar]
 Roudier, T., & Muller, R. 1986, Sol. Phys., 107, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Ruchti, G. R., Fulbright, J. P., Wyse, R. F. G., et al. 2010, ApJ, 721, L92 [NASA ADS] [CrossRef] [Google Scholar]
 Ruiz Cobo, B., & del Toro Iniesta, J. C. 1992, ApJ, 398, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Skartlien, R. 2000, ApJ, 536, 465 [NASA ADS] [CrossRef] [Google Scholar]
 Smalley, B., Gardiner, R. B., Kupka, F., & Bessell, M. S. 2002, A&A, 395, 601 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sneden, C. A. 1973, The University of Texas at Austin, Ph.D. Thesis [Google Scholar]
 SocasNavarro, H. 2011, A&A, 529, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Steffen, M. 1993, in Inside the Stars, eds. W. W. Weiss, & A. Baglin, IAU Colloq. 137, ASP Conf. Ser., 40, 300 [Google Scholar]
 Steffen, M., Ludwig, H.G., & Kruess, A. 1989, A&A, 213, 371 [NASA ADS] [Google Scholar]
 Steffen, M., Ludwig, H.G., & Caffau, E. 2009, Mem. Soc. Astron. Italiana, 80, 731 [NASA ADS] [Google Scholar]
 Stein, R. F., & Nordlund, A. 1998, ApJ, 499, 914 [Google Scholar]
 Stein, R. F., & Nordlund, Å. 2001, ApJ, 546, 585 [NASA ADS] [CrossRef] [Google Scholar]
 Stein, R. F., Georgobiani, D., Schafenberger, W., Nordlund, Å., & Benson, D. 2009, 1094, 764 [Google Scholar]
 Stein, R. F., Lagerfjärd, A., Nordlund, Å., & Georgobiani, D. 2011, Sol. Phys., 268, 271 [NASA ADS] [CrossRef] [Google Scholar]
 Stempels, H. C., Piskunov, N., & Barklem, P. S. 2001, 223, 878 [Google Scholar]
 Tanner, J. D., Basu, S., & Demarque, P. 2013, ApJ, 767, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Trampedach, R. 2001, J. Astron. Data, 7, 8 [Google Scholar]
 Trampedach, R. 2007, AIP Conf. Proc. 948, 141 [Google Scholar]
 Trampedach, R., & Stein, R. F. 2011, ApJ, 731, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Trampedach, R., Stein, R. F., ChristensenDalsgaard, J., & Nordlund, Å. 1999, 173, 233 [Google Scholar]
 Trampedach, R., Asplund, M., Collet, R., Nordlund, Å., & Stein, R. F. 2013, ApJ, 769, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Vernazza, J. E., Avrett, E. H., & Loeser, R. 1976, ApJS, 30, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weiss, A., & Schlattl, H. 2008, Ap&SS, 316, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Wende, S., Reiners, A., & Ludwig, H.G. 2009, A&A, 508, 1429 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Williamson, J. H. 1980, J. Comput. Phys., 35, 48 [Google Scholar]
All Tables
Stellar parameters: effective temperature, surface gravity (Cols. 1 and 2 in [K] and [dex]).
All Figures
Fig. 1 Kiel diagram (T_{eff} − log g diagram) showing the targeted Staggergrid parameters for the 217 models, comprising seven different metallicities (colored circles). Four additional standard stars (see text) are also indicated (squares). In the background, the evolutionary tracks for stellar masses from 0.7 to 1.5 M_{⊙} and for solar metallicity are shown (thin grey lines). 

In the text 
Fig. 2 Top panel: nonequidistant vertical spacing Δz of the depth scale as a function of geometrical depth in our solar model (solid line). The zscale is optimized to resolve the flows and thermal structure, which naturally results in the highest spatial resolution around the photosphere. Furthermore, we also show the maximum of the vertical gradient of the absorption coefficient max⟨dlnα_{Ross}/dz⟩ as a function of depth (dotted line). Bottom panel: aspect ratio Δz/Δx (solid line) and we also indicated its lower allowed limit with 1:4 (dashed line). The actual verticaltohorizontal aspect ratio ranges from 0.26 at the photosphere to 1.18 at the bottom of the simulation domain. 

In the text 
Fig. 3 Twelve opacity bins selected for the solar simulation by plotting the opacity strength (or, more precisely, the formation height) against wavelength for all sampled wavelength points. The individual bin elements are indicated by colored symbols. For clarity, we plotted only a subset of the wavelength points considered for the opacity binning procedure. In the background, we included the smoothed histogram of the opacity strength distribution (blue contour). This shows how the majority of λpoints are mostly concentrated close to the continuumforming layers and only a smaller fraction contributes to lines. 

In the text 
Fig. 4 Comparison of the radiative heating and cooling resulting from monochromatic computations q_{λ} (filled dots) and the opacity binning method q_{bin} (solid line) for the solar model mean stratification. In the top panel we show both q_{rad} vs. optical depth, while in the bottom panel, we compare the two against each other. 

In the text 
Fig. 5 Overview of the constant entropy value of the adiabatic convection zone, which is given by the fixed entropy of the inflowing plasma at the bottom, s_{bot}, (circles) as well as the atmospheric entropy minimum, s_{min}, (squares) for two metallicities ([Fe/H] = −0.5 and 0.0, blue and black respectively) against T_{eff}. The jump in entropy, Δs, is indicated with vertical dotted lines. Simulations with the same gravity are connected with functional fits for s_{bot} and s_{min} (solid and dashed lines respectively; see Appendix B), while similar simulations with different [Fe/H] are connected with short solid black lines. 

In the text 
Fig. 6 Lines of constant entropy of the adiabat s_{bot} from 1.3 to 4.0 × 10^{13}erg/g/K in steps of 0.3 × 10^{13}erg/g/K and for [Fe/H] = −3.0 and 0.0 (blue solid and black dashed lines respectively) in the Kiel diagram. Additionally, we show evolutionary tracks for 0.7 to 1.5 M_{⊙} with solar metallicity (thin grey lines). 

In the text 
Fig. 7 Comparison of the entropy jump Δs against the constant entropy value of the adiabatic convection zone s_{bot} for all grid models (filled circles; colorcoding explained in the legend). Furthermore, we show also hyperbolic tangent functional fits (see Eq. (B.4) and Appendix B) between T_{eff} = 4000 and 6000 K (red lines; top panel) and between log g = 1.5 and 4.5 (grey lines; bottom panel). 

In the text 
Fig. 8 Overview of the emergent (bolometric) intensity for a selection of stars, namely mainsequence (MS), turnoff (TO), Kgiant and Kdwarf (from left to right, respectively) at a given time instant. For each star, we show four metallicities [Fe/H] = 0.0, − 1.0, − 2.0 and − 3.0 (from top to bottom, respectively). To facilitate comparisons between the different metallicity of each star, the intensity scale and the horizontal geometrical size of the metalpoor simulations are identical to [Fe/H] = 0.0, and the individual intensity contrasts [%] are indicated in each box. 

In the text 
Fig. 9 Temporally averaged histograms of the (bolometric) intensity (solid lines) for various stellar parameters (the binsize is 100). To ease comparisons between the different stellar parameters, we normalized the individual intensity scales with its mean value. Furthermore, we have indicated the respective FWHM of the intensity histograms, I_{FWHM}, (dotted lines), a measure of the intensity contrast. Note the different ordinate scale in the top panel. The bimodal distributions are given due to the asymmetry and the large contrast in the up and downflows. 

In the text 
Fig. 10 Overview of the (bolometric) intensity contrast ΔI_{rms} against T_{eff} for [Fe/H] = −2.0 and 0.0 (blue and black respectively). Models with the same gravity, but different T_{eff} are connected. The enhanced naked or hidden granulation at lower metallicity can be extracted in the larger range of intensity contrasts (see text for details). 

In the text 
Fig. 11 Overview of the granule diameter d_{gran} derived from the maximum of the mean 2D spatial power spectrum of the bolometric intensity against T_{eff} for [Fe/H] = −2.0 and 0.0 (blue and black respectively). Models with the same gravity are connected by their respective functional fits (solid lines; see Appendix B). 

In the text 
Fig. 12 Top panel: comparison of the granule size approximated with the pressure scale height d_{NSR} (see Eq. (18)) against the mean granule diameter d_{gran} (same as Fig. 11) for all models. The individual stellar parameters are indicated (T_{eff}, log g and [Fe/H] with red, gray, and blue respectively). We indicated the position of d_{NSR} = d_{gran} the solid diagonal line. Bottom panel: relative residuals. We indicated also the mean residual (dotted line), which amounts to ~ 30%. Here, models with the same gravity are connected (solid grey lines). 

In the text 
Fig. 13 ⟨3D⟩ stratifications of the temperature vs. optical depth for various stellar parameters (solid lines) and 1D models with α_{MLT} = 1.5 (dashed lines). Furthermore, we have marked the positions of entropy minimum (plus), vertical peak velocity (diamond), and maximum in ∇_{sad} = ∇ − ∇_{ad} (triangle). 

In the text 
Fig. 14 ⟨3D⟩ temperature stratifications with the variation of one stellar parameter at a time, while the two others are fixed (T_{eff}, log g, and [Fe/H], from top to bottom, respectively). We show our standard averages on constant optical depth ⟨3D⟩ (solid line) and on constant geometrical depth ⟨3D⟩_{z} (dotted line). 

In the text 
Fig. 15 Overview of the maximum temperature gradient ∇_{peak} (top panel) and Rosseland opacity κ_{Ross} taken at the height τ_{Ross} ≈ 3.0 (bottom panel) against T_{eff} for [Fe/H] = −2.0 and 0.0 (blue and black respectively). Models with the same surface gravity are connected by their respective functional fits in the top panel (solid lines; see Appendix B). 

In the text 
Fig. 16 Mean internal energy against mean density for dwarf models (log g = 4.5) with [Fe/H] = 0.0 and − 2.0 (left and right respectively). The specific isocontours for the entropy s_{bot} (green solid) and s_{min} (green dotted), Rosseland opacity per volume, ρκ_{Ross}, (blue) and temperature T (orange) are underlayed. Moreover, the positions of entropy minimum s_{min} (plus), optical surface (large square), vertical peak velocity (diamonds) and maximum in ∇_{sad} (triangle) are marked respectively. The amplitude of is indicated with horizontal bars with markings in 1 km s^{1}. We also included the 1D models with α_{MLT} = 1.5 (blue dashed lines). The range in optical depth is shown from log τ_{Ross} = −5.0 to + 5.0 for each dex (small squares). However, we note that our simulations boxes are much deeper (⟨log τ_{Ross}⟩ ≈ + 7.5). 

In the text 
Fig. 17 Vertical rmsvelocity, v_{z,rms}, from the ⟨3D⟩ stratifications (solid lines) and convective velocity v_{MLT} from the corresponding 1D MLT models (dashed lines) as function of optical depth log τ_{Ross} for various stellar parameters. 

In the text 
Fig. 18 Overview of the entropy jump (top panel), maximal vertical rmsvelocity below the surface (middle panel) and the density at the same height ρ_{peak} (bottom panel) vs. T_{eff} for [Fe/H] = −2.0 and 0.0 (blue and black respectively). Models with the same gravity are connected with their respective functional fits (solid lines; see Appendix B). Note the inverse correlation between density and velocity. 

In the text 
Fig. 19 Comparison of the entropy jump Δs (top panel) and the density at the optical depth of the maximal vertical velocity below the optical surface ρ_{peak} (bottom panel) vs. maximal vertical rmsvelocity for all models. A linear fit is also indicated in both panels (red lines). 

In the text 
Fig. 20 ⟨3D⟩ stratification of the absolute value of the vorticity ω vs. optical depth is shown for various stellar parameters. 

In the text 
Fig. 21 Fraction of turbulent pressure to total pressure p_{turb}/p_{tot} vs. optical depth log τ_{Ross} is displayed for various stellar parameters. 

In the text 
Fig. 22 ⟨3D⟩ total pressure p_{tot} against the optical depth log τ_{Ross} for various stellar parameters (solid lines). For comparison we also plot the 1D models with dashed lines. Note the different ordinate scale in the top panel. 

In the text 
Fig. 23 ⟨3D⟩ stratifications of the local electron number density n_{el} against optical depth log τ_{Ross} for various stellar parameters (solid lines). For comparison, we also show the corresponding stratifications from 1D models (dashed lines). 

In the text 
Fig. 24 ⟨3D⟩ and ⟨3D⟩_{z} mean stratifications (solid and dashed respectively) of the entropy s vs. the total pressure normalized to the pressure at the optical surface log p_{tot}/p_{surf} for various stellar parameters. We show also the sstratifications of the 1D models (dotted lines). Note the different ordinate scale in the top panel. 

In the text 
Fig. 25 ⟨3D⟩_{Ross} (solid lines) and ⟨3D⟩_{z} (dashed lines) superadiabatic gradient ∇_{sad} vs. optical depth log τ_{Ross} for various stellar parameters. Furthermore, for comparison, we also show the corresponding values from 1D models (dotted lines). Note the different ordinate scale in the top panels. 

In the text 
Fig. 26 Behavior of the normalized energy fluxes F/F_{tot} against the total pressure normalized to the pressure at the optical surface log p_{tot}/p_{surf} as a function of variations in the individual stellar parameters (T_{eff}, log g, and [Fe/H], from top to bottom, respectively). In each panel, the various curves are shown varying one of the parameters while keeping the other two fixed. The individual normalized components of F_{tot} are F_{enth} (dashed), F_{kin} (dotted), F_{ion} (dashtripledotted), F_{th} (long dashed) and F_{rad} (solid). Averages are evaluated at constant geometrical height. 

In the text 
Fig. 27 Normalized cooling and heating rates vs. optical depth log τ_{Ross} for various stellar parameters. The shown averages are retrieved on constant geometrical height ⟨3D⟩_{z}. Note the different ordinate scale in the top panel. 

In the text 
Fig. 28 Comparison of the temperature stratification for the ⟨3D⟩ and the 1D models by showing the difference ⟨1D⟩ − ⟨3D⟩ (colored bars) between log τ_{Ross} = −5.0 and 2.0 in the Kieldiagram. We present four different metallicity ([Fe/H] = −3.0, − 2.0, − 1.0 and 0.0). 

In the text 
Fig. 29 Similar as Fig. 28, however here we compare the ⟨3D⟩ with the MARCS models for [Fe/H] = −3.0, − 2.0, − 1.0 and 0.0. For a better comparison, we applied the same temperature range as given Fig. 28. 

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