© 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 semi-empirical 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 late-type 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 line-blanketed atmosphere models for late-type stars appeared with the publication of MARCS (Gustafsson et al. 1975, 2008) and ATLAS models (Kurucz 1979; Castelli & Kurucz 2004). Subsequently, other one-dimensional (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 mixing-length theory (MLT, see Böhm-Vitense 1958), which is characterized by several free parameters, the most commonly known being the mixing-length 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 mixing-length theory has in total four free parameters (see Böhm-Vitense 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 plane-parallel 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⁵ 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 semi-empirical models. In these models, the temperature stratification is inferred from observations (e.g. from lines forming at different heights or continuum center-to-limb variations). Often-used semi-empirical 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 three-dimensional (3D) atmosphere structures using inversion techniques (Ruiz Cobo & del Toro Iniesta 1992; Socas-Navarro 2011). Semi-empirical 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, time-dependent, non-local, 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 center-to-limb 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 macro-turbulence that have hampered stellar spectroscopy for many decades. Instead, in 3D simulations, convection emerges naturally, by solving the time-dependent hydrodynamic equations for mass-, momentum- and energy-conservation, coupled with the 3D radiative transfer equation in order to account correctly for the interaction between the radiation field and the plasma. Also, the non-thermal 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 self-organize naturally to form a distinct flow pattern that exhibits the characteristic granulation at the surface. Furthermore, additional spectral observables such as limb-darkening 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 metal-poor late-type 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 re-absorption of the continuum-radiation from below and adiabatic cooling due to the expansion of upflowing gas. In metal-poor 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 metal-poor star HE1327-2326, 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 stellar-structure models, and which are usually referred to as surface effects (Rosenthal et al. 1999). With classical 1D stellar structures, higher frequency p-modes 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 p-mode 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 micro-variability 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 planet-hosting 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 Stagger-code, which is developed specifically to run efficiently on the massively parallel machines available today (Nordlund & Galsgaard 1995 ¹; Kritsuk et al. 2011). The Bifrost-code is an Oslo derivative of the Stagger-code (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⁵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⁵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 time-scale. 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 high-resolution spectroscopic and asteroseismic observations.

In this paper, we present a new large grid of 3D model atmospheres for late-type 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 post-process 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 ². 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 Stagger-grid.

2. Methods

2.1. The STAGGER-code

The 3D model atmospheres presented here were constructed with a custom version of the Stagger-code, a state-of-the-art, multipurpose, radiative-magnetohydrodynamics (R-MHD) code originally developed by Nordlund & Galsgaard (1995), and continuously improved over the years by its user community. In pure radiation-hydrodynamics mode, the Stagger-code solves the time-dependent 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 ³, 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 Stagger-code uses a sixth-order explicit finite-difference scheme for numerical derivatives and the corresponding fifth-order interpolation scheme. The solution of the hydrodynamic equations is advanced in time using an explicit third-order Runge-Kutta integration method (Williamson 1980). The code operates on a staggered, Eulerian, rectangular mesh: the thermodynamic variables, density and internal energy per volume, are cell-centered, 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 flux-conservative formulation of the magnetohydrodynamic equations, at the same time ensuring that the magnetic field remains divergence-free. The solution of the discretized equations is stabilized by hyper-viscosity which aims at minimizing the impact of numerical diffusion on the simulated flow, while providing the necessary diffusion for large-eddy simulations with finite-difference schemes (see also Nordlund & Galsgaard 1995, for further details). The values of the numerical viscosity parameters ⁴ are empirically tuned for the solar surface-convection 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 Stagger-code we used for this work is fully MPI-parallel. 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 so-called box-in-a-star 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 box-in-a-star setup is still valid because the bulk up-flows 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³, since a resolution of about 200³ − 250³ was found to be adequate (see Asplund et al. 2000). Five layers at the bottom and the top are reserved for the so-called ghost-zones: these extra layers serve to enforce boundary conditions for the high-order 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 well-resolved 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₂ 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 OPAL-EOS. 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 bi-cubic 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)⁵ 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 ⁶ (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 VALD-2 database ⁷ (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 (domain-decomposed) 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 multi-group 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 time-series 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 free-streaming limits, i.e. between the optical thick and optical thin regimes, below and above the photospheric transition zone, respectively. The mean bin-opacity κ_i is calculated as a Rosseland-like average (7)in the optical thick regime, and as a mean-intensity-weighted 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 streaming-regime”), 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 STAGGER-grid models

The Stagger-grid 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 ⁸. 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 non-solar metallicity analogs, and four additional standard stars, namely HD 84937, HD 140283, HD 122563 and G 64-12 that are presented in Bergemann et al. (2012). For metal-poor chemical compositions with [Fe/H] ≤  − 1.0 we applied an α-enhancement of [α/Fe] =  + 0.4   dex, in order to account for the enrichment by core-collapse 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 main-sequence (MS) over the turnoff (TO) up to the red-giant branch (RGB) for low-mass 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 work-flow of procedures for generating our grid. More specifically, we developed a large set of IDL-tools 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 π ₀ of the radial p-mode with the largest amplitude, then damp these modes with an artificial exponential-friction term with period π ₀ 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.

• Re-compute the opacity tables with 12 bins for the relaxed simulation.

• Evolve the simulations for at least ~7 periods of the fundamental p-mode, roughly corresponding to ~ 2 convective turnover times, typically, a few thousand time-steps, 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, p-mode 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 quasi-stationary equilibrium state. Also, when the resulting effective temperature of an otherwise relaxed simulation deviated more than 100 K from the targeted T_eff, we re-scaled the simulation to the targeted value of T_eff and started over from the top of our list of relaxation steps.

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 upwards-shifts 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 depth-dependent 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 depth-scale 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).

After the initial scaling, we construct the geometrical depth scale z for the new simulation by enforcing the same (quasi-)hydrostatic-equilibrium condition as in the starting simulation, but with the newly scaled pressure and density. The resulting new z-scale is usually not smooth, therefore we generate a new z-scale, 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 z-scale 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 optical-depth scale. With such an optimized z-scale we can efficiently resolve the same features with fewer grid-points, compared to an equidistant vertical mesh. Furthermore, a limiting vertical-to-horizontal 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 zero-point 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 sound-crossing 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 sound-crossing 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 opacity-binning approximation is to reproduce the radiative heating and cooling rates as accurately as possible with a small number of opacity-bins, 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⁵ wavelength points in the opacity-sampling (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 L-shaped 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 non-linear 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.

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.

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 Stefan-Boltzmann law T_eff = [F_rad/σ]^1/4, with σ being the Stefan-Boltzmann 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_effs, 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

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 ⁹ the first law of thermodynamics in the form (14)adding an arbitrary integration constant in order to shift the zero-point 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 ¹⁰ (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 ¹¹ 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 entropy-jump (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 non-trivial, 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 non-local radiative transfer depends non-linearly 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).

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 metal-poor stars feature a broader range in s_bot compared to metal-rich 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 ¹², 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 iso-s 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 quasi-exponentially 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öhm-Vitense 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.

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 color-coded the T_eff- and log g-values 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 mixing-length 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

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 main-sequence (MS) simulation (the Sun), a turnoff (TO) simulation, a K-giant, and a K-dwarf 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 metal-poor 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¹⁰, see SN98), which is the by far the dominant continuum opacity source in the visible for late-type stellar photospheres. Since the temperature difference between the granules and the intergranular lanes is very large (>10³  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 self-similarity of the granulation patterns despite the large variations in size-scales. 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.

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, low-contrast dwarfs, to hotter, high-contrast 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 spectro-interferometric 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

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 metal-rich ones.

We find a remarkable validation for the approximation of the maximal horizontal extent of a granule based on mass-conservation 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²] ρ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 rms-velocity 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 scale-height 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 sub-solar metallicity ([Fe/H] = 0.0 and − 2.0, solar and metal-poor, 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 non-adiabatic 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 late-type 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

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_effs, 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 metal-poor ⟨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 solar-metallicity models is largely controlled by radiative equilibrium, while for low-metallicity models this is not generally the case: for metal-poor models, the absorption features become considerably weaker, therefore, the radiative heating by spectral line re-absorption (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, metal-poor 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 non-linear, and we find no simple systematic trends within our grid models.

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_effs 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_effs, 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_effs 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.

In the sub-photospheric region (log τ_Ross> 0.5), where convection dominates, the temperature gradients ∇ = dlnT/dlnp_tot¹³ 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 metal-poor 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 metal-poor simulations have flatter and hot metal-poor simulations have steeper temperature stratifications than the metal-rich 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 back-warming.

• 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 line-blanketing. This in turn leads to a steepening of the temperature gradient and to additional heating of the sub-surface layers, also known as back-warming (see Mihalas 1970; Nordlund et al. 2009). In 1D models it is straightforward to quantify the effect of back-warming, 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, line-blanketing and back-warming effects are also naturally included through our opacity-binning method. Isolating the radiative back-warming 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 enthalpy-jump Δh in Eq. (22)). The ε-jump coincides with the sub-photospheric 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 bound-free 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 metal-poor 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 metal-poor ⟨3D⟩ models also display higher densities at the optical surface, thereby spanning a larger ρ-range for different T_effs. 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 metal-poor dwarfs than for the solar metallicity dwarfs. Therefore, the density ranges covered above the optical surface by the individual atmospheres is small for metal-poor models (min[log ρ] ~  − 1.0 and − 2.0, respectively; see Fig. 16).

3.2.2. Velocity field

Next, we consider the velocity field in our simulations, which arise self-consistently from the solution of the hydrodynamic equations. In Fig. 17 we show the rms-velocity 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 sub-photospheric 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 p-modes, excited in the SAR but leaking out of the acoustic cavity as they have frequencies above the acoustic cut-off. The metal poor simulations show a slightly smaller increase in v_z,rms, since their density gradients are shallower due to steeper T-gradients. 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 metal-poor simulations display a flatter profile, and the position of the minimum is increasingly shifted towards higher layers, especially for extreme metal-poor 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 non-local 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.

The peak vertical rms-velocity, (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 irradiation-duration, 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 non-trivial and beyond the scope of the present paper.

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 power-law 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 log-log scale (red lines in Fig. 19), exhibiting the slopes of and , hence a scaling with the respective slopes.

In 3D RHD simulations, the non-thermal, macroscopic velocity fields arising from convective instabilities are computed self-consistently 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 free-parameter-dependent velocity field v_MLT is derived for the convective flux in MLT. Also, for radiative transfer and spectral line formation calculations, two ad-hoc Gaussian velocity distributions – the so-called micro- and macroturbulence (ξ_turb and χ_turb, respectively) – are usually introduced to model Doppler broadening of spectral lines associated with non-thermal (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 depth-independent value of the microturbulence ξ_turb and one global value of the macroturbulence χ_turb are applied in theoretical spectrum syntheses with 1D model atmospheres. Full-3D line formation calculations using 3D models similar to those described here, have demonstrated that in late-type stars the required non-thermal 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 non-thermal velocity field is clearly depth-dependent, while micro- and macroturbulence are almost always assumed to be non-varying 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öhm-Vitense 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.

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 tube-like 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

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_effs.

For hotter stars, in particular metal-rich 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 density-weighted 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 ¹⁴ξ_turb, since v_z,rms is a horizontal average of the velocity field. In 1D line formation calculations, a depth-independent ξ_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 non-locality of the radiation field, which is additionally impeded by non-linear atomic physics. Therefore, this correlation has to be determined empirically by comparing the results of 3D line-formation 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 p-modes by affording them a larger cavity, and hence lowering their frequencies. This is part of the seismic near-surface effect which has plagued helio- and asteroseismology.

3.2.4. Total pressure and density

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 metal-poor dwarfs and the lowest pressures in the hottest metal-rich giants. In the upper layers of hot metal-poor 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 ¹⁵ρ_peak, which is the density at the height of the maximum rms-vertical 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_elk_BT follows similar trends as the electron density in terms of variations with stellar parameters and depth.

3.2.6. Entropy

Local, box-in-a-star, 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 mixing-length parameter α_MLT with 1.5 for all stellar parameters.

3.2.7. Superadiabatic temperature gradient

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 sub-adiabatic, 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 ∇_ads 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 sub-photospheric 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_thv_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 metal-poor 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 ¹⁶ (20   Mm) that F_enth and F_kin increase with depth, while their sum remains constant (see Stein et al. 2009).

The majority of the total energy flux F_tot in the convection zone is carried in form of ionized hydrogen¹⁷ 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 metal-poor 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 self-consistent velocity field. The energy fluxes from 3D RHD simulations, on the other hand, arise self-consistently 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 ¹⁸ Δ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 ¹⁹, 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 (non-local) 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).

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 metal-poor giants (see top right panel in Fig. 27). While for cool dwarfs the width is typically Δ_ph ≈ 3.0   dex, for hot metal-poor 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 plane-parallel, 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 well-resolved 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 metal-poor 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 metal-poor ⟨3D⟩ models, the effect of the non-vanishing 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 metal-poor 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 mixing-length parameter is kept constant with α_MLT = 1.5, which is the commonly applied value (see Gustafsson et al. 2008). However, it is well-known 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

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 metal-poor 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, state-of-the-art, three-dimensional (3D), time-dependent, radiative-hydrodynamic (RHD) stellar atmosphere models for late-type 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).

We discussed in great detail the depth-dependent 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 mixing-length theory (MLT) and the assumption of radiative equilibrium. The latter leads to an overestimation of the temperature stratification in metal-poor 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 non-negligible 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 Stagger-grid. 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 so-called 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 HR-diagram 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 limb-darkening 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 high-resolution spectral lines. These will serve for deriving 3D effects on line-shifts, line-asymmetries 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 ab-initio 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 large-eddy 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 non-LTE 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 Stagger-grid 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 Stagger-grid models solves the coupled equations of hydrostatic equilibrium and energy flux conservation in 1D plane-parallel geometry. The 1D models use the same EOS and opacity package in order to allow consistent 3D-1D 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 plane-parallel 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 depth-integral 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 Stefan-Boltzmann 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 well-known 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 super-adiabaticity ∇_sad = ∇ − ∇_ad of the atmospheric stratification (see also Sect. 3.2.7), by a convective efficiency factor with the optical thickness τ_e = κ_Rossl_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 Newton-Raphson 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 Newton-Raphson iteration using the integral method, based on a second-order 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 Newton-Raphson 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 non-linear 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 least-squares minimization with an automated Levenberg-Marquardt 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₂ and f₃ 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.u-strasbg.fr.

---

Acknowledgments

We acknowledge access to computing facilities at Rechen Zentrum Garching (RZG) through Max-Planck 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 line-opacities. 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

All Tables

All Figures