© 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 compute-intensive 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 granulation-like 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 Stefan-Boltzmann 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 Kramers-like 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 two-dimensional 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 one-dimensional 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 two-dimensional 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 rate-of-strain 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 rate-of-strain tensor S) is given by the frequency-integrated Planck function, so S = (σ_SB/π)T⁴, where σ_SB is the Stefan-Boltzmann 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 frequency-integrated 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 Kramers-like 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 plane-parallel 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ⁿ, 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 Kramers-like 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 Schwarzschild-stable 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 stress-free, 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 volume-averaged 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¹, which solves the hydrodynamic differential equations with a high-order finite-difference 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 non-local with respect to each processor. There are two local parts that are compute-intensive and one that is non-local 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, non-local contribution and an intrinsic, local one. An analogous calculation is done for calculating τ along the geometric coordinate as , where l₀ 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₀ e^{τ₀ − τ}, where I₀ and τ₀ 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₀ = I(τ₀) 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 round-off 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.

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² 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⁷ 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 ρ = ρ₀exp( − z/H_p), where H_p = ℛT/μg is the pressure scale height, and ρ₀ is a constant that we set to ρ₀ = 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.

It is convenient to express κ in the form (16)where is a rescaled opacity and is related to κ₀ by ; where T₀ = T_bot is used. (By contrast, ρ₀ 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³ g^-1 (=10^-8 cm² g^-1). For each value of n, we choose 4 different values of , 10⁵, 10⁶, and 10⁷ Mm^-1 cm³ g^-1. We note that the actual Kramers opacity for free–free and bound–free transitions with a = 1 and b = −7/2 has κ₀ between 6.6 × 10²² and 4.5 × 10²⁴ cm⁵ g^-2K^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 one-dimensional 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⁴, 10⁵, 10⁶ and 10⁷ Mm^-1 cm³ 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.

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.

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.

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⁴K, whereas in Run D4 the value is T ≈ 2.9 × 10⁴K. A similar behavior can also be seen for the other values of .

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 (∇_zs> 0) for Set A, while for Set B it is constant and equal to zero (∇_zs ≈ 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, ∇_zs< 0. This means that their atmospheres are convectively unstable. (Convection will of course not occur in our one-dimensional model, but we will obtain the so-called 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⁴ 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.

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.

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.

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 ∇_zT = 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.

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 .

The actual polytropic index can be computed from the resulting super-adiabatic (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₀ = 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 ∇_zT = 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.

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 .

The quantity is a measure of the ratio of Kelvin-Helmholtz 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 three-dimensional 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.

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.

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∇²T. Here, we have assumed K to be constant, so we shall from now on refer to its value as K₀, so our solutions will be polytropes with constant polytropic index n = dlnρ/ dlnT and constant double-logarithmic 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₀. This is analogous to the case with radiative transfer, where n and are specified, and is related to K₀ via Eq. (18). Here, it is convenient to define a non-dimensional 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₀ 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.

We emphasize that the only place where the choice of density enters our calculation is in Eq. (18) when we convert into K₀. As already indicated at the end of Sect. 2.5, an increase of ρ₀ by some factor is equivalent to an increase of by the same factor. We note here that ρ₀ 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 ρ₀.

3.13. Convection

We now consider two-dimensional 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, non-convecting 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 mid-layer 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).

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/ρ₀, 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⁴, d = 2.70 Mm, ν = 1.8 Mm km s^-1 (corresponding to ν = 1.8 × 10¹³ cm² 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.

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 quasi-stable until about 350 ks (≈ 4 days) and then switches into a single-cell 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 two-dimensional 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³ s^-3 K^-4 by eight orders of magnitude, keeping however K₀ fixed. We have chosen here K₀ = 1.3 × 10^-7 = g cm^-3 km³ 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₀ 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.

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 Kramers-like 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 Stefan-Boltzmann “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 two-dimensional 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.

---

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 621-2011-5076 and 2012-5797. 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

Online material

Appendix A: Cooling rate and correction factor

The purpose of this appendix is to show that for one-dimensional 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 two-dimensional 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 one-dimensional approximation is not useful. To explain this in more detail, we begin by considering first the general case of three-dimensional 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 three-dimensional 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 mean-free 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 well-known cooling rate χk² 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 one-dimensional, 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 two-dimensional perturbations such as in two-dimensional convection considered in Sect. 3.13, only 2/3 of the terms under the sum in Eq. (15) would contribute. However, in a two-dimensional radiative transfer calculation, the additional 1/3 would be absent, which explains the D/ 3 correction factor with D = 2 in this case.

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 mean-free path varies from 0.025 to 25 Mm as is decreased from 10⁵ to 10² Mm^-1 cm³ 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 one-dimensional 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², which corresponds to the diagonal in Fig. A.1 and agrees with the red solid line for ℓk ≲ 0.5.

For three-dimensional perturbations, the correct cooling rate in the optically thin regime is three times faster than for one-dimensional perturbations. This is because now the radiation goes in all three directions. Solutions to three-dimensional perturbations clearly cannot be reproduced in less than three dimensions. However, for one-dimensional perturbations, the correct cooling rate is now obtained with a one-dimensional 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³ s^-3, T_top = 12 320 K, and d = 2.70 Mm. Next, given that the temperature varies linearly, we compute the mid-layer temperature as . This allows us to compute ρ_mid = ρ_bot (T_mid/T_bot)ⁿ = 2.2 × 10^-4 g cm^-3, where ρ_bot = 3.3 × 10^-4 g cm^-3 is smaller than ρ₀ by a factor ρ_bot/ρ₀ = 0.83; see Appendix C. Thus, χ_mid = F_rad/ (ρ_midg∇/∇_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] ⁿ and the surface density is (C.2)

Here, ρ_bot is the bottom density of the final state, which is different from the initial value ρ₀, as explained in Sect. 2.6. Integrating Eq. (C.2) and using dz/ dT = K₀/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.

All Tables

All Figures