© ESO, 2011

1. Introduction

The internetwork solar chromosphere is continuously pervaded by acoustic waves and shocks with periods of around 3 min. A large literature exists on this subject. Some examples are: Krijger et al. (2001), who studied upper-photospheric and low-chromospheric oscillations observed in the UV continuum around 150 nm with the Transition Region And Coronal Explorer (TRACE). Rutten et al. (2008) studied the response of chromospheric diagnostics to photospheric events such as granular buffeting and exploding granules as observed with the Dutch Open Telescope (DOT). Somewhat longer ago, Lites et al. (1993) studied solar chromospheric oscillations using spectrograms in Ca II  H obtained with the Vacuum Tower Telescope at NSO/Sacramento Peak. Wikstøl et al. (2000) studied upper-chromospheric and transition region oscillations using spectral time series obtained around 103 nm with the SUMER spectrograph aboard SOHO. These and other studies confirm the picture of acoustic wave excitation by granular dynamics in the photosphere, the waves then propagate predominantly in the vertical direction and evolve into shocks. Depending on the magnetic field topology they occasionally retain their identity all the way up into the transition zone (McIntosh et al. 2001). The shock fronts are high temperature disturbances in an otherwise cool background atmosphere.

Observations of CO lines in quiet sun areas (e.g., Noyes & Hall 1972; Ayres & Testerman 1981; Ayres et al. 2006) indicate temperatures low enough to form significant amounts of CO molecules. Static 1D semi-empirical models designed to reproduce such observations can have temperatures down to 2750 K (Ayres & Rabin 1996). These low temperatures in such static models correspond to the cool inter-shock phases of the quiet chromosphere.

This scenario was confirmed in the one-dimensional (1D) simulations of Carlsson & Stein (1997) that explained the formation of Ca ii K_2V and H_2V bright grains as an effect of upward-propagating shock waves excited in the photosphere. The first attempt to model such wave-response in the chromosphere with multidimensional numerical simulations was made by Skartlien et al. (2000). They studied the response of the chromosphere to collapsing granules using a 3D radiation-hydrodynamic (RHD) simulation spanning from the upper convection zone to 1.2 Mm above ⟨τ₅₀₀⟩ = 1, and again confirmed the picture of the internetwork chromosphere as a cool background state pervaded by waves and shocks. Wedemeyer et al. (2004) performed a 3D RHD simulation with higher spatial resolution but less sophisticated radiative transfer. The higher resolution allowed them to distinguish different shock front shapes depending on the excitation mechanism: spherical fronts excited by pressure fluctuations in intergranular lanes; planar fronts, excited by simulation-box oscillations, the simulated counterpart of solar p-modes; and the previously described irregular front shapes caused by collapsing granules. The simulations employed simplified radiative transfer and an equation of state (EOS) based on Saha ionization equilibrium for hydrogen, an assumption that fails spectacularly in the chromosphere (Carlsson & Stein 2002). Thus, these 3D simulations properly model the shock compression and post-shock expansion, but inaccurately incorporate radiative heating, and give erroneous temperatures from the mass density and internal energy because of the assumption of Saha ionization equilibrium.

The 1D simulations performed by Carlsson & Stein (1997) modeled the radiative losses in great detail, including the effect of slow hydrogen ionization/recombination on the equation of state. However, the 1D geometry affected the temperatures in the chromosphere, leading to too high temperatures in the shocks as well as in the inter-shock phases. Leenaarts et al. (2007) performed a 2D radiation-magneto-hydrodynamic (RMHD) simulation with an EOS that took non-equilibrium hydrogen ionization into account using approximations for the radiation field based on Sollum (1999). This simulation, however, still computed the radiative losses based on Saha equilibrium, and required an ad-hoc heating term that prevented the temperature to drop below 2400 K.

All of the above simulations were limited in such a way that an accurate prediction of the minimum temperature occurring in the quiet solar chromosphere could not be made.

In this paper we discuss the temperature structure in a simulation that tries to remedy the various limitations of the previous models: it is two-dimensional, includes a non-equilibrium equation of state, includes radiative losses without the assumption of Saha equilibrium and includes heating through absorption of coronal radiation. We do not include a significant magnetic field, so that our simulation serves as a baseline against which simulations with stronger field can be compared.

In Sect. 2 we discuss the RMHD model assumptions in detail and discuss some of the basic physical processes that occur in the quiet chromosphere. In Sect. 3 we discuss the results of our run, especially paying attention to the accuracy of the radiative heating. We finish with a discussion and our conclusions in Sect. 4.

2. The model

We model the chromosphere using the RMHD code Bifrost (Gudiksen et al. 2011).

MHD equations.

The Bifrost code solves the equations of resistive MHD including the effects of heat conduction and radiation: $\begin{array}{ccc}& & \frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{u}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\mathit{\partial \rho }{u}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{(}\mathit{\rho }{u}{u}\mathrm{-}{\mathit{\tau }}_{\mathrm{=}}\mathrm{)}\mathrm{=}\mathrm{-}\mathrm{\nabla }\mathit{P}\mathrm{+}{j}\mathrm{\times }{B}\mathrm{+}\mathit{\rho }{g}\mathit{,}\\ & & \frac{\mathit{\partial e}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{e}{u}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathit{P}\mathrm{(}\mathrm{\nabla }\mathrm{·}{u}\mathrm{)}\mathrm{-}\mathrm{\nabla }\mathrm{·}{F}\mathrm{c}\mathrm{-}\mathrm{\nabla }\mathrm{·}{F}\mathrm{r}\mathrm{+}{\mathit{Q}}_{\mathrm{visc}}\mathrm{+}{\mathit{Q}}_{\mathrm{Joule}}\mathit{,} \\ & & {\mathit{\mu }}_{\mathrm{0}}{j}\mathrm{=}\mathrm{\nabla }\mathrm{\times }{B}\mathit{,}\\ & & \frac{\mathit{\partial }{B}}{\mathit{\partial t}}\mathrm{=}\mathrm{\nabla }\mathrm{\times }\mathrm{(}{u}\mathrm{\times }{B}\mathrm{)}\mathrm{-}\mathrm{\nabla }\mathrm{\times }{\mathit{\eta }}_{\mathrm{=}}{j}\mathit{,}\end{array}$with ρ the mass density, u the velocity, ${\mathit{\tau }}_{{\mathrm{=}}_{}}$ the viscous stress tensor, P the gas pressure, j the current density, B the magnetic field, g the gravitational acceleration, e the internal energy per volume, F_c the energy flux due to heat conduction, F_r the radiative energy flux, Q_visc the viscous energy dissipation, Q_Joule the Joule heating, μ₀ the permeability of free space and ${\mathit{\eta }}_{{\mathrm{=}}_{}}$ the electrical resistivity tensor. Bifrost employs a staggered grid and uses sixth-order operators to compute spatial derivatives. The time-stepping is done using a third-order predictor-corrector scheme by Hyman (1979), modified for variable timestep.

Equation of state.

Conservation of internal energy is expressed as $\begin{array}{ccc}\mathit{e}& \mathrm{=}& \frac{\mathrm{3}\mathit{kT}}{\mathrm{2}}\left({\mathit{n}}_{\mathrm{e}}\mathrm{+}{\mathit{n}}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{+}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathrm{l}}}{\mathit{n}}_{\mathit{i}}\mathrm{+}{\mathit{n}}_{\mathrm{o}}\right)\mathrm{+}{\mathit{n}}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{(}{\mathit{e}}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{+}{\mathit{\chi }}_{\mathrm{H}\mathrm{2}}\mathrm{)}\\ & & \mathrm{+}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathrm{l}}}{\mathit{n}}_{\mathit{i}}\hspace{0.17em}{\mathit{\chi }}_{\mathit{i}}\mathrm{+}{\mathit{e}}_{\mathrm{o}}\mathrm{\left(}{\mathit{n}}_{\mathrm{e}}\mathit{,T}\mathrm{\right)}\mathit{.}\end{array}$(6)Here, k, T, n_e, n_H2, n_l, n_i, n_o are Boltzmann’s constant, the gas temperature, the electron density, the H₂ molecule density, the number of levels in the hydrogen model atom, the density of atomic hydrogen in excitation or ionization state i and the number density of all other atoms and molecules that are not, or do not contain, hydrogen, respectively. The rotational and vibrational energy per H₂ molecule is e_H2, χ_H2 is the dissociation energy of H₂, χ_i is the excitation or ionization energy of atomic hydrogen and e_o is the internal energy of all other atoms and molecules that are not, or do not contain, hydrogen. The last quantity depends on the temperature and electron density and is computed in LTE.

Equation (6) is solved together with the equations for charge conservation, hydrogen nucleus conservation, evolution equations for the atomic hydrogen level populations n_i (5 bound levels plus the continuum): $\frac{\mathit{\partial }{\mathit{n}}_{\mathit{i}}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}{\mathit{n}}_{\mathit{i}}{u}\mathrm{\right)}\mathrm{=}\sum _{\mathit{j,j}\ne \mathit{i}}^{{\mathit{n}}_{\mathrm{l}}}{\mathit{n}}_{\mathit{j}}{\mathit{P}}_{\mathit{ji}}\mathrm{-}{\mathit{n}}_{\mathit{i}}\sum _{\mathit{j,j}\ne \mathit{i}}^{{\mathit{n}}_{\mathrm{l}}}{\mathit{P}}_{\mathit{ij}}\mathit{,}$(7)and an equation for non-equilibrium formation of H₂ molecules: $\frac{\mathit{\partial }{\mathit{n}}_{{\mathrm{H}}_{\mathrm{2}}}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}{\mathit{n}}_{{\mathrm{H}}_{\mathrm{2}}}{u}\mathrm{\right)}\mathrm{=}{\mathit{R}}_{\mathrm{f}}{\mathit{n}}_{\mathrm{1}}^{\mathrm{3}}\mathrm{-}{\mathit{R}}_{\mathrm{b}}{\mathit{n}}_{\mathrm{1}}{\mathit{n}}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{,}$(8)where n₁ is the ground state population of atomic hydrogen.

The radiative part of the rate coefficients P_ij is computed using the approximations given by Sollum (1999), the temperature-dependent rate coefficients R_f and R_b are taken from the UMIST database (Woodall et al. 2007, www.udfa.net).

These equations yield values for T, n_e, n_H2, and n_i. The gas pressure is then given by $\mathit{P}\mathrm{=}\mathit{kT}\left({\mathit{n}}_{\mathrm{e}}\mathrm{+}{\mathit{n}}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{+}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathrm{l}}}{\mathit{n}}_{\mathit{i}}\mathrm{+}{\mathit{n}}_{\mathrm{o}}\right)\mathit{,}$(9)and used in the momentum and energy equations (Eqs. (2), (3)). This non-equilibrium equation-of-state requires advection of 6 atomic and 1 molecular hydrogen population in addition to the 8 MHD variables, for a total of 15 advected quantities. See Leenaarts et al. (2007), Golding (2010) and Gudiksen et al. (2011) for more details.

Radiative heating.

Formally, the radiative flux divergence is given by $\mathrm{\nabla }\mathrm{·}{F}\mathrm{r}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathrm{\int }}_{\mathrm{\Omega }}\mathit{\alpha }\mathrm{\left(}\mathit{\nu ,}{n̂}\mathrm{\right)}\hspace{0.17em}\left(\mathit{S}\mathrm{\left(}\mathit{\nu ,}{n̂}\mathrm{\right)}\mathrm{-}\mathit{I}\mathrm{\left(}\mathit{\nu ,}{n̂}\mathrm{\right)}\right)\hspace{0.17em}\mathrm{d\Omega }\hspace{0.17em}\mathrm{d}\mathit{\nu ,}$(10)with S, α and I the source function, opacity and intensity at frequency ν in direction into the solid angle Ω. In practice, this double integral is too computationally expensive to compute, and various approximations are made.

The integral over frequency is replaced by four representative radiation bins, with their own associated bin-averaged opacity per mass unit κ_i, photon destruction probability ϵ_i and thermal emission E_i, assuming isotropic scattering and isotropic source function and opacity (Nordlund 1982; Skartlien 2000). The radiative heating of such a so-called multi-group scheme is then given by ${\mathit{Q}}_{\mathrm{rad}}\mathrm{=}\mathrm{4}\mathit{\pi \rho }\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{4}}{\mathit{\kappa }}_{\mathit{i}}\mathrm{(}{\mathit{ϵ}}_{\mathit{i}}{\mathit{J}}_{\mathit{i}}\mathrm{-}{\mathit{E}}_{\mathit{i}}\mathrm{)}\mathit{.}$(11)The implementation of this scheme in the Bifrost code is discussed in detail by Hayek et al. (2010). The quantities κ_i, ϵ_i and E_i are precomputed and tabulated assuming LTE populations and line-scattering in the 2-level Van Regemorter approximation. They are tabulated as function of the electron density and gas temperature. This is an improvement over previous work by Leenaarts et al. (2007) who used tables as function of the internal energy and the mass density, implicitly assuming LTE values for the electron density and the temperature.

The multi-group scheme has the drawback that to properly approximate the effects of the strongest spectral lines one needs to add several bins just to cover the largest opacities. Otherwise the cooling and heating effects of the few strong spectral lines in a bin otherwise comprised of spectral features with much lower opacity are washed out in the construction of the bin-averaged radiation quantities (see e.g., Eqs. (9)–(11) of Hayek et al. 2010). Furthermore, the assumptions of an opacity in LTE and a source function from the 2-level Van Regemorter approximation do not work for the strongest chromospheric lines. In the case of simulations of the chromosphere this means that the radiative heating and cooling from the mid-chromosphere upwards is badly modeled with a multi-group scheme only.

Therefore, the MHD model includes additional cooling in such lines in the following way:

Optically thin radiative losses from the corona and transition region are included through a frequency-integrated loss function Λ(T) based on the coronal approximation for the level populations: ${\mathit{Q}}_{\mathrm{thin}}\mathrm{=}\mathrm{-}\mathrm{\Lambda }\mathrm{\left(}\mathit{T}\mathrm{\right)}\hspace{0.17em}{\mathit{n}}_{\mathrm{e}}\hspace{0.17em}{\mathit{n}}_{\mathrm{H}}\mathit{.}$(12)The function Λ(T) is only significantly different from zero above T = 15   000 K.

Optically thick parameterized radiative losses and gains from H i, Mg ii and Ca ii are included through: ${\mathit{Q}}_{\mathrm{\left[}\mathrm{H}\mathit{,}\mathrm{Ca}\mathit{,}\mathrm{Mg}\mathrm{\right]}}\mathrm{=}\mathit{C}\mathrm{(}\mathit{T}{\mathrm{)}}_{\mathrm{\left[}\mathrm{H}\mathit{,}\mathrm{Ca}\mathit{,}\mathrm{Mg}\mathrm{\right]}}\hspace{0.17em}{\mathit{n}}_{\mathrm{e}}\hspace{0.17em}\mathit{\rho }\hspace{0.17em}\begin{array}{}\min \end{array}\left(\frac{{\mathrm{e}}^{\mathrm{-}{\mathit{k}}_{\mathrm{\left[}\mathrm{H}\mathit{,}\mathrm{Ca}\mathit{,}\mathrm{Mg}\mathrm{\right]}}{\mathit{m}}_{\mathrm{c}}}}{{\mathit{m}}_{\mathrm{c}}^{\mathrm{0.3}}}\mathit{,}\mathrm{1}\right)\mathit{.}$(13)Here the constant k and the temperature-dependent coefficient C are determined from detailed radiative transfer computations with the RADYN (see e.g., Carlsson & Stein 1995) and Multi3D codes (Leenaarts & Carlsson 2009), and m_c is the column mass. These functions include hydrogen lines and the Lyman continuum and the pertinent lines and continua from Ca ii and Mg ii (Carlsson & Leenaarts, in prep.).

Radiative heating through absorption of coronal radiation in UV continua is modeled through the representative bound-free absorption cross-section σ of He i at 25 nm: ${\mathit{Q}}_{\mathrm{He}}\mathrm{=}{\mathit{\sigma }}_{\mathrm{He}\mathit{,}\mathrm{25}\hspace{0.17em}\mathrm{nm}}\hspace{0.17em}{\mathit{n}}_{\mathrm{He} \mathrm{I}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{He}\mathit{,}\mathrm{25}\hspace{0.17em}\mathrm{nm}}}\hspace{0.17em}{\mathit{J}}_{\mathrm{thin}}$(14)with J_thin the angle-averaged radiation field resulting from Q_thin (Carlsson & Leenaarts, in prep.).

In addition we include an ad-hoc heating term that drives the temperature up to T₀ = 1660 K once it drops below that value. It is given by: ${\mathit{Q}}_{\mathrm{w}}\mathrm{=}{\mathit{C}}_{\mathrm{w}}{\mathit{n}}_{\mathrm{H}}^{\mathrm{2}}{\left[\begin{array}{}\max \end{array}\left(\mathrm{0}\mathit{,}\mathrm{(}{\mathit{T}}_{\mathrm{0}}\mathrm{-}\mathit{T}\mathrm{)}\right)\right]}^{\mathrm{2}}\mathit{,}$(15)with C_w a constant that sets the heating rate. The simulation thus still has an artificially set minimum temperature, but it is 740 K lower than in the simulation of Leenaarts et al. (2007). We add this ad-hoc term to prevent our simulation from cooling to too low temperatures where the radiation tables we employ become inaccurate.

All the different contributions are added so that the radiative heating term in Eq. (3) becomes $\mathrm{-}\mathrm{\nabla }\mathrm{·}{F}\mathrm{r}\mathrm{=}{\mathit{Q}}_{\mathrm{rad}}\mathrm{+}{\mathit{Q}}_{\mathrm{thin}}\mathrm{+}{\mathit{Q}}_{\mathrm{H}}\mathrm{+}{\mathit{Q}}_{\mathrm{Ca}}\mathrm{+}{\mathit{Q}}_{\mathrm{Mg}}\mathrm{+}{\mathit{Q}}_{\mathrm{He}}\mathrm{+}{\mathit{Q}}_{\mathrm{w}}\mathit{.}$(16)

Simulation setup.

We performed a 2D simulation with a grid size of 512 × 325, spanning from 1.5 Mm below ⟨τ₅₀₀⟩ = 1 to 14 Mm above it in the lower corona. The horizontal grid spacing is 32.5 km, the vertical grid spacing is 28 km from the convection zone up to the low corona, and increased to 150 km higher up in the corona.

The lower boundary is open, allowing fluid to freely enter and leave the box, while specifying the entropy of the inflowing gas to maintain sufficient energy flux into the computational domain. The upper boundary uses the methods of characteristics to extrapolate the MHD variables, letting waves exit the domain with almost no reflections. The upper boundary is set to strive to a temperature of 250 000 K to prevent the corona from cooling as a 2D weakly-magnetic simulation cannot sustain a corona self-consistently. We choose this rather low temperature because we model a weakly magnetic part of the atmosphere. Its exact value has no influence on the behaviour of the chromosphere. As formulated in the code, the thermal conduction operator requires a magnetic field to be present, so a weak magnetic field (average unsigned flux density in the photosphere of 0.3 G) is added, but is too weak to have any effect on the atmosphere.

The simulation was started from a previously relaxed snapshot computed without non-equilibrium hydrogen ionization and ran for 1 h of solar time. We discard the first 10 min to remove start-up effects.

Fig. 1 Snapshot of the simulated atmosphere. Top panel: gas temperature, with the column shown in Figs. 2, 3 and 7–12 indicated by a dashed line. Bottom panel:  − ∇·u, positive values indicate the gas is locally compressed in the co-moving frame.

Fig. 2 Properties of the atmosphere along the cut indicated in Fig. 1. a) Mass density. b) The solid black curve indicates the internal energy per mass unit, with the zero point at neutral atoms in the ground state. Red: ionization energy of hydrogen; black dotted: excitation energy of atomic hydrogen; blue: kinetic contribution of all atoms, molecules and electrons; green: contribution of the rotational, vibrational and dissociation energy of H2 molecules. c) Vertical gas velocity, positive means upward. d) Gas temperature. e) Gas pressure. f) Compression rate  − ∇·u.

3. Results

Figure 1 shows the chromosphere in a snapshot from the simulation after 34 min of solar time have elapsed. The upper panel shows the gas temperature, with granules visible below z = 0 Mm, the high-temperature transition region and corona in white at the top. In between lies the chromosphere, visible as a cool background state pervaded by waves and shock-fronts with peak temperatures increasing with height. The bottom panel shows  −∇·u, a measure for the local compression rate of the gas. Intergranular lanes and their extension into the upper photosphere appear as weakly compressing purple stalks at the bottom of the panel. The scene is dominated by strongly compressive shocks in the chromosphere and corona with expanding shock-wakes, such as the one around (x,z) = (13,2) Mm. Note that some of the shock fronts are co-spatial with the chromosphere-corona interface, and push the corona upward. An extreme example of this is seen along the dashed line at x = 13.1 Mm in the upper panel. The shock front at z = 3.7 Mm has pushed the corona up, leaving a large pocket of cool expanding gas in its wake. At z = 0.9 Mm a new shock front is propagating up through the cold gas, compressing it again.

Figure 2 shows various quantities along the cut indicated by the dashed line in Fig. 1. Panel a shows the density profile. It decreases with height everywhere except around z = 3.8 Mm at the site of the shock front that pushes the corona upward. Panel b shows the total internal energy in black with its distribution over various contributions. It has a peak at the shock front at z = 0.9 Mm and a smooth increase with height above it. The main contributors to the internal energy are the energy of the random motions of the gas particles (blue) and the ionization energy of hydrogen (red). The latter is the largest contribution above 1.6 Mm. There is a small but important contribution of H₂ molecules at 0.5 Mm height. Panel c shows the vertical velocity. It exhibits a rough sawtooth shape common to shock-propagation in a stratified atmosphere. Panel d shows the temperature, with a high-temperature shock front at 0.9 Mm, and a plateau at 1660 K between 1 and 2.5 Mm. Panel e shows the gas pressure and panel f shows the compression rate, again indicating the presence of the two shock fronts at 0.9 and 3.8 Mm, weak compression at 2.8 Mm and expansion in the rest of the chromosphere.

Fig. 3 Energy balance along the cut indicated in Fig. 1. The upper part shows heating, the lower part cooling. Black: total heating rate (right-hand side of Eq. (3)), excluding the ad hoc term QW. Green: compression work  − P(∇·u). Red: total radiative heating (Eq. (16)), excluding QW. Blue: ad-hoc heating QW (Eq. (15)).

Figure 3 displays the energy balance in the chromosphere (the right-hand-side of Eq. (3)) for the same column as Fig. 2. Viscous energy dissipation, Joule heating and thermal conduction are all negligible in the chromosphere and are not shown. The energy balance is set by compression work and radiation. As expected, the shock fronts are heated compressively and cool radiatively. The cool area between 1 and 3.4 Mm is cooling through expansion and is heated by radiation except at 2.7 Mm, where a horizontally propagating wave causes some compression. The expansion cooling is stronger than the radiation heating by an order of magnitude. This imbalance ultimately leads to the ad-hoc Q_W contribution becoming active to prevent the atmosphere from cooling below 1660 K. The instant in time shown in the figure shows a scene with large ad-hoc heating. However, this term is only intermittently active and switches off as soon as the temperature is above 1660 K.

Fig. 4 Time evolution of the temperature at different heights in the column indicated by the dashed line in Fig. 1. The time moment shown in that figure is indicated by the long tick marks at t = 8.4 min. Blue: z = 1 Mm; green: 1.5 Mm; red: 2 Mm.

Figure 4 shows the time evolution of the temperature at different heights along the dashed line of Fig. 1. The blue curve between t = 16 and t = 24 min show a clean example of the temperature variations of passing shocks, with a rapid increase in temperature and a more gradual cooling phase after passage of the shock front. In general the time evolution is more irregular because of slanted shock propagation and interference. The corona intermittently dips down to below 1.5 Mm, indicated by the red and green curves running off the temperature scale.

Fig. 5 Diagram of the temperature occurrence as function of atmospheric height in the simulation. The red curve shows the temperature in the column marked in Fig. 1. The blue curves indicate the minimum and maximum temperatures as a function of height during the whole simulation run. The black curves shows the location of 50% He i (top curve) and H (bottom curve) ionization assuming Saha equilibrium and the average run of the electron density with height.

The above discussion and Figs. 2–4 again confirm the picture of the quiet chromosphere as a layer that undergoes quick compression during shock passages followed by a longer phase of expansion cooling. The time scale of this expansion cooling t_PdV = e/(P(∇·u)) varies from 200 s to 500 s assuming typical values of e = 10¹² erg g^-1 and 2 × 10⁹ < P   (∇·u) < 5 × 10⁹ erg s^-1 g^-1. Here e is roughly the thermal energy of solar gas at 12   000 K. The ionization energy remains nearly constant, so the gas in the chromosphere behaves approximately as an ideal gas. So, expansion cooling alone can, in cases of strong expansion, cool chromospheric gas to very low temperatures between 2 shock passages, as the radiative heating rarely exceeds 2 × 10⁸ erg s^-1 g^-1 in the chromosphere.

This rough estimate is confirmed by Fig. 5. It shows a histogram of the occurrence of gas temperature values as a function of atmospheric height. The range of temperatures increases with increasing height when going up from the photosphere. The maximum temperature curve shows the increase of peak shock temperature with height up to 1 Mm. The maximum temperature at 1 Mm is above 15   000 K, indicating that the corona occasionally reaches this far down. The thick dark band at 10   000 K between 1 and 4 Mm is caused by a combination of shock fronts and the layer of 10   000 K gas just below the corona (see the upper panel of Fig. 1). This band is discussed in Sect. 3.2. Between z = 0.8 and 2.2 Mm the histogram peaks along a narrow dark band below T ≈ 2000 K, indicating that such low temperatures are common in the simulated chromosphere. Above 2.2 Mm such low temperatures occur less frequently, but the atmosphere can be as cold as the ad-hoc heating threshold temperature of 1660 K up to 3.4 Mm height.

In the following subsections we investigate the accuracy of this result by discussing the various physical processes that determine the minimum temperature in the chromosphere and the accuracy with which they are represented in the simulation.

Fig. 6 Diagram of the occurrence of temperature as function of mass density in the simulation. The red contours show the fraction of hydrogen atoms bound in H2 molecules, 2nH2/(nH   i + 2nH2), as function of temperature and mass density, assuming ICE and all hydrogen neutral. The contours are spaced a factor 100 apart, with the labels indicating the exponent a in 10a. The blue line at T = 1.66 kK specifies the threshold for the ad hoc heating (Eq. (15)).

3.1. H₂ thermostat action

When the low-chromospheric temperature reaches down to 2000–3000 K, H₂ molecules can form in significant amounts. This process is exothermic, releasing 4.48 eV per molecule formed, compared to a thermal energy of kT = 0.19 eV per particle at 2000 K. Thus, as H₂ begins to form it releases a large amount of energy, which by Eq. (6) is predominantly converted to thermal energy, counteracting the expansion cooling. This effect is illustrated in Fig. 6. It shows the simultaneous occurrence rate of combinations of temperature and mass density in the simulation, with overplotted contours of the fraction of hydrogen atoms bound in H₂ assuming instantaneous chemical equilibrium.

The minimum temperature in the chromosphere follows the curve that represents 10% of all hydrogen bound in H₂ between ρ = 10^-9 g cm^-3 and 10^-7.5 g cm^-3. At lower mass densities (higher up in the chromosphere) the formation rate of H₂ becomes so low that it cannot form in large enough amount to prevent the expanding parts of the atmosphere from cooling further. This is indicated by the turnoff of the bottom of the grey cloud away from the red curve at ρ = 10^-9 g cm^-3. Once the temperature drops below 1660 K the ad-hoc heating becomes active, preventing the atmosphere from cooling even further.

3.2. H and He as thermostats

Figure 1 shows a layer of 10   000 K hanging as a ragged skirt below the corona. This layer is one of the causes of the dark band at the same temperature in Fig. 5. The overlaid He  ionization curve demonstrates that He i works similarly to H₂ as a thermostat in this simulation. This is unphysical, however, since the actual ionization-recombination balancing is likely to be even more out of instantaneous equilibrium for helium than in the case of hydrogen. Thus, the Saha equilibrium assumed in our simulation for helium is erroneous. Relaxing this assumption is likely to remove helium’s thermostat action, just as for non-equilibrium hydrogen ionization (Carlsson & Stein 2002; Leenaarts et al. 2007). Indeed, hydrogen does not cause any such thermostat clustering of temperature values (which would be at the lower black curve in Fig. 5) since it is treated properly in non-equilibrium. The erroneous assumption of LTE balancing for He i does not affect the present analysis, however, since it affects only the hottest phases of the chromospheric gas and not the coolest ones addressed here.

Fig. 7 Various contributions to the total radiative heating rate (Eq. (16)) for the column indicated in Fig. 1. Black solid: total radiative heating. Black dashed: radiative heating from the multi-group scheme (Qrad, see Eq. (11)). Red: QH. Green: QCa. Blue: Qthin. Red dashed: QHe. Grey with black dots: QMg.

3.3. Radiative heating

The chromosphere is heated by radiation during its cool phases. Fig. 7 shows the various contributions that make up the total radiative heating rate in our simulation along the column indicated in Fig. 1.

The radiative heating due to the multi-group scheme (Eq. (11)) dominates the heating below z = 1 Mm as indicated by the near-equality of the solid and dashed black curves. It also cools strongly in the upper chromosphere. Absorption of coronal radiation (red dashed curve, Eq. (14)) adds heating above 1.3 Mm and is dominant between 1.8 and 3.2 Mm. The Ca ii lines (green, Eq. (13)) cool in the shock front at 0.9 Mm and above 2.8 Mm, they are the dominant contribution in the cool gas between 1 and 1.8 Mm. Hydrogen (red, Eq. (13)) cools in the shock front at 0.9 Mm and above 2.8 Mm, where the Lyman lines and continua become effectively optically thin. Optically thin losses from the corona (blue) do not play a role in the chromosphere. The Mg ii (grey-black dotted) lines only cool significantly in the upper chromosphere.

We will now discuss the different contributions in detail and show them for the representative column from our simulation displayed in Figs. 2–4.

Fig. 8 Comparison of the multi-group radiative heating rate in the first bin (black solid, Q1 see Eq. (11)) compared to the heating rate due to H− in NLTE (red) and LTE (green).

Negative hydrogen ion H⁻.

Absorption of radiation in the H⁻ continuum is a potentially large source of heating. This absorption is included in the first radiation bin of our model. We computed the H⁻ heating in detail for the snapshot of Fig. 1 assuming LTE, i.e., the source function is the Planck function and the H⁻ density is set by its Saha-equilibrium with neutral hydrogen. For comparison, we computed the H⁻heating rate in NLTE including the following reactions: $\begin{array}{ccc}\mathrm{H}\mathrm{-}\mathrm{+}\mathit{h\nu }& \mathrm{\leftrightarrow }& \mathrm{H}\mathrm{+}\mathrm{e}\mathit{,}\\ \mathrm{H}\mathrm{-}\mathrm{+}\mathrm{H}& \mathrm{\leftrightarrow }& \mathrm{2}\mathrm{H}\mathrm{+}\mathrm{e}\mathit{,}\\ \mathrm{H}\mathrm{-}\mathrm{+}\mathrm{e}& \mathrm{\leftrightarrow }& \mathrm{H}\mathrm{+}\mathrm{2}\mathrm{e}\mathit{,}\\ \mathrm{H}\mathrm{-}\mathrm{+}\mathrm{p}& \mathrm{\leftrightarrow }& \mathrm{2}\mathrm{H}\mathit{,}\\ \mathrm{H}\mathrm{-}\mathrm{+}\mathit{H}& \mathrm{\leftrightarrow }& {\mathrm{H}}_{\mathrm{2}}\mathrm{+}\mathrm{e}\mathit{,}\end{array}$where we kept the neutral hydrogen, proton, electron and, H₂ density constant. This constancy is justified because of the high number density of these particles relative to the H⁻density. The reaction rates for the latter four reactions are from Lambert & Pagel (1968).

Figure 8 shows a comparison with the first bin of the multi-group heating rate (black, Eq. (11)) and the heating rate caused by H⁻bound-free and free-free transitions, assuming LTE (green) and NLTE (red). The LTE assumption generates the largest heating, caused by both the low Planck function compared to the average radiation field, and the high H⁻number density. The NLTE heating rate is much smaller. This is caused by a combination of photo-ionization of H⁻, leading to a NLTE departure coefficient much smaller than unity, and strong scattering, leading to smaller J_ν − S_ν splitting. Absorption of radiation by H⁻ is therefore not efficient in heating the chromosphere.

The first bin in the multi-group scheme closely matches the heating due to H⁻ where H⁻ dominates the opacity (up to 0.8 Mm). It deviates in the shock around 0.9 Mm and above. The first bin underestimates the NLTE H⁻ heating in the cool area between 1 and 2.5 Mm. Above 3 Mm it erroneously overestimates the NLTE cooling by several orders of magnitude.

Fig. 9 Comparison of the parameterized radiative heating rate due to Ca ii (QCa, see Eq. (13)) to the heating computed based on a detailed radiative transfer computation. Black: QCa. Blue: heating due to Ca ii lines. Green: heating due to Ca ii continua. Red: total Ca ii heating. The blue curve is nearly equal to the red curve.

Ca II.

The lines of Ca ii are other strong heating agents in the chromosphere. We computed the Ca ii heating rates for the snapshot of Fig. 1 using the radiative transfer code Multi3d. This detailed computation was performed treating each column in the simulation as a plane-parallel atmosphere and included the effects of partial redistribution (PRD) in the H & K lines. The results for the representative column are shown in Fig. 9. The Ca ii continua (green curve) have a negligible effect on the total heating rate. The lines cool strongly in the shock front at 0.8 Mm and above 3.4 Mm. They heat the atmosphere in the cool area between 0.9 and 3.4 Mm.

The parameterized Ca ii cooling employed in our simulation is in reasonable agreement with the detailed computation above 0.5 Mm height. The largest differences are the overestimation of the cooling between 3.5 and 4.2 Mm and shifted location around 3 Mm where the heating term switches sign. The difference below 0.5 Mm is caused by a cutoff to avoid doubling the heating rate as heating at low heights is also included in the multi-group scheme.

Manual inspection of the heating rate in many different columns of the simulation shows that the parameterized heating is accurate most of the time.

Fig. 10 Comparison of the parameterized radiative heating rate due to Mg ii (QMg, see Eq. (13)) to the losses computed based on a detailed radiative transfer computation. Black: QMg. Blue: heating due to Mg ii lines. Green: heating due to Mg ii continua. Red: total Mg ii heating. The blue curve is nearly equal to the red curve.

Mg II.

Another potentially large contributor to the total chromospheric heating is Mg ii. We computed the radiative losses in detail including the effects of PRD in the h&k lines. The result for the representative column is given in Fig. 10. The h&k lines cool in the shock front at 0.9 Mm, heat the cool area above the shock front and cool the upper chromosphere. The continua play a minor role.

The cutoff in the parameterized heating pushes the heating to zero at 1.5 Mm. The heating is much lower than the detailed computation in the cool area between 1.5 and 3 Mm and it shows the same shift as for Ca ii in the location where the heating switches sign. The cooling peak around 3.8 Mm is reproduced correctly by Q_Mg.

Manual inspection shows that the parameterized Mg ii heating is accurate for gas temperatures above 6000 K. Below this temperature the heating rate is modeled qualitatively correct, but can deviate because the effects of the precise atmospheric structure, in particular the velocity field are not taken into account in the simple parameterization.

Hydrogen.

Hydrogen has only a small effect on the heating rate of the lower and middle chromosphere in the semi-empirical time-independent VAL C model atmosphere (Vernazza et al. 1981). In this model the Lyα transition is so optically thick that it does not play a role in the energy balance, and the losses in Hα are approximately offset by the absorption in the Balmer continuum. This effect remains valid in the time-dependent models of Carlsson & Stein (2002) computed with the RADYN code, except in shock fronts. Unfortunately it is currently not possible to compute the detailed time-dependent heating rate due to hydrogen in multi-dimensional simulations, so we cannot rigorously prove the same effect holds in our 2D simulation. However, the small hydrogen heating rate is mostly an effect of the atomic structure of hydrogen, and does not depend strongly on the exact chromospheric structure. Therefore, our 2D simulation should behave in a similar manner. This is confirmed in Fig. 11.

Fig. 11 Comparison of the parameterized radiative heating rate due to hydrogen (QH, see Eq. (13)) to the losses implied by the non-equilibrium hydrogen radiative rates. Black: QH. Blue: implied heating due to hydrogen lines. Green: implied heating due to hydrogen continua. Red: total implied heating.

This figure shows the parameterized heating due to hydrogen (Eq. (13)) in black. The shock front at 0.9 Mm is cooling, there is little heating between 1.2 and 2.5 Mm, and slowly increasing cooling above 2.8 Mm, peaking at 4.1 Mm where the Lyα line becomes optically thin.

Figure 11 also shows the heating implied by the radiative transition rates computed from Eq. (7). We define the implied heating rate in a transition between levels i and j (i < j) by ${\mathit{Q}}_{\mathrm{imp}}\mathrm{=}\mathit{h}{\mathit{\nu }}_{\mathrm{0}}\mathrm{(}{\mathit{R}}_{\mathit{ij}}\mathrm{-}{\mathit{R}}_{\mathit{ji}}\mathrm{)}\mathit{,}$(17)with ν₀ the line center or ionization edge frequency, and R_ij the radiative rate from level i to level j. This simple definition makes an error for bound-free transitions due to the non-negligible width of the ionization edges, but this error is small. Note that the implied heating does not appear as a source term in Eq. (16), the assumed radiation field only serves to define the radiative rates in Eq. (7). It reproduces the near-zero heating in the cool area between 1 Mm and 2.8 Mm. The shock-front and the upper chromosphere are heated.

We conclude that both the parameterized cooling and the implied heating rate reproduce the lack of hydrogen heating in the cool phases of the chromosphere. Any errors caused be the approximations have negligible effect on the minimum temperature in our model due to the small heating rate compared to other heating agents.

Absorption of coronal radiation.

Nearly all coronal radiation emitted towards the chromosphere is eventually absorbed (a small part is backscattered into space). Exactly where this radiation is absorbed depends on its detailed spectral energy distribution and the absorption coefficient in the chromosphere. However, these details are of minor importance to the overall heating rate as long as all energy is absorbed. Our simulation correctly models this behavior. The finite extent of the corona in the simulation and its relatively low temperature might result in a too low amount of radiation absorbed in the chromosphere. The coronal loss function (Q_thin, Eq. (12)) peaks strongly in the lower corona due to its quadratic density dependence, so we expect this error to be small.

4. Discussion and conclusions

Robustness of the minimum temperature

The formation of H₂ acts as a thermostat in the lower chromosphere, preventing the temperature from dropping below 2000 K. ICE is not valid in the middle and upper chromosphere. There, the slow formation rate of H₂ prevents thermostat action, allowing solar gas to cool well below 2000 K. Our simulation correctly models this behavior.

In Sect. 3.3 we discussed the various lines and bound-free transitions that can radiatively heat the chromosphere. We showed heating in lines of Ca ii and absorption of coronal radiation are the dominant processes. We showed that H⁻ is inefficient in heating cool pockets of chromospheric gas owing to strong scattering and low H⁻population due to photo-ionization. Mg ii and hydrogen play a minor role in the heating of the chromosphere. Our simulation accurately models the heating due to Ca ii and coronal radiation. The other processes are less accurately modeled but these have, in sum, only a small effect on the energy balance.

The heating due to acoustic waves is self-consistently included in the simulation since we include the upper convection zone and the stochastic excitation of acoustic waves there. The limited resolution may lead to an underestimate of the excitation of high-frequency waves. A wave at a frequency of 20 mHz is represented with 11 grid-points so the effect will only be important for really high frequency waves where simulations and observations indicate the effect for the heating of the chromosphere is minimal (Fossum & Carlsson 2005; Carlsson et al. 2007).

The ad-hoc heating term acts intermittently in the largest expansion bubbles, thus keeping our minimum chromospheric temperature at 1660 K.

We do not expect that a 3D simulation with similar physics will change this result. Three-dimensional simulations, such as the one reported on by Leenaarts et al. (2010) show temperatures down to 2000 K even with LTE hydrogen ionization and instantaneous H₂ formation. It is very likely that such 3D simulations with a non-equilibrium EOS would show the same low temperatures as reported here (see Fig. 4 of Leenaarts et al. 2007, for an illustration of the resulting temperature differences between these equations-of-state).

We therefore conclude that our simulation provides an upper bound to the minimum temperature of the non-magnetic chromosphere. The simulation predicts the occurrence of gas temperatures as low as 1660 K up to 3.4 Mm height. This low temperature is common between 0.7 and 2.5 Mm height. Our ad-hoc heating prevents our model from cooling further. This means that a non-magnetic chromosphere can have temperatures even lower than this value. It seems impossible to avoid such low temperatures using only radiation and hydrodynamic processes.

Is the chromosphere really this cool?

The sun has a turbulent photospheric magnetic field of the order of 130 G (Trujillo Bueno et al. 2004). This turbulent field is likely to extend into the quiet chromosphere and might generate enough Joule heating to prevent the quiet chromosphere from cooling down to temperatures as low as in our simulation. In addition the presence of a relatively strong magnetic field will also change the propagation properties of the generated acoustic waves, reducing their amplitudes by confining them to flux tubes and/or by conversion to other magneto-hydrodynamic modes in the vicinity of the β = 1 surface.

The simulations of a solar surface dynamo by Vögler & Schüssler (2007) and Schüssler & Vögler (2008) suggest that the strength of the turbulent field declines rapidly with height, which would limit the amount of mid and high chromospheric heating the field could cause. On the other hand, the work of Schrijver & Title (2003) suggests that part of the network magnetic fields can connect back in the internetwork photosphere. Such a multiscale magnetic carpet would lead to an increase of the field strength in the quiet chromosphere relative to a purely locally generated field, with a corresponding increase in the potential for Joule heating.

Magnetic fields in the low ionization-degree chromosphere may also lead to the generation of electric currents through neutral-ion drag (Krasnoselskikh et al. 2010).

Our simulation serves as a baseline case to test the effect of such magnetic heating processes. It should be compared against 3D models like the one by Vögler & Schüssler (2007) but extended up into the corona, and including all the physics discussed in this paper to study the effect of the turbulent field. Such a simulation does not require a large horizontal extent of the computational domain and, while computationally very demanding, can in principle be done on the largest currently available supercomputers.

Fig. 12 Molecular densities in the exemplary column indicated by a dashed line in Fig. 1. The densities of CO (green), OH (blue) and TiO (red) are based on instantaneous chemical equilibrium. The H2 density (black dashed) is based on our EOS. The total number density of hydrogen atoms is plotted in solid black for comparison.

Simulation of the effect of the magnetic carpet requires a larger horizontal extent to include some network magnetic field and a large enough area of quiet sun. The high spatial resolution (≈ 10 km grid spacing) needed to get local dynamo action combined with the large size of a network cell (≈ 30 Mm) makes such a simulation challenging.

At this moment the minimum temperature of the quiet sun chromosphere remains an open question. Our simulation shows a non-magnetic chromosphere will get colder than our ad-hoc limit of 1660 K, but it is unknown whether regions uninfluenced by magnetic fields occur in the chromosphere. However, if these low temperature areas exist, they can in principle be observed.

Suggestions for observations

Direct observation of cool pockets of chromospheric gas is hard. Atomic spectral lines with sufficient opacity form generally in NLTE, making it hard to derive temperatures. The Atacama Large Millimeter Array (ALMA) would be the ideal instrument to probe chromospheric temperatures because of its high resolution, the LTE source function and relatively well-defined formation height range of the sub-millimeter continua (Loukitcheva et al. 2004; Wedemeyer-Böhm et al. 2007).

Alternatively, off-limb spectroscopy in molecular lines with sub-arc-second resolution and low stray light could be used. As an example we show the number densities of some abundant molecules as a function of height in our exemplary column in Fig. 12. Because the chemical reaction timescales become very large in the chromosphere the number densities for the species computed assuming ICE should be taken as a maximum possible value only, and are likely orders of magnitude too high. The molecules of choice for such observations seem to be H₂ and CO. The solar spectrum shows H₂ UV emission lines (Jordan et al. 1978) that can be observed by the SUMER instrument aboard the Solar and Heliospheric Observatory (SOHO) spacecraft (Kuhn & Morgan 2006). The mere presence of molecular lines at several arc-seconds above the limb would prove the existence of low temperature gas at high-chromospheric heights. Unfortunately the reverse is not true, the absence of lines might merely indicate the absence of molecules and not the absence of cool gas.

Conclusions.

We have performed a 2D radiation-MHD simulation of the non-magnetic solar chromosphere, including the effects of non-equilibrium ionization of hydrogen and non-equilibrium formation of H₂ molecules and NLTE radiative transfer. We find that our simulation contains pockets of cool gas with temperatures down to 1660 K from 0.8 Mm up to 3.4 Mm height. We discuss the physical mechanisms that set the minimum temperature and find that our simulation likely overestimates the minimum temperature. We conclude it is impossible to avoid such low temperatures with radiation-hydrodynamical processes only.

Such cool pockets of chromospheric gas might be even cooler in the real sun than in our simulation, provided a quiet chromosphere without significant magnetic heating exists. We suggest off-limb molecular spectroscopy to look for such pockets, and suggest 3D simulations with sufficient resolution to support a local dynamo with and without network-like magnetic fields to investigate whether the assumption of negligible magnetic heating is justified.

Acknowledgments

We thank N. Vitas for compiling data on molecular partition functions and chemical equilibrium constants. We thank R. J. Rutten for illuminating discussions and improvements to the manuscript. J.L. recognizes support from the Netherlands Organization for Scientific Research (NWO). This research was supported by the Research Council of Norway through the grant “Solar Atmospheric Modelling” and through grants of computing time from the Programme for Supercomputing.

References

All Figures

	Fig. 1 Snapshot of the simulated atmosphere. Top panel: gas temperature, with the column shown in Figs. 2, 3 and 7–12 indicated by a dashed line. Bottom panel:  − ∇·u, positive values indicate the gas is locally compressed in the co-moving frame.
In the text

Fig. 2 Properties of the atmosphere along the cut indicated in Fig. 1. a) Mass density. b) The solid black curve indicates the internal energy per mass unit, with the zero point at neutral atoms in the ground state. Red: ionization energy of hydrogen; black dotted: excitation energy of atomic hydrogen; blue: kinetic contribution of all atoms, molecules and electrons; green: contribution of the rotational, vibrational and dissociation energy of H2 molecules. c) Vertical gas velocity, positive means upward. d) Gas temperature. e) Gas pressure. f) Compression rate  − ∇·u.

In the text

	Fig. 3 Energy balance along the cut indicated in Fig. 1. The upper part shows heating, the lower part cooling. Black: total heating rate (right-hand side of Eq. (3)), excluding the ad hoc term QW. Green: compression work  − P(∇·u). Red: total radiative heating (Eq. (16)), excluding QW. Blue: ad-hoc heating QW (Eq. (15)).
In the text

	Fig. 4 Time evolution of the temperature at different heights in the column indicated by the dashed line in Fig. 1. The time moment shown in that figure is indicated by the long tick marks at t = 8.4 min. Blue: z = 1 Mm; green: 1.5 Mm; red: 2 Mm.
In the text

Fig. 5 Diagram of the temperature occurrence as function of atmospheric height in the simulation. The red curve shows the temperature in the column marked in Fig. 1. The blue curves indicate the minimum and maximum temperatures as a function of height during the whole simulation run. The black curves shows the location of 50% He i (top curve) and H (bottom curve) ionization assuming Saha equilibrium and the average run of the electron density with height.

In the text

Fig. 6 Diagram of the occurrence of temperature as function of mass density in the simulation. The red contours show the fraction of hydrogen atoms bound in H2 molecules, 2nH2/(nH   i + 2nH2), as function of temperature and mass density, assuming ICE and all hydrogen neutral. The contours are spaced a factor 100 apart, with the labels indicating the exponent a in 10a. The blue line at T = 1.66 kK specifies the threshold for the ad hoc heating (Eq. (15)).

In the text

	Fig. 7 Various contributions to the total radiative heating rate (Eq. (16)) for the column indicated in Fig. 1. Black solid: total radiative heating. Black dashed: radiative heating from the multi-group scheme (Qrad, see Eq. (11)). Red: QH. Green: QCa. Blue: Qthin. Red dashed: QHe. Grey with black dots: QMg.
In the text

	Fig. 8 Comparison of the multi-group radiative heating rate in the first bin (black solid, Q1 see Eq. (11)) compared to the heating rate due to H− in NLTE (red) and LTE (green).
In the text

	Fig. 9 Comparison of the parameterized radiative heating rate due to Ca ii (QCa, see Eq. (13)) to the heating computed based on a detailed radiative transfer computation. Black: QCa. Blue: heating due to Ca ii lines. Green: heating due to Ca ii continua. Red: total Ca ii heating. The blue curve is nearly equal to the red curve.
In the text

	Fig. 10 Comparison of the parameterized radiative heating rate due to Mg ii (QMg, see Eq. (13)) to the losses computed based on a detailed radiative transfer computation. Black: QMg. Blue: heating due to Mg ii lines. Green: heating due to Mg ii continua. Red: total Mg ii heating. The blue curve is nearly equal to the red curve.
In the text

	Fig. 11 Comparison of the parameterized radiative heating rate due to hydrogen (QH, see Eq. (13)) to the losses implied by the non-equilibrium hydrogen radiative rates. Black: QH. Blue: implied heating due to hydrogen lines. Green: implied heating due to hydrogen continua. Red: total implied heating.
In the text

	Fig. 12 Molecular densities in the exemplary column indicated by a dashed line in Fig. 1. The densities of CO (green), OH (blue) and TiO (red) are based on instantaneous chemical equilibrium. The H2 density (black dashed) is based on our EOS. The total number density of hydrogen atoms is plotted in solid black for comparison.
In the text