Issue 
A&A
Volume 571, November 2014



Article Number  A68  
Number of page(s)  16  
Section  The Sun  
DOI  https://doi.org/10.1051/00046361/201322461  
Published online  13 November 2014 
Nearpolytropic stellar simulations with a radiative surface^{⋆}
^{1} Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden
^{2} Department of Astronomy, AlbaNova University Center, Stockholm University, 10691 Stockholm, Sweden
^{3} MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
email: barekat@mps.mpg.de
Received: 7 August 2013
Accepted: 30 May 2014
Context. Studies of solar and stellar convection often employ simple polytropic setups using the diffusion approximation instead of solving the proper radiative transfer equation. This allows one to control separately the polytropic index of the hydrostatic reference solution, the temperature contrast between top and bottom, and the Rayleigh and Péclet numbers.
Aims. Here we extend such studies by including radiative transfer in the gray approximation using a Kramerslike opacity with freely adjustable coefficients. We study the properties of such models and compare them with results from the diffusion approximation.
Methods. We use the Pencil code, which is a highorder finite difference code where radiation is treated using the method of long characteristics. The source function is given by the Planck function. The opacity is written as κ = κ_{0}ρ^{a}T^{b}, where a = 1 in most cases, b is varied from −3.5 to + 5, and κ_{0} is varied by four orders of magnitude. We adopt a perfect monatomic gas. We consider sets of onedimensional models and perform a comparison with the diffusion approximation in one and twodimensional models.
Results. Except for the case where b = 5, we find onedimensional hydrostatic equilibria with a nearly polytropic stratification and a polytropic index close to n = (3 − b)/(1 + a), covering both convectively stable (n> 3/2) and unstable (n< 3/2) cases. For b = 3 and a = −1, the value of n is undefined a priori and the actual value of n depends then on the depth of the domain. For large values of κ_{0}, the thermal adjustment time becomes long, the Péclet and Rayleigh numbers become large, and the temperature contrast increases and is thus no longer an independent input parameter, unless the StefanBoltzmann constant is considered adjustable.
Conclusions. Proper radiative transfer with Kramerslike opacities provides a useful tool for studying stratified layers with a radiative surface in ways that are more physical than what is possible with polytropic models using the diffusion approximation.
Key words: radiative transfer / hydrodynamics / Sun: atmosphere
Appendices are available in electronic form at http://www.aanda.org
© ESO, 2014
1. Introduction
Convection in stars and accretion disks is a consequence of radiative cooling at the surface. Pioneering work by Nordlund (1982, 1985) has shown that realistic simulations of solar granulation can be performed with not too much extra effort and the required computing resources are comparable to the mandatory costs for solving the hydrodynamics part. Yet, many studies of hydrodynamic and hydromagnetic convection today ignore the effects of proper radiative transfer, sometimes even at the expense of using computeintensive implicit solvers to cope with a computationally stiff problem in the upper layers where the radiative conductivity becomes large (e.g., Cattaneo et al. 1991; Gastine & Dintrans 2008). Therefore, the main reason for ignoring radiation cannot be just the extra effort, but it is more likely a reduced flexibility in that one is confined to a single physical realization of a system and the difficulty in varying parameters that are in principle fixed by the physics. With only a few exceptions (e.g., Edwards 1990), radiation hydrodynamics simulations of stratified convection also employ realistic opacities combined with a realistic equation of state. In the case of the Sun this means that one can only simulate for the duration of a few days solar time (Stein & Nordlund 1989, 1998, 2012).
There are other types of realistic simulations that are able to cover longer time scales by simulating only deeper layers, so they ignore radiation. However, these simulations still need to pose an upper boundary condition, where the gas is cooled (Miesch et al. 2000). This leads to a granulationlike pattern at a depth where the flow topology is known to consist of individual downdrafts rather than a connected network of intergranular lanes. This compromises the realism of such simulations. Other types of simulations give up the ambition for realism altogether and try to model a “toy Sun” in which the broad range of time and length scales is compressed to a much narrower range (Käpylä et al. 2013). This can be useful if one wants to understand the physics of the solar dynamo, where we are not even sure about the possible importance of the surface (Brandenburg 2005), or the physics of sunspots, where so far only models of a toy Sun have produced spontaneous magnetic flux concentrations similar to those of sunspots (Brandenburg et al. 2013). It is therefore important to know how to manipulate the parameters to accommodate the relevant physics, given certain numerical constraints such as the number of mesh points available.
In the present paper we include radiation, which introduces the StefanBoltzmann constant, σ_{SB}, as a new characteristic quantity into the problem. It characterizes the strength of surface cooling, or, conversely, the temperature needed to radiate the flux that is transported through the rest of the domain. Earlier simulations that ignored radiation have specified the surface temperature in an ad hoc manner so as to achieve a certain temperature contrast across the domain. An example are the simulations of Brandenburg et al. (1996), who specified a parameter ξ as the ratio of pressure scale height at the surface, which is proportional to the temperature at the top, and the thickness of the convectively unstable layer. Alternatively, one can use a radiative surface boundary condition. It involves σ_{SB} and couples therefore the surface temperature T_{top} to the lower part of the system, so T_{top} is then no longer a free parameter, unless one chooses an effective value of σ_{SB} so as to achieve the desired temperature contrast. This was done in recent simulations by Käpylä et al. (2012), who kept the aforementioned parameter ξ as the basic control parameter, which then determines the effective value of σ_{SB} in their simulations.
The goal of the present work is to explore the physics of models that introduce radiation without being confined to just one realization. We do this by using a Kramerslike opacity law, but with freely adjustable parameters. It turns out that it is possible in some cases to imitate polytropic models with any desired polytropic index and Rayleigh number. This then eliminates any restrictions to a single setup, allowing one to perform parameter surveys, just like with earlier polytropic models. To compare radiative transfer models with those in the diffusion approximation, we consider twodimensional convection simulations. An ultimate application of this work is to study the formation of surface magnetic flux concentrations through the negative effective magnetic pressure instability (Brandenburg et al. 2013), which has been shown to produce bipolar region (Warnecke et al. 2013; Mitra et al. 2014), and to investigate the relation to the magnetic cooling instability of Kitchatinov & Mazur (2000), which could favor sunspot formation in the presence of radiative cooling. This will be discussed again at the end of the paper.
We would like to point out that, in view of more general applications, we cannot assume the effective temperature to be given or fixed. Thus, unlike the case usually considered in the theory of stellar atmospheres, the dependence of temperature on optical depth is not known a priori. Therefore, it is more convenient to fix instead the temperature at the bottom of the domain and obtain the effective temperature, and thus the flux, as a result of the calculation.
We begin by presenting first the governing equations and then describe the basic setup of our model. Next we compare a set of onedimensional simulations with the associated polytropic indices that correspond to Schwarzschild stable or unstable solutions. Finally, we explore the effect of including radiative transfer instead of using the diffusion approximation combined with a radiative boundary condition by comparing one and twodimensional simulations.
2. The model
2.1. Governing equations
We solve the hydrodynamics equations for logarithmic density lnρ, velocity u, and specific entropy s, in the form where p is the gas pressure, g is the gravitational acceleration, ν is the viscosity, is the traceless rateofstrain tensor, I is the unit tensor, T is the temperature, and is the radiative flux. For the equation of state, we assume a perfect gas with p = (ℛ /μ)Tρ, where ℛ is the universal gas constant and μ is the mean molecular weight. The pressure is related to s via p = ρ^{γ}exp(s/c_{v}), where the adiabatic index γ = c_{p}/c_{v} is the ratio of specific heats at constant pressure and constant volume, respectively, and c_{p} − c_{v} = ℛ /μ. To obtain the radiative flux, we adopt the gray approximation, ignore scattering, and assume that the source function S (not to be confused with the rateofstrain tensor S) is given by the frequencyintegrated Planck function, so S = (σ_{SB}/π)T^{4}, where σ_{SB} is the StefanBoltzmann constant. The divergence of the radiative flux is then given by (4)where κ is the opacity per unit mass (assumed independent of frequency) and is the frequencyintegrated specific intensity corresponding to the energy that is carried by radiation per unit area, per unit time, in the direction , through a solid angle dΩ. We obtain by solving the radiative transfer equation, (5)along a set of rays in different directions using the method of long characteristics.
2.2. Opacity
For our work it is essential that we can control the value and functional form of the opacity. We therefore choose a Kramerslike opacity given by (6)where a and b are free parameters that characterize the relevant radiative processes. It is useful to consider the radiative conductivity K(ρ,T), which is given by (7)We note that, in a planeparallel polytropic atmosphere, T(z) varies linearly with height z and in the stationary state, K(ρ,T) is constant in the optically thick part. This implies that ρ is proportional to T^{n}, where (8)is the polytropic index (not to be confused with the direction of a ray ). This relation was also used by Edwards (1990), but the author regarded those solutions as “a little contrived”. This is perhaps the case if such solutions are applied throughout the entire domain. It should also be noted that Edwards (1990) included thermal conduction along with radiative transfer. This meant that one had to pose a boundary condition for the temperature at the top also, which will not be necessary in our case, where, unless stated otherwise, no thermal conductivity is included. Indeed, as we shall show, with a Kramerslike opacity, nearly polytropic solutions are a natural outcome in the lower optically thick part of the domain, while in the upper optically thin part of the domain the stratification tends to become approximately isothermal.
For a perfect gas, the specific entropy gradient is related to the gradients of the other thermodynamic variables via (9)and vanishes when n = 1/(γ − 1). For a monatomic gas where γ = 5/3, the stratification is Schwarzschildstable for n> 3/2.
2.3. Boundary conditions
We consider a slab with boundary conditions in the z direction at z_{bot} and z_{top}, where we assume the gas to be stressfree, i.e., (10)We assume zero incoming intensity at the top, and compute the incoming intensity at the bottom from a quadratic Taylor expansion of the source function, which implies that the diffusion approximation is obeyed; see Appendix A of Heinemann et al. (2006) for details. To ensure steady conditions, we fix temperature at the bottom, (11)while the temperature at the top is allowed to evolve freely. There is no boundary condition on the density, but since no mass is flowing in or out, the volumeaveraged density is automatically constant. Since most of the mass resides near the bottom, the density there will not change drastically and will be close to its initial value at the bottom.
2.4. The radiation module
We use for all simulations the Pencil code^{1}, which solves the hydrodynamic differential equations with a highorder finitedifference scheme. The radiation module was implemented by Heinemann et al. (2006). It solves the transfer equation in the form (12)where dτ = κρ dl is the differential of the optical depth along a given ray and l is a coordinate along this ray.
The code is parallelized by splitting the calculation into parts that are local and nonlocal with respect to each processor. There are two local parts that are computeintensive and one that is nonlocal and fast, so it does not require any computation. Since S is assumed independent of I (scattering is ignored), we can write the solution of Eq. (12) as an integral for I(τ), which is thus split into two parts, (13)where the subscripts “extr” and “intr” indicate respectively an extrinsic, nonlocal contribution and an intrinsic, local one. An analogous calculation is done for calculating τ along the geometric coordinate as , where l_{0} is the geometric end point on the previous processor. In the first step, we calculate I_{intr}(τ), which can be evaluated immediately on all processors in parallel, while the first integral is written in the form I_{extr}(τ) = I_{0} e^{τ0 − τ}, where I_{0} and τ_{0} are already being computed as part of the I_{intr} calculation on neighboring processors and the results included in the last step of the computation.
In the second step, the values of and I_{0} = I(τ_{0}) are communicated from the end point of each ray on the previous processor, which cannot be done in parallel, but this does not require any computational time. In the final step one computes (14)and constructs the final intensity as I(τ) = I_{extr}(τ) + I_{intr}(τ).
Instead of solving the radiative transfer equation directly for the intensity, the contribution to the cooling term Q(τ) = I(τ) − S(τ) is calculated instead, as was done also by Nordlund (1982). This avoids roundoff errors in the optically thick part. For further details regarding the implementation we refer to Heinemann et al. (2006). To avoid interpolation, the rays are chosen such that they go through mesh points. The angular integration in Eq. (4) is discretized as (15)where i enumerates the N rays with directions and D/ 3 is a correction factor that is relevant when the number of dimensions, D, of the calculation is less than three. It does not affect the steady state, but it affects the cooling rate both in the optically thick and thin regimes; see Appendix A for details. In one dimension with D = 1, we have N = 2 rays, which are , while for D = 2 we can either have N = 4 with and , or N = 8 with the additional 4 combinations . In three dimensions, the correction factor is D/ 3 = 1, so the angular integral is just 4π times the average of the intensity over all directions.
2.5. Parameters and initial conditions
In the following, we measure length in Mm, speed in km s^{1}, density in g cm^{3}, and temperature in K.
Units used in this paper and conversion into cgs units.
This implies that time, for example, is measured in ks (=1000 s). The advantage of using this system of units is that it avoids extremely large or small values of various quantities by using units that are commonly used in solar physics such as Mm and km/s. A summary of our units and the conversion of various quantities between cgs and our units is given in Table 1.
For the gravitational acceleration, we take g = (0,0, − g) with g = 274 km^{2} s^{2} Mm^{1} being the solar surface value (Stix 2002). Instead of prescribing T_{bot}, we prescribe the sound speed c_{s}, where , and fix c_{s} = c_{s0} = 30 km s^{1} at z_{bot} = 0. With ℛ = 8.314 × 10^{7} erg K^{1} mol^{1} and μ = 0.6 g mol^{1}, this choice corresponds to T_{bot} = 38 968 K. We found it instructive to start with an isothermal solution that is in hydrostatic equilibrium, but not in thermal equilibrium, so the upper parts will gradually cool until a static solution is reached. Thus, we use ρ = ρ_{0}exp( − z/H_{p}), where H_{p} = ℛT/μg is the pressure scale height, and ρ_{0} is a constant that we set to ρ_{0} = 4 × 10^{4} g cm^{3}. This value was chosen based on values from a solar model at a depth of approximately 7 Mm below the surface. However, this particular choice is quite uncritical and just corresponds to renormalizing the opacity. In other words, instead of making a calculation with a ten times larger value of ρ, we can just use an otherwise equivalent calculation with a ten times larger value of κ.
2.6. Simulation strategy
We choose the exponents a and b such that they correspond to five different values of n. In the case a = −1, b = 3, we have K(ρ,T) = const, but the value of n is undefined. Our choice of parameters is summarized in Table 2.
Summary of used a and b.
It is convenient to express κ in the form (16)where is a rescaled opacity and is related to κ_{0} by ; where T_{0} = T_{bot} is used. (By contrast, ρ_{0} is only approximately equal to the density at the bottom – except initially.) With this choice, the units of are independent of a and b, and always Mm^{1} cm^{3} g^{1} (=10^{8} cm^{2} g^{1}). For each value of n, we choose 4 different values of , 10^{5}, 10^{6}, and 10^{7} Mm^{1} cm^{3} g^{1}. We note that the actual Kramers opacity for free–free and bound–free transitions with a = 1 and b = −7/2 has κ_{0} between 6.6 × 10^{22} and 4.5 × 10^{24} cm^{5} g^{2}K^{7/2}, respectively (Kippenhahn & Weigert 1990). This corresponds to and , which are respectively four and six orders of magnitude larger than the largest value considered in this paper.
3. Results
We perform onedimensional simulations with a resolution of 512 equally spaced grid points using five sets of values for the exponents a and b in the expression for the Kramers opacity; see Eq. (16). Each set of runs is denoted by a letter A–E. In the first four sets of runs, we keep a = 1 and change the value of b from −7/2, to 0, 1, and 5. For each of these sets, we perform four runs that differ only in the values of . The numeral on the label of each run refers to a different value of . In Set A, we use a = 1 and b = −7/2. Runs A4, A5, A6 and A7 correspond to equal to 10^{4}, 10^{5}, 10^{6} and 10^{7} Mm^{1} cm^{3} g^{1}, respectively. All the other designations follow the same scheme. All runs have been started with the same isothermal initial condition. However, the size of the domain is changed so as to accommodate the upper isothermal part by a good margin. If the domain is too big, one needs a large number of meshpoints to resolve the resulting strong stratification, and if it is too small, the solution changes in the top part, as will be discussed in Sect. 3.9.
Fig. 1 Vertical temperature profile at five different times t = 0, 3, 30, 120, and 1578 ks for Run A6 with . Squares, circles and crosses represent different optical depths τ = 0.1, τ = 1 and τ = 10, respectively. The arrow represents the direction of the time evolution of the temperature profile. 
After a sufficient amount of running time, a unique equilibrium state is reached and the resulting profiles of temperature, density and entropy have a nearly polytropic stratification in the lower part of the domain and a nearly isothermal stratification in the upper part of the domain. An exception are the runs of Set E where the polytropic index is undefined (n = 0/0). This will be discussed in more detail in Sect. 3.8. We summarize the important quantities obtained from all runs in Table 3. These quantities are calculated in the equilibrium state. All runs show a similar evolution of density, temperature and entropy. In the next sections we describe the resulting profiles in more detail.
Summary of the runs.
3.1. Approach toward the final state
As mentioned above, we find it convenient to obtain equilibrium solutions by starting from an isothermal state. The upper layers begin to cool fastest, and eventually an equilibrium state is reached. We plot the evolution of the temperature profile of Run A6 in Fig. 1 as an exemplary case with . Already after a short time of t = 3 ks (1 h), the temperature has decreased by more than half its initial value at the top and follows a polytropic solution in most of the domain, where the temperature gradient has a similar value than in the equilibrium state. At t = 30 ks (8 h), close to the top boundary, an isothermal part is seen to emerge. However, it takes more than t = 1500 ks (17 days) until the equilibrium solution is reached with a prominent isothermal part of T ≈ 7000 K. The locations of three different optical depths, τ = 0.1, 1 and 10, are shown in Fig. 1. Here, (17)is the optical depth with respect to the surface of the domain. If the domain is tall enough, one can see that an initially isothermal stratification cools down first near the location where τ(z) ≈ 1, which is where the cooling is most efficient. As an example we plot in Fig. 2 the early stages of the temperature evolution at t = 0, 0.01, 0.1, and 0.2 ks for a taller variant of Run A6 with z_{top} = 12 Mm using 1024 equally spaced grid points.
Fig. 2 Temperature over height for Run A6 at different times t = 0, 0.01, 0.1, and 0.2 ks, plotted as solid, dotted, dashed and dasheddotted lines, respectively. The red dots correspond to the location of τ ≈ 1. 
At t = 0.01 ks, the temperature starts to decrease at the height where τ ≈ 1, while in the upper part, which is far enough from the surface, the temperature remains at first unchanged. Only at a somewhat later time (t = 0.1 ks) does the temperature at z = 12 Mm start to decrease. This is explained by the fact that the radiative cooling rate (or inverse cooling time) is largest near τ = 1 (Spiegel 1957; Unno & Spiegel 1966; Edwards 1990); see also Appendix A.
3.2. Temperature stratification
For all runs, the temperature reaches an equilibrium state after a certain time; see Fig. 1. The temperature profile can be divided into two distinguishable parts, a nearly polytropic part which starts from the bottom of the domain and extends to a certain height, and a nearly isothermal part which starts from this height and extends to the top of the domain. The transition of the temperature from the initial state to the equilibrium state follows a specific pattern, which is the same for all the runs. The higher the value of , the lower the temperature is in the isothermal part and the longer it takes to reach this state. Increasing the normalized opacity by three orders of magnitude results in a decrease in the temperature by a factor of five for Set A and a factor of three for Set D. As the exponent b changes from the smallest value in Set A to the largest one in Set D, the slope of the temperature decreases with height. This means that the polytropic part of the atmosphere is smaller for larger values of b. We note that the size of the domain is chosen larger for smaller values of b. For Sets A, B and C, the temperature in the polytropic part is almost the same for different values of , although for the lowest value of the temperature deviates somewhat. However, in Set D, for different values of , the slope of the temperature is different for each value of . This has to do with the fact that in this case with (b − 3)/(1 + a) = −1 the stratification is no longer a polytropic one. (A polytrope with n = −1 would have constant pressure, which is unphysical.)
The temperatures in the isothermal part also show a dependency on b. For the surface temperature in Run A4 is T ≈ 2.2 × 10^{4}K, whereas in Run D4 the value is T ≈ 2.9 × 10^{4}K. A similar behavior can also be seen for the other values of .
Fig. 3 Density, temperature and entropy of the equilibrium state versus height, from left to right, for Sets A, B, C and D, from top to bottom. The four different lines in each plot corresponds to the value of the rescaled opacity , 10^{5}, 10^{6}, 10^{7} Mm^{1} cm^{3} g^{1}. The dots in all plots represent the surface τ ≈ 1. The red dashed lines represent the initial profile of each set. 
Next, we calculate the optical depth for all runs. We find that the transition point from the polytropic part to the isothermal part coincides with the τ ≈ 1 surface. We indicate the surface τ ≈ 1 by dots in all plots in Fig. 3. The polytropic part corresponds to the optically thick part with τ> 1 and the isothermal part corresponds to the optically thin part with τ< 1. For each set, the transition point depends on the value of . As we go from smaller to larger values of , the surface shifts to larger heights and becomes cooler. This is because the radiative heat conductivity K is inversely proportional to and directly proportional to the flux. Therefore, by increasing the value of , K decreases and, as a consequence, the radiative flux also decreases. By decreasing the flux, the effective temperature decreases as . This means that the temperature at the surface is smaller for larger values of . For Set A, the τ ≈ 1 surfaces lie on the polytropic part of the temperature profile. However, by increasing the value of b, the locations of the τ = 1 surfaces shift toward the lower boundary and the optically thick part becomes narrower. This is particularly severe for the solutions with a small value of , especially for , when the boundary condition T = T_{bot} at z = 0 becomes unphysical and the temperature drops between the first two meshpoints in a discontinuous fashion; see Fig. 3.
3.3. Entropy stratification
We plot the entropy profiles for all sets of runs in the equilibrium state in the last column of Fig. 3. For Runs C6−7 and D5−7, the entropy decreases in the polytropic part and starts to increase in the isothermal part. All runs show a positive vertical entropy gradient in the isothermal part. In the lower part, the entropy gradient is positive (∇_{z}s> 0) for Set A, while for Set B it is constant and equal to zero (∇_{z}s ≈ 0). This shows that for Set B, the atmospheres are isentropic. In Sets C and D, except for the case , the entropy gradient is negative, ∇_{z}s< 0. This means that their atmospheres are convectively unstable. (Convection will of course not occur in our onedimensional model, but we will obtain the socalled hydrostatic reference solution that is used to compute the Rayleigh number, as will also be done later in this paper.) In Set D the entropy gradients are larger than in case C where their atmospheres are marginally stable. In the isothermal part of Set C, the entropy gradient is much larger than in Set D. For each set of runs, as we go from smaller values of to larger ones, the entropy profiles have larger gradients.
3.4. Incoming and outgoing intensity profiles
It is instructive to inspect the vertical profiles of the intensity for rays pointing in the up and downward directions, , denoted in the following by I^{±}. If we have just these two rays, the energy flux is given by F_{rad} = (2π/ 3)(I^{+} − I^{−}). In thermal equilibrium, the difference between I^{+} and I^{−} must be constant. This is indeed the case; see Fig. 4, where we plot I^{+} and I^{−} and compare with S. The vertical lines in this figure represents the difference between I^{+} and I^{−}, where I^{+} − I^{−} ≈ 10^{4} erg cm^{2} s^{1} in the whole domain as we have radiative equilibrium . Radiative equilibrium also demands that J = S, so S has to be equal to the average of I^{+} and I^{−}, which is indeed the case.
Fig. 4 Source function and vertical profiles of incoming and outgoing intensity near the surface for Run A7. The dashed line represents the source function and the solid lines represent the incoming intensity I^{+} (red) and outgoing intensity I^{−} (blue). The vertical lines represent the (constant) difference between I^{+} and I^{−}. 
3.5. Radiative heat conductivity
In Sets A, B, and C, the value of radiative heat conductivity K turns out to be constant in the optically thick part of the atmosphere, but not for Set D. The value of K at the bottom of the optically thick part of the domain is denoted in Table 3 by K_{bot}, and agrees roughly with (18)Indeed, for a fixed value of , it has even the same order of magnitude for Sets A, B and C, independently of the value of b, but it is one order of magnitude larger for Set D. This is because in this set the density is lower in the optically thick part compared to the other sets. Moreover, as we go from higher values of to the lower ones, the radiative heat conductivity increases. This can be explained by the inverse proportionality of K with opacity as . For smaller values of , K is larger and vice versa. As an example, we plot in Fig. 5 the resulting vertical profiles of radiative heat conductivity for Set C. We note that K is constant in the optically thick part and starts to increase in the optically thin part. In the optically thin part, κρ decreases, so K increases as K ∝ 1 /κρ. To maintain , the modulus of ∇T has to decrease. As K increases even further, a thermostatic equilibrium can be satisfied if ∇T comes close to zero.
Fig. 5 Radiative heat conductivity K versus height for Set C. K is plotted for different values of where is shown by dotteddashed line, dashed line, dotted line and solid line. 
Fig. 6 Effective temperature T_{eff} versus rescaled opacity for Sets A, C, and D. The crosses, circles and stars show the values of T_{eff} for different values of for Sets A, C, and D, respectively. Different lines correspond to line fit of T_{eff} with normalized opacity . 
3.6. Effective temperature
The effective temperature T_{eff} of all runs is calculated from the z component of the radiative flux , (19)The values of T_{eff} of all sets of runs are summarized in Table 3. By increasing the value of b, T_{eff} also increases. The value of T_{eff} decreases as we go from lower to higher opacities for each set. We plot T_{eff} versus in Fig. 6 for Sets A, C and D where the values of T_{eff} are represented by crosses, circles and stars, respectively. For each set of runs, we fit a line to T_{eff} versus . We find that T_{eff} has a power law relation with . The power of , which is the slope of the plot, depends on the polytropic index and therefore on b. For larger values of b, the power is smaller than for smaller values of b. Additionally, the offset shows also a weak dependence on b. A power law relation between T_{eff} and the opacity of roughly 1/4 can be expected, because of the linear relation of the radiative flux and the opacity. Toward larger b, this dependency is no longer accurate. We also calculate for each run the corresponding optical depth where T = T_{eff}. For all runs, T_{eff} corresponds to the optical depth τ ≈ 1/3. This is less than what is expected for a gray atmosphere, where T_{eff} = T at τ ≈ 2/3. This is presumably because in our integration of τ we have not included the contribution between ∞ and z = z_{top}.
3.7. Thermal adjustment time
In our simulations we define a thermal adjustment time τ_{adjust} as the time it takes for each run to reach its numerically obtained final equilibrium temperature in the isothermal part to within 1% (see Fig. 7). The unit of τ_{adjust} is ks (see Sect. 2.5). The value of τ_{adjust} for all runs is summarized in Table 3. As we can see in Table 3, the thermal adjustment time becomes smaller for larger b and smaller n. For each set of runs, τ_{adjust} grows approximately linearly with , although for , the dependency is more shallow. For larger values of , τ_{adjust} seems to have a stronger dependency on b. We speculate that the reason for increasing the value of τ_{adjust} for higher values of is that by increasing the opacity the energy transport via radiation becomes less efficient as the mean free path of the photon decreases. But it seems that there exists a threshold of efficiency, leading to a comparatively long adjustment time for the lowest values of .
We plot the vertical dependence of the mean free path of the photons ℓ = 1 /κρ normalized by the size of the domain L for Set C.
Fig. 7 Temperature T at z_{top} = 4 Mm versus time t for Run C5. The two horizontal lines mark the 1% margin around the final value of T and the vertical line marks the time τ_{adjust} ≈ 20 ks after which T lies within these margins. 
Fig. 8 Normalized mean free path of photons ℓ/L versus height for Set C. The dots represent the surface τ ≈ 1. 
As we can see in Fig. 8, the mean free path increases by several orders of magnitude from the bottom of the domain to the top. Furthermore, ℓ is larger for smaller . In the optically thick part, the difference in ℓ is one order of magnitude, which is equal to a corresponding change in . In the optically thin part, the difference in the values of ℓ becomes smaller, as we reach the top of the domain. For , ℓ is the smallest and, at the bottom of the domain, three orders of magnitude smaller than for . Furthermore, ℓ is 10 times the size of the domain for , which makes the cooling more efficient. We would have expected to see a large change in the mean free path as we go through the surface. Nevertheless, the exponential growth seems to be roughly the same throughout the domain, at least for the smallest value of .
3.8. Properties of an atmosphere with undefined n
By choosing a = −1 and b = 3, we have a constant heat conductivity K that is independent of density and temperature as the heat conductivity is given by Eq. (7). The nominal value of n is given by (20)In this case, since K = const, we expect to have only a polytropic solution which satisfies the thermostatic equilibrium if ∇_{z}T = const, but it is then unclear how ρ varies. In Fig. 9, we plot the profiles of density, temperature and entropy for all three runs of Set E.
Fig. 9 Density, temperature, and entropy profiles for Set E (n = 0/0). Dashed, dotted, and solid lines represent E4, E5, and E6, respectively. The red dots present the surface of the model where τ = 1. 
The first panel shows that in the optically thick part the density is nearly the same for the three values of , while in the optically thin part it decreases with increasing . For higher values of , the density drops faster than in the case of smaller . In all cases, K is constant in both the optically thick and thin parts, but an interesting aspect is that its bottom value K_{bot} is of the same order of magnitude as in Sets A, B, and C. In the second panel of Fig. 9, we plot the temperature profiles. As expected, there is no isothermal part. The slope of temperature decreases approximately linearly as we go to higher values of , because K is proportional to . Although the solutions show no transition from the polytropic part to an isothermal one, the atmosphere has still a layer where τ = 1, which is shown as red dots in all panels of Fig. 9. In contrast to the other sets, A, B, C, and D, the temperature profiles look qualitatively different. As in Sets A, B, and C, in the optically thick part, the different temperature profiles have nearly the same gradient, while in Set E, the gradient is different for the three values of . This is because in this case, thermostatic equilibrium is obeyed with a constant value of K (independent of z).
In the third panel of Fig. 9, we plot entropy profiles for the three values of . In all cases the entropy increases with a slope that depends on .
Fig. 10 b versus a for different values of the polytropic index n. The black dots represents the combination of a = 1 with different values of b which are used in the main simulation; see Table 3. The stars represent the combination of a = 2 and squares represent a = 0.5 with different values of b; see Table 4. 
The actual polytropic index can be computed from the resulting superadiabatic (or entropy) gradient, (21)where ∇_{ad} = 1 − 1 /γ is the adiabatic gradient. This gives ∇, which is related to the actual n via (22)which follows from the perfect gas relation p ∝ ρT. We note also that convection is not possible in one dimension, so we obtain directly the hydrostatic solution, which may be unstable.
By solving Eq. (8) for b, we obtain b = 3 − n(1 + a); see Fig. 10. We see that for different values of n, the graphs of b versus a intersect each other at one common point. This corresponds to K = K_{0} = const.; see Eq. (7). This means that the solution for constant K can belong to any of these polytropic indexes. For a = −1 and b = 3, in which case n is undefined, the solutions have a value of n_{actual} that depends on the height of the domain. Using Eq. (22) together with hydrostatic equilibrium, dp/ dz = −ρg, we find (23)so n_{actual} increases as z_{top} is increased. However, this increase is partially being compensated by a small simultaneous decrease of T_{top}, which reduces the increase of n_{actual} by about 10% as z_{top} increases. When the domain is sufficiently thin, the value of n_{actual} drops below the critical value (γ − 1)^{1} = 3/2, so the system would be unstable to the onset of convection. We return to the relation between n and the height of the domain in Sect. 3.12, where we consider solutions using the optically thick approximation with a radiative boundary condition at the top.
3.9. Dependence on the size of the domain
In our model, the size of the domain plays an important role in getting the polytropic and isothermal solutions for the temperature profile. The domain has to be big enough so that the transition point lies inside the domain. In Fig. 11, we show the vertical dependence of temperature for six domain sizes for Run A5. If the size of the domain is z< 7 Mm, it is too small to yield the isothermal part where ∇_{z}T = 0 and a boundary layer is produced. The opacity is then too large to let the heat be radiated away. A size of around z = 8 Mm is sufficient to get the isothermal part. However, a domain size that is too large (z = 10 Mm) leads to numerical difficulties near the top boundary, especially if the resolution is too low. For all the runs shown in Table 3, we have always started by performing several test simulations to find a suitable domain size.
Fig. 11 Temperature gradient (upper panel) and temperature profile (lower panel) of seven different sizes of the domain z = 3, 4, 5, 6, 7, 8, and 10 Mm of Run A5. We note that the lines for z ≥ 6 all fall on top of each other. 
3.10. Radiative diffusivity
In numerical simulations, the radiative diffusivity χ is an important parameter, and has the same dimension as the kinematic viscosity ν. Both χ and ν determine whether the results of numerical turbulence simulations are reliable or not and whether they are able to dissipate all the energy within the mesh. In a numerical simulation we are restricted to a certain number of grid points. If the diffusion of the temperature in a simulation is very small, it can happen that the changes in the temperature are too large over the distance of neighboring grid points. Hence, the changes in the temperature cannot be resolved in such a simulation. Therefore, it is important to measure how large are the thermal diffusivity in our models of a radiative atmosphere. The Péclet number is a dimensionless number that quantifies the importance of advective and diffusive term, which is here defined as (24)where H_{p} is a pressure scale height and u_{rms} is rms velocity. The radiative diffusivity is defined as (25)where K is evaluated using Eq. (7). As we do not solve for a velocity equation in our model, we use instead the sound speed, which can be related to u_{rms} via the Mach number Ma = u_{rms}/c_{s}. The normalized Péclet number in our simulation is then given by (26)As an example, we plot for Set C in Fig. 12. As we can see in Fig. 12, is a large number for the optically thick part and it decreases as we go toward the optically thin part. This can be explained with Eq. (25), where χ is proportional to K. In the optically thin part, K increases, so χ also increases. As a result, decreases. Furthermore, is larger for the larger value of .
Fig. 12 versus z for Set C using different values of : (dotteddashed line), (dashed line), (dotted line) and (solid line). 
The quantity is a measure of the ratio of KelvinHelmholtz time to the sound travel time, τ_{sound} = d/c_{s}. In our case, τ_{sound} ≈ 0.1 ks. Looking at Fig. 12, one sees that is proportional to and thus proportional to the thermal adjustment time. depends on z, but in the middle of the layer at z = 2 Mm we have . This time is much longer than the response time to general threedimensional disturbances (Spiegel 1987).
3.11. The same polytropic index with different a and b
As we can see in Fig. 10, for a certain value of the polytropic index, we can choose different combinations of a and b. For each value of n that we have in Table 3, we choose two different other combinations of a and b with the same value of . For example for the polytropic index n = 1 we choose two other combinations as a = 0.5 and b = 1.5 for one set and a = 2 and b = 0 for another one (see Table 4). We run eight more simulations with the same initial conditions as in previous runs and we obtain a similar equilibrium solution for the same polytropic index n. We calculate the effective temperature and the position where τ ≈ 1 as reference parameters with our old runs. The results are summarized in Table 4. For each set of runs with the same polytropic index, we labeled the runs similarly to those in Table 3.
Summary of the results for different values of a and b with the same polytropic index n.
As we see in Table 4, for each set of runs the effective temperature does not vary strongly, but there is a systematic behavior. By increasing the value of a, the effective temperature increases when n< 3/2 and decreases when n ≥ 3/2, but the surface is shifted to the lower part of the domain for all sets. The stratification of temperature and other important properties of these atmospheres can be explained analogously to those of Sets A, B, C and D. As an example, we plot in Fig. 13 the temperature profiles (upper panel) and radiative diffusivity χ (lower panel) for Set H. In the optically thick part, χ is the same for different combinations of a and b for the same n. However, in the optically thin part, χ becomes larger for larger values of a. This can be explained using Eqs. (7) and (25) to show that in the upper isothermal part, χ increases with decreasing ρ like ρ^{− (a + 2)}. Thus, we can conclude that, even though χ is the same and the solution similar to the optically thick part, there are differences in the optically thin part.
Fig. 13 Profile of T and χ for Set H. In the upper panel, the red dots denote the locations of τ = 1. 
3.12. Optically thick case with radiative boundary
To compare our results with those in the optically thick approximation, we adopt the radiative boundary condition, (27)and keep all other conditions the same as in the radiative transfer calculation, except that in Eq. (3) is replaced by K∇^{2}T. Here, we have assumed K to be constant, so we shall from now on refer to its value as K_{0}, so our solutions will be polytropes with constant polytropic index n = dlnρ/ dlnT and constant doublelogarithmic temperature gradient ∇ = dlnT/ dlnp.
The value of ∇ = 1/(1 + n) can be computed from the equations governing hydrothermal equilibrium, (28)which yields (29)Such a model is characterized by choosing values for n and K_{0}. This is analogous to the case with radiative transfer, where n and are specified, and is related to K_{0} via Eq. (18). Here, it is convenient to define a nondimensional radiative conductivity as (30)The radiative flux is then given by (31)so we get the temperature at the top immediately as (32)Since the temperatures at top and bottom are now known, the thickness of the layer cannot be chosen independently and is instead given by (33)which is equivalent to Eq. (23) used before. Again, the fact that the value of d cannot be chosen independently is analogous to the case with radiative transfer, where the thickness of the optically thick layer with nearly constant K emerges as a result of the calculation. In Table 5 we present models for the same parameters as in Table 3, where T_{eff} = T_{top} in the optically thick model. In agreement with our radiative transfer calculations, we have here treated (instead of ) as our main input parameter (in addition to n). We have used Eq. (18) to convert into K_{0} and then used Eq. (30) to compute . It turns out that there is good agreement regarding the values of d, ρ_{top}, and T_{top} between the optically thick approximation using a radiative upper boundary condition and the radiative transfer calculations. However, unlike Fig. 6, the data in Table 5 show power law dependence of T_{eff} versus with the same exponent of 1/4 in all cases. To characterize the strength of density and temperature stratification, we also list the ratios ρ_{bot}/ρ_{top} and the ratio of the pressure scale height at the top to the thickness of the layer, . As expected, smaller values of ξ are reached by increasing the value of , but even for the smallest values of ξ are 0.03 for n = 3.25 and 0.08 for n = 1.
Summary of model parameters as a function of n and as obtained from the optically thick approximation with radiative upper boundary condition.
Fig. 14 Comparison of velocity and entropy distribution in twodimensional convection using the optically thick approximation with a radiative upper boundary condition (upper panel) and radiative transfer (lower panel). In both cases we have Pr = 100 and , corresponding to Ra = 3.6 × 10^{4}. In the lower panel, the dashed line gives the contour τ = 1. In order to compare similar structures in the two plots, we have extended the color table of the lower panel to slightly more negative values of s and have clipped it at high values, which are dominated by the strong increase of s above the τ = 1 surface. 
We emphasize that the only place where the choice of density enters our calculation is in Eq. (18) when we convert into K_{0}. As already indicated at the end of Sect. 2.5, an increase of ρ_{0} by some factor is equivalent to an increase of by the same factor. We note here that ρ_{0} enters both as the initial density at the bottom and in the definition of opacity through Eq. (16). The latter ensures that the opacity only changes through changes in , and not also through changes in ρ_{0}.
3.13. Convection
We now consider twodimensional convection and compare again results from the optically thick approximation using a radiative upper boundary condition with a calculation using radiative transfer. The control parameter characterizing onset and amplitude of convection is the Rayleigh number, (34)where is the superadiabatic gradient of the unstable, nonconvecting hydrostatic reference solution, and is the pressure scale height in the middle of the optically thick layer. Furthermore, we define the Prandtl number as Pr = ν/χ_{mid}, where ν is assumed to be constant and χ_{mid} is the radiative diffusivity in the middle of the optically thick layer; see Appendix B for details and an example. In Table 6 we list the values of the product Pr Ra, as well as d and χ_{mid} for models with n = 1 and different values of . We adopt periodic boundary condition in the x direction over a domain with side length L_{x}. When we adopt the diffusion approximation we take z_{top} = d, where d is calculated from Eq. (33) and given in Table 6. The midlayer is then at z = d/ 2, which is also the case when using radiative transfer, where the value of z_{top} (>d) was chosen to be sufficiently large and a = b = 1 is chosen to yield n = 1 (see Table 3).
Summary of model parameters as a function of for n = 1.
We determine the critical value for the onset of convection by calculating the rms velocity in the domain, u_{rms}, for different values of and extrapolate to u_{rms} → 0. For the final to initial bottom density ratio, ρ_{bot}/ρ_{0}, we use Eq. (C.5) derived in Appendix C. Since is proportional to χ and χ_{mid}, we can compute the product Pr Ra using Eq. (34). It turns out that for Pr = 1, the critical value is at (corresponding to Ra ≲ 710) in the optically thick approximation. This corresponds to , which is too small to obtain a proper polytropic lower part. Therefore we choose in the following Pr = 100, in which cases the critical value for the onset of convection is at (corresponding to Ra ≲ 6600).
In the following, we take , so Ra = 3.6 × 10^{4}, d = 2.70 Mm, ν = 1.8 Mm km s^{1} (corresponding to ν = 1.8 × 10^{13} cm^{2} s^{1}), and ; see the sixth row of Table 6. We choose L_{x} = 14 Mm, which is large enough to accommodate two convection cells into the domain; see Fig. 14. The τ = 1 surface in the radiative transfer calculation agrees approximately with the height expected from the optically thick models using a radiative upper boundary condition. The flow is only weakly supercritical and therefore not very vigorous, which is also reflected by the fact that the τ = 1 surface is nearly flat. In the radiative transfer calculation, the specific entropy increases sharply with height above the τ = 1 surface. We note also that the characteristic narrow downdrafts of the optically thick calculation are now much broader when radiation transfer is used. Furthermore, the expected entropy minimum near the surface is virtually absent in the latter case; see also the middle panel of Fig. 15. This is because near z = d, the local value of χ is rather large (see the lower panel of Fig. 13), so the thickness of the thermal boundary layer becomes comparable to d itself.
Fig. 15 Comparison of vertical temperature, entropy, and velocity profiles at different x positions for the model shown in Fig. 14. The black dotted lines refer to the run with diffusion approximation and radiative upper boundary condition and the red solid lines to the run with radiative transfer. 
Fig. 16 Similar to the lower panel of Fig. 14, but at a later time (t = 500 ks), when the solution has switched into a singlecell configuration. 
In the model using the diffusion approximation, the temperature variations are much larger than in the model with radiative transfer; see Fig. 15. In the latter case, the velocities overshoot into the upper stably stratified layer, and the downflows are much broader and therefore slower than in the diffusion approximation.
It turns out that the solution with radiative transfer shown in Fig. 14 is quasistable until about 350 ks (≈ 4 days) and then switches into a singlecell configuration with a fairly isolated updraft; see Fig. 16. We have seen similar behavior in other cases with radiation too, and it is possible that this is a consequence of our setup. Firstly, the restriction to twodimensional convection is a serious artifact. Secondly, the assumption of a fixed temperature at the bottom was a mathematical convenience, but it is not physical motivated.
4. Increasing the density contrast
As alluded to in the introduction, the inclusion of the physics of radiation implies the occurrence of σ_{SB} as an additional physical constant that couples the resulting temperature and density contrasts to changes in the Rayleigh number. Given that we have already made other simplifications such as the negligence of hydrogen ionization, we end up with rather small density contrasts of less than ten when n = 1, as seen in Table 5. By considering σ_{SB} an adjustable parameter, we can alleviate this constraint. We demonstrate this in Table 7, where we increase the value of σ_{SB} from its physical value of 5.67 × 10^{20} g cm^{3} km^{3} s^{3} K^{4} by eight orders of magnitude, keeping however K_{0} fixed. We have chosen here K_{0} = 1.3 × 10^{7} = g cm^{3} km^{3} s^{3} Mm K^{1}, which corresponds to the model in the sixth row of Table 6. Since σ_{SB} enters the definition of , its value is now no longer the same, even though K_{0} is. The Rayleigh numbers change slightly, because they depend on the values of density and temperature in the middle of the domain, which do of course change.
Density contrast and other model parameters as a function of σ_{SB} for n = 1 and K_{0} = 1.3 × 10^{7}.
Large density contrasts are one of the important ingredients in modeling the physics of sunspot formation by surface effects such as the negative effective magnetic pressure instability; see Brandenburg et al. (2013) for a recent model. Including radiation into such still rather idealized models and to study the relation to other competing or corroborating mechanisms such as the one of Kitchatinov & Mazur (2000) was indeed an important motivation behind the work of the present paper.
5. Conclusions
The inclusion of radiative transfer in a hydrodynamic code provides a natural and physically motivated way of placing an upper stably stratified layer on top of an optically thick layer that may be stably or unstably stratified, which of the two depends on the opacity. Using a Kramerslike opacity law with freely adjustable exponents on density and temperature yields polytropic solutions for certain combinations of the exponents a and b. The prefactor in the opacity law determines essentially the values of the Péclet and Rayleigh numbers. However, in contrast to earlier studies of convection in polytropic layers, the temperature contrast is no longer a free parameter and increases with increasing Rayleigh number – unless one considers the StefanBoltzmann “constant” as an adjustable parameter. The physical values of the prefactor on the opacity are much larger than those used here, but larger prefactors lead to values of the radiative diffusivity that become eventually so small that temperature fluctuations on the mesh scale cannot be dissipated by radiative diffusion. In previous work (Nordlund 1982; Steffen et al. 1989; Vögler et al. 2005; Heinemann et al. 2007; Freytag et al. 2012), this problem has been avoided by applying numerical diffusion or using numerical schemes that dissipate the energy when and where needed. However, this may also suppress the possibility of physical instabilities that we are ultimately interested in. This motivates the investigation of models with prefactors in the Kramers opacity law that are manageable without the use of numerical procedures to dissipate energy artificially.
It turns out that in all cases with a and b such that n> − 1, the stratification corresponds to a polytrope with index n below the photosphere and to an isothermal one above it. This was actually expected given that such a solution has previously been obtained analytically in the special case of constant κ (corresponding to a = b = 0); see Spiegel (2006). On the other hand, the isothermal part was apparently not present in the simulations of Edwards (1990).
Contrary to the usual polytropic setups (e.g., Brandenburg et al. 1996), the temperature contrast is now no longer an independent parameter but it is tied essentially to the Rayleigh number. Significant temperature contrasts can only be achieved at large Rayleigh numbers, which corresponds to a large prefactor in front of the opacity. While this aspect can be reproduced already with polytropic models using the radiative boundary conditions, there are some surprising differences between the two. Most important is perhaps the fact that the specific entropy must increase above the unstable layer, while with a radiative boundary condition the specific entropy always decreases. Although the differences in the resulting temperature profiles are small, there are major differences in the flows speeds in the two cases. We also find that at late times the convection cells in simulations with full radiative transport tend to merge into larger ones. Whether or not this is an artefact of our restriction to twodimensional flows remains open. In this connection, we should also point out the presence of a geometric correction factor in front of the radiative heating and cooling term in Eq. (15) that is needed to reproduce the correct cooling rate, but it does not affect the steady state solution.
Comparing with realistic simulations of the Sun, there is not really an isothermal part, but a pronounced sudden drop in temperature followed by a continued decrease in temperature (see, e.g., Stein & Nordlund 1998). On the other hand, in our simulations there is no jump in the temperature profile near the surface and the atmosphere changes smoothly from polytropic to isothermal. We suspect that the reason for this difference is that in our models ionization effects are ignored, while in the solar atmosphere the degree of ionization of hydrogen increases with depth. In the Sun, the density decreases significantly from the upper part of the convection zone as we go to the photosphere. This makes the opacity smaller and the atmosphere in the photosphere becomes transparent. At the height where the ionization temperature of hydrogen is reached, the H^{−} opacity becomes important, which is not included in our simulations. The radiative heat conductivity in our simulations is found to be constant throughout the optically thick part and then increases sharply in the optically thin part. Solving this in the optically thick approximation, which has sometimes been done, becomes computationally expensive and even unphysical, so radiative transfer becomes a viable alternative for studying layers that are otherwise polytropic in the lower part of the domain.
Online material
Appendix A: Cooling rate and correction factor
The purpose of this appendix is to show that for onedimensional temperature perturbations, the correct cooling rates are obtained with just two rays if the D/ 3 correction factor in Eq. (15) is applied. Similar considerations apply also to the case of twodimensional problems. Cooling rates are important for understanding temporal aspects such as the approach to the final state (Sect. 3.1) or the thermal adjustment time (Sect. 3.7). Thus, the equilibrium solutions discussed in the other sections are not affected by the following considerations.
The source of the problem lies in the fact that the 4π angular integration in Eq. (4) becomes inaccurate in one dimension and dependent on the optical thickness. In the optically thick regime, the diffusion approximation holds, so the cooling rate is proportional to K, which has a 1/3 factor in Eq. (7). In one dimension, one uses only the two rays in the vertical direction, so one misses the 1/3 factor and has to apply it afterwards to account for the “redundant” rays in the other two coordinate directions that show no variation. This is what is done in Eq. (15). However, in the optically thin limit, the mean free path becomes infinite and cooling is now possible in all three directions. In that case, the onedimensional approximation is not useful. To explain this in more detail, we begin by considering first the general case of threedimensional perturbations with wavevector k (Spiegel 1957). In that case, one can use the Eddington approximation to solve the transfer equation for the mean intensity, J = ∫I dΩ/4π, (A.1)so the cooling rate (for threedimensional perturbations) is (Unno & Spiegel 1966; Edwards 1990) (A.2)It is convenient to introduce here a photon diffusion speed as (A.3)and to write Eq. (A.2) in the form (A.4)where c_{γ}ℓ/ 3 = χ is the radiative diffusivity, as defined in Eq. (25), and ℓ = 1 /κρ is the local meanfree path of photons.
Solving Eq. (5) for two rays corresponds to solving Eq. (A.1) without the 1/3 factor. We would then obtain Eq. (A.4) without the two 1/3 factors. This would evidently violate the wellknown cooling rate χk^{2} in the optically thick limit, but in the optically thin limit it would be in agreement with Eq. (A.4), because the two 1/3 factors would cancel for large values of ℓ. However, we have to remember that temperature perturbations are here assumed onedimensional, so the intensity can only vary in the z direction, while the rays still go in all three directions. This means that under the sum in Eq. (15) only one third of the I − S terms give a contribution, and that the cooling rate is therefore (A.5)which has now only a single 1/3 factor. Likewise, if we had twodimensional perturbations such as in twodimensional convection considered in Sect. 3.13, only 2/3 of the terms under the sum in Eq. (15) would contribute. However, in a twodimensional radiative transfer calculation, the additional 1/3 would be absent, which explains the D/ 3 correction factor with D = 2 in this case.
Fig. A.1 Dependence of the cooling rates computed from models with different values of (from 10^{2} to 10^{5} Mm^{1} cm^{3} g^{1}) and k_{z} (=1, 2, and 4 indicated by diamonds, triangles, and squares, respectively). 2D models with four and eight rays are indicated by crosses and circles, respectively, while 3D models with six rays are shown as plus signs. The red solid line corresponds to Eq. (A.5), the dashed blue line to Eq. (A.4), and the dotted line with open circles to the case without correction factor. 
We have verified that with the correction factor in place, the code now yields the same cooling rates in both the optically thick and thin regimes, regardless of the numbers of rays used. This is shown in Fig. A.1, where we plot cooling rates for different values of in a domain of size 2π (in Mm), so the smallest wavenumber is 1 Mm^{1}. With ρ = 4 × 10^{4} g cm^{3} the photon meanfree path varies from 0.025 to 25 Mm as is decreased from 10^{5} to 10^{2} Mm^{1} cm^{3} g^{1}. For the Kramers opacity, we use the exponents a = 1 and b = 0. (No gravity is included here, so there would be no convection.) The temperature is 38,968 K, as before, which yields c_{γ} = 3.87 km s^{1} for the photon diffusion speed. There is excellent agreement between 1D cases with correction factor and the 3D calculation (with onedimensional perturbation). However, the 2/3 correction factor in the 2D calculation (both with four and with eight rays) seems to be systematically off and should instead by around 0.8 for better agreement. However, as discussed before, the correction factor does not affect the steady state and therefore also not the results presented in Sect. 3.13. The diffusion approximation would imply λ = (c_{γ}k/ 3)ℓk = χk^{2}, which corresponds to the diagonal in Fig. A.1 and agrees with the red solid line for ℓk ≲ 0.5.
For threedimensional perturbations, the correct cooling rate in the optically thin regime is three times faster than for onedimensional perturbations. This is because now the radiation goes in all three directions. Solutions to threedimensional perturbations clearly cannot be reproduced in less than three dimensions. However, for onedimensional perturbations, the correct cooling rate is now obtained with a onedimensional calculation both in the optically thin and thick regimes.
Appendix B: Expressions for Pr Ra and χ_{mid}
In Table 6 we listed the values of Pr Ra and χ_{mid} in the middle of the layer. The purpose of this appendix is to give the explicit expressions and to demonstrate the calculation with the help of an example. Since n = 1 was assumed, we have ∇ = (1 + n)^{1} = 1/2. Considering the case , Eqs. (31)−(33) yield F_{rad} = 0.00131 g cm^{3} km^{3} s^{3}, T_{top} = 12 320 K, and d = 2.70 Mm. Next, given that the temperature varies linearly, we compute the midlayer temperature as . This allows us to compute ρ_{mid} = ρ_{bot} (T_{mid}/T_{bot})^{n} = 2.2 × 10^{4} g cm^{3}, where ρ_{bot} = 3.3 × 10^{4} g cm^{3} is smaller than ρ_{0} by a factor ρ_{bot}/ρ_{0} = 0.83; see Appendix C. Thus, χ_{mid} = F_{rad}/ (ρ_{mid}g∇/∇_{ad}) = 0.0175 Mm km s^{1}, as well as . This yields , where ∇ − ∇_{ad} = 0.1.
Appendix C: Final to initial bottom density ratio
Initially, the stratification is isothermal, so the density is given by and the initial surface density is (C.1)In the final state, the stratification is polytropic, so the density is given by ρ(z) = ρ_{bot} [T(z) /T_{bot}] ^{n} and the surface density is (C.2)
Here, ρ_{bot} is the bottom density of the final state, which is different from the initial value ρ_{0}, as explained in Sect. 2.6. Integrating Eq. (C.2) and using dz/ dT = K_{0}/F_{rad} from Eq. (28) yields (C.3)Using Eq. (29) together with ∇ = 1/(1 + n) and , we have (C.4)Using mass conservation, we have Σ_{fin} = Σ_{ini}, so we obtain from Eqs. (C.1) and (C.4) (C.5)for the final to initial bottom density ratio.
Acknowledgments
We are indebted to Tobi Heinemann for his engagement in implementing the radiative transfer module into the Pencil code. We also thank the referee for many helpful comments that have led to improvements of the manuscript. This work was supported in part by the European Research Council under the AstroDyn Research Project No. 227952, and by the Swedish Research Council under the project grants 62120115076 and 20125797. We acknowledge the allocation of computing resources provided by the Swedish National Allocations Committee at the Center for Parallel Computers at the Royal Institute of Technology in Stockholm and the National Supercomputer Centers in Linköping, the High Performance Computing Center North in Umeå, and the Nordic High Performance Computing Center in Reykjavik.
References
 Brandenburg, A. 2005, ApJ, 625, 539 [Google Scholar]
 Brandenburg, A., Jennings, R. L., Nordlund, Å., et al. 1996, J. Fluid Mech., 306, 325 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Brandenburg, A., Kleeorin, N., & Rogachevskii, I. 2013, ApJ, 776, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Cattaneo, F., Brummell, N. H., Toomre, J., Malagoli, A., & Hurlburt, N. E. 1991, ApJ, 370, 282 [NASA ADS] [CrossRef] [Google Scholar]
 Edwards, J. M. 1990, MNRAS, 242, 224 [NASA ADS] [CrossRef] [Google Scholar]
 Gastine, T., & Dintrans, B. 2008, A&A, 484, 29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Freytag, B., Steffen, M., Ludwig, H.G., et al. 2012, J. Comput. Phys., 231, 919 [Google Scholar]
 Heinemann, T., Dobler, W., Nordlund, Å., & Brandenburg, A. 2006, A&A, 448, 731 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heinemann, T., Nordlund, Å., Scharmer, G. B., & Spruit, H. C. 2007, ApJ, 669, 1390 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2012, ApJ, 755, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J., Mantere, M. J., Cole, E., Warnecke, J., & Brandenburg, A. 2013, ApJ, 778, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Kippenhahn, R., & Weigert, A. 1990, Stellar structure and evolution (Berlin: Springer) [Google Scholar]
 Kitchatinov, L. L., & Mazur, M. V. 2000, Solar Phys., 191, 325 [Google Scholar]
 Miesch, M. S., Elliott, J. R., Toomre, J., et al. 2000, ApJ, 532, 593 [NASA ADS] [CrossRef] [Google Scholar]
 Mitra, D., Brandenburg, A., Kleeorin, N., & Rogachevskii, I. 2014, MNRAS, 445, 761 [NASA ADS] [CrossRef] [Google Scholar]
 Nordlund, Å. 1982, A&A, 107, 1 [Google Scholar]
 Nordlund, Å. 1985, Solar Phys., 100, 209 [Google Scholar]
 Spiegel, E. A. 1957, ApJ, 126, 202 [NASA ADS] [CrossRef] [Google Scholar]
 Spiegel, E. A. 1987, in The internal solar angular velocity, eds. B. R. Durney, & S. Sofia (Dordrecht: Reidel), 321 [Google Scholar]
 Spiegel, E. A. 2006, EAS Pub. Ser., 21, 127 [CrossRef] [EDP Sciences] [Google Scholar]
 Steffen, M., Ludwig, H.G., & Krüß, A. 1989, A&A, 123, 371 [NASA ADS] [Google Scholar]
 Stein, R. F., & Nordlund, Å. 1989, ApJ, 342, L95 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Stein, R. F., & Nordlund, Å. 1998, ApJ, 499, 914 [Google Scholar]
 Stein, R. F., & Nordlund, Å. 2012, ApJ, 753, L13 [Google Scholar]
 Warnecke, J., Losada, I. R., Brandenburg, A., Kleeorin, N., & Rogachevskii, I. 2013, ApJ, 777, L37 [NASA ADS] [CrossRef] [Google Scholar]
 Unno, W., & Spiegel, E. A. 1966, Publ. Astron. Soc. Japan, 18, 85 [Google Scholar]
 Vögler, A., Shelyag, S., Schüssler, M., Cattaneo, F., & Emonet, T. 2005, A&A, 429, 335 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stix, M. 2002, The Sun: an introduction (Berlin: SpringerVerlag) [Google Scholar]
All Tables
Summary of the results for different values of a and b with the same polytropic index n.
Summary of model parameters as a function of n and as obtained from the optically thick approximation with radiative upper boundary condition.
Density contrast and other model parameters as a function of σ_{SB} for n = 1 and K_{0} = 1.3 × 10^{7}.
All Figures
Fig. 1 Vertical temperature profile at five different times t = 0, 3, 30, 120, and 1578 ks for Run A6 with . Squares, circles and crosses represent different optical depths τ = 0.1, τ = 1 and τ = 10, respectively. The arrow represents the direction of the time evolution of the temperature profile. 

In the text 
Fig. 2 Temperature over height for Run A6 at different times t = 0, 0.01, 0.1, and 0.2 ks, plotted as solid, dotted, dashed and dasheddotted lines, respectively. The red dots correspond to the location of τ ≈ 1. 

In the text 
Fig. 3 Density, temperature and entropy of the equilibrium state versus height, from left to right, for Sets A, B, C and D, from top to bottom. The four different lines in each plot corresponds to the value of the rescaled opacity , 10^{5}, 10^{6}, 10^{7} Mm^{1} cm^{3} g^{1}. The dots in all plots represent the surface τ ≈ 1. The red dashed lines represent the initial profile of each set. 

In the text 
Fig. 4 Source function and vertical profiles of incoming and outgoing intensity near the surface for Run A7. The dashed line represents the source function and the solid lines represent the incoming intensity I^{+} (red) and outgoing intensity I^{−} (blue). The vertical lines represent the (constant) difference between I^{+} and I^{−}. 

In the text 
Fig. 5 Radiative heat conductivity K versus height for Set C. K is plotted for different values of where is shown by dotteddashed line, dashed line, dotted line and solid line. 

In the text 
Fig. 6 Effective temperature T_{eff} versus rescaled opacity for Sets A, C, and D. The crosses, circles and stars show the values of T_{eff} for different values of for Sets A, C, and D, respectively. Different lines correspond to line fit of T_{eff} with normalized opacity . 

In the text 
Fig. 7 Temperature T at z_{top} = 4 Mm versus time t for Run C5. The two horizontal lines mark the 1% margin around the final value of T and the vertical line marks the time τ_{adjust} ≈ 20 ks after which T lies within these margins. 

In the text 
Fig. 8 Normalized mean free path of photons ℓ/L versus height for Set C. The dots represent the surface τ ≈ 1. 

In the text 
Fig. 9 Density, temperature, and entropy profiles for Set E (n = 0/0). Dashed, dotted, and solid lines represent E4, E5, and E6, respectively. The red dots present the surface of the model where τ = 1. 

In the text 
Fig. 10 b versus a for different values of the polytropic index n. The black dots represents the combination of a = 1 with different values of b which are used in the main simulation; see Table 3. The stars represent the combination of a = 2 and squares represent a = 0.5 with different values of b; see Table 4. 

In the text 
Fig. 11 Temperature gradient (upper panel) and temperature profile (lower panel) of seven different sizes of the domain z = 3, 4, 5, 6, 7, 8, and 10 Mm of Run A5. We note that the lines for z ≥ 6 all fall on top of each other. 

In the text 
Fig. 12 versus z for Set C using different values of : (dotteddashed line), (dashed line), (dotted line) and (solid line). 

In the text 
Fig. 13 Profile of T and χ for Set H. In the upper panel, the red dots denote the locations of τ = 1. 

In the text 
Fig. 14 Comparison of velocity and entropy distribution in twodimensional convection using the optically thick approximation with a radiative upper boundary condition (upper panel) and radiative transfer (lower panel). In both cases we have Pr = 100 and , corresponding to Ra = 3.6 × 10^{4}. In the lower panel, the dashed line gives the contour τ = 1. In order to compare similar structures in the two plots, we have extended the color table of the lower panel to slightly more negative values of s and have clipped it at high values, which are dominated by the strong increase of s above the τ = 1 surface. 

In the text 
Fig. 15 Comparison of vertical temperature, entropy, and velocity profiles at different x positions for the model shown in Fig. 14. The black dotted lines refer to the run with diffusion approximation and radiative upper boundary condition and the red solid lines to the run with radiative transfer. 

In the text 
Fig. 16 Similar to the lower panel of Fig. 14, but at a later time (t = 500 ks), when the solution has switched into a singlecell configuration. 

In the text 
Fig. A.1 Dependence of the cooling rates computed from models with different values of (from 10^{2} to 10^{5} Mm^{1} cm^{3} g^{1}) and k_{z} (=1, 2, and 4 indicated by diamonds, triangles, and squares, respectively). 2D models with four and eight rays are indicated by crosses and circles, respectively, while 3D models with six rays are shown as plus signs. The red solid line corresponds to Eq. (A.5), the dashed blue line to Eq. (A.4), and the dotted line with open circles to the case without correction factor. 

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