Radiation hydrodynamics in simulations of the solar atmosphere

Abstract

Nearly all energy generated by fusion in the solar core is ultimately radiated away into space in the solar atmosphere, while the remaining energy is carried away in the form of neutrinos. The exchange of energy between the solar gas and the radiation field is thus an essential ingredient of atmospheric modeling. The equations describing these interactions are known, but their solution is so computationally expensive that they can only be solved in approximate form in multi-dimensional radiation-MHD modeling. In this review, I discuss the most commonly used approximations for energy exchange between gas and radiation in the photosphere, chromosphere, and corona.

1 Introduction

The interaction of matter and radiation is an indispensable ingredient of solar and stellar atmospheric modeling. Around 98% of the energy generated in the solar core is transported outwards first by radiative diffusion, then convection, and ultimately escapes into space in the photosphere of the Sun where the overlying solar material becomes transparent. The remainder of the energy escapes the Sun in the form of neutrinos (Turck-Chièze and Couvidat 2011). The layers above the photosphere (chromosphere and corona) are hotter than radiative equilibrium models predict. Therefore deposition of non-thermal energy and conversion into heat must occur in these layers. The radiative energy losses (ignoring the $\sim 10^{-6}$ fraction of energy carried away by the solar wind, Le Chat et al. 2012) from the chromosphere and corona must balance this energy deposition in a time and space-averaged sense. Radiation is thus an essential ingredient in setting the structure of the outer solar atmosphere and the response of the atmosphere to non-thermal energy deposition.

Modeling the interaction of radiation and matter in the solar atmosphere is in general a difficult problem. The specific intensity, the fundamental quantity used to characterise the radiation field, depends on seven parameters: three space dimensions, two angles describing direction, frequency, and time. In addition, the radiative transfer problem is non-linear and non-local: local changes to the intensity through absorption and emission depend on the intensity itself, and radiation emitted at one place in the atmosphere can influence other locations.

Writing general equations that describe the interaction of radiation and matter is not so difficult. Because of limitations in computing time these equations have so far been solved rather completely in one-dimensional geometry only. The Radyn code (Carlsson and Stein 1992, 2002), is perhaps the most well-known example, but other codes are used too (e.g., Flarix, see Kašparová et al. 2009).

However, solving these general equations in 2D and 3D is still beyond current computational capabilities. This is problematic because modeling convection requires at least 2D geometry, and fully modeling the rich physics caused by the interaction of the magnetic field with matter requires 3D geometry.

Since the pioneering simulations of Nordlund (1982), a plethora of codes that aim to model solar and/or stellar atmospheres in 3D have been developed. The ones that I am aware of are: Stagger (e.g., Stein and Nordlund 1998, but many different versions of this code exist), MURaM (Vögler 2004; Rempel 2017), Bifrost (Gudiksen et al. 2011), CO5BOLD (Freytag et al. 2012), StellarBox (Wray et al. 2015), Ramens (Iijima and Yokoyama 2015; Iijima 2016), Mancha3D (Khomenko et al. 2017, 2018), Antares (Leitner et al. 2017), and RadMHD (Abbett and Fisher 2012). The latter is the only code that does not use the radiation treatment developed by Nordlund (1982), and instead uses a simpler but much faster method.

Radiation-hydrodynamics in the solar atmosphere is a vast topic. In this review, I limit myself to describing only the most commonly used approximations for radiative transfer in the photosphere, chromosphere, and corona in these multi-dimensional radiation-MHD codes. I do not discuss results obtained with any of these codes. An excellent review of solar magnetoconvection as studied using radiation-MHD simulations is given in Stein (2012). The field of radiation-MHD modeling of the combined photosphere, chromosphere, transition region, and/or corona has developed tremendously during the last 15 years. To my knowledge no recent review covering this development exists. Example starting points for studying applications are Carlsson et al. (2016), Martínez-Sykora et al. (2017), Khomenko et al. (2018), and Cheung et al. (2019).

In the convection zone the radiation diffusion approximation holds to a high degree of precision and the energy exchange between radiation and the solar gas is close to zero. In the photosphere the gas loses large amounts of energy in the form of radiation which then largely escapes into space. The LTE assumption for the source function and opacity is still not too bad and even the grey approximation is still reasonably accurate. I devote a large fraction of this review discussing approximations for the photosphere and low chromosphere in Sect. 3.

The situation in the chromosphere and transition region is more complex. The chromosphere is optically thin for optical continuum radiation and most radiative energy exchange takes place in a few strong spectral lines. Radiation scattering is important so that non-LTE effects must be taken into account, the ionisation balance of hydrogen and helium is out of equilibrium so that assuming LTE or statistical equilibrium to compute opacities or source functions is in general no longer accurate. Modeling energy exchange taking into account these complexities is discussed in Sects. 5 and 6.

In the corona the physics of radiative energy losses and gains becomes again somewhat less complex. Most coronal structures are optically thin for all frequencies except in the radio regime. The radio regime lies in the far tail of the Planck function and does not contribute significantly to the radiative losses. For all other frequencies it can be assumed that photon absorption does not take place and that all photons emitted by the gas escape from the corona. It is discussed in Sect. 4.

2 Fundamentals

Excellent books that review the fundamentals of radiation hydrodynamics and modeling of stellar atmospheres are Mihalas and Mihalas (1984) and Hubeny and Mihalas (2014). The first book focusses on fundamental theory, while the second one discusses modeling of 1D static and moving atmospheres. Below, I briefly touch upon some general aspects that are relevant for the methods discussed in Sect. 3–7.

2.1 The MHD equations including radiation

The dominant paradigm for multidimensional modeling of the solar atmosphere has been the magnetohydrodynamics (MHD) approximation. The MHD equations for the density, momentum and internal energy including radiation terms can be written as

$$\begin{aligned} \frac{\partial \rho }{\partial t} = - {{\varvec{\nabla }}}\cdot \left( \rho {\mathbf{v}} \right) , \end{aligned}$$

(1)

$$\begin{aligned} \frac{\partial {\mathbf{p}} }{\partial t} = - {{\varvec{\nabla }}}\cdot \left( {\mathbf{v}} \otimes {\mathbf{p}} - {{\varvec{\tau }}} \right) - {{\varvec{\nabla }}}P + {\mathbf{J}} \times {\mathbf{B}} + \rho {\mathbf{g}} - {{\varvec{\nabla }}}\mathbf{P} _{{{\rm rad}}}, \end{aligned}$$

(2)

$$\begin{aligned} \frac{\partial e}{\partial t} = - {{\varvec{\nabla }}}\cdot \left( e {\mathbf{v}} \right) - P {{\varvec{\nabla }}}\cdot {\mathbf{v}} + Q + Q_{{\rm rad}}, \end{aligned}$$

(3)

with $\rho$ the mass density, ${\mathbf{v}}$ the velocity vector, ${\mathbf{p}}$ the momentum density, ${{\varvec{\tau }}}$ the stress tensor, P the gas pressure, ${\mathbf{J}}$ the current density, ${\mathbf{B}}$ the magnetic field vector, ${\mathbf{g}}$ the acceleration due to gravity, $\mathbf{P} _{{{\rm rad}}}$ the radiation pressure tensor, e the internal energy, $Q_{{\rm rad}}$ the heating or cooling owing to radiation, while Q expresses energy exchange by any other processes such as dissipation of currents, heat conduction, and viscosity.

The assumptions under which MHD is valid tend to be fulfilled in the photosphere and convection zone as well as the corona, but break down in the chromosphere and transition region where the frequencies of collisions between particles become smaller than the cyclotron frequencies of ions and electrons (e.g., Khomenko et al. 2014).

Several radiation-MHD codes have been extended to include some effects beyond MHD, through including the ambipolar diffusion term, the Hall term, and/or Biermann’s battery term, to the induction equation (Martínez-Sykora et al. 2012; Cheung and Cameron 2012; Khomenko et al. 2017, 2018). These inclusions retain the single fluid description. This greatly simplifies the treatment of the radiative term $Q_{{\rm rad}}$, because energy lost or gained by the gas modifies the internal energy of the gas as a single entity, instead of modifying the energy of the electrons and different species of atoms, ions, and molecules separately.

Efforts are underway to move beyond single-fluid radiation-MHD to a multi-fluid description, treating neutrals, ions, and/or electrons as separate fluids each with their own temperatures and velocities. Radiative transitions can modify the internal energy and change the ionisation state of atoms, and in a multi-fluid description these must be computed in detail (e.g., Khomenko et al. 2014). This review does not discuss radiation-hydrodynamics in multi-fluid models.

2.2 Energy density of radiation and matter

It is useful to compare the energy density and flux of the radiation field to the energy density of the gas in the solar atmosphere. To get an estimate we assume that the radiation field is isotropic and given by the Planck function at $T_{{\rm rad}}=5777$ K. Then the energy density is

$$\begin{aligned} E_{{\rm rad}}= \frac{4\sigma }{c} T_{{\rm rad}}^4 = 0.84\,{{\rm J}}\,{{\rm m}}^{-3}, \end{aligned}$$

(4)

with $\sigma$ the Stefan–Boltzmann constant. Assuming that the solar surface is a radiating blackbody at the same temperature gives the radiative flux at the surface as

$$\begin{aligned} F_{{\rm rad}}= \sigma T_{{\rm rad}}^4 = 6.3\times 10^7 \,{{\rm W}}\,{{\rm m}}^{-2}, \end{aligned}$$

(5)

Table 1 compares the radiative energy density to the internal energy density e of the solar gas for typical photospheric, chromospheric and coronal values assuming the solar gas has an ideal gas equation of state ($e=n k_{{\rm B}} T$). The radiative energy density in the chromosphere and corona is corrected for the non-isotropy of the radiation above the photosphere. The energy density of the radiation is much lower than the energy density of the gas in the photosphere and chromosphere, but in the corona they are about equal. The corona is however optically thin for most photospheric radiation and only little absorption or scattering of radiative energy by the gas occurs.

Table 1 Comparison of energy densities and radiative flux

2.3 Radiation pressure and force

The force exerted by the radiation pressure tensor is typically ignored under normal solar conditions because it is small compared to other forces. To illustrate this one can compare the radiation pressure of isotropic black body radiation to the gas pressure in the photosphere. This pressure is

$$\begin{aligned} \frac{4 \sigma }{3c} T^4 = 0.27\,{{\rm Pa}} \end{aligned}$$

(6)

at a photospheric temperature of $T=5700$ K, while the gas pressure in the photosphere is roughly $10^4$ Pa. Similarly, one can compute a rough estimate of the upward acceleration of photospheric material by radiation:

$$\begin{aligned} a_{{{\rm rad}}} \approx \frac{\kappa F}{c} = 3 \times 10^{-3}\,{{\rm m}}\,{{\rm s}}^{-2}, \end{aligned}$$

(7)

with $\kappa$ the Rosseland opacity per mass unit and F the frequency-integrated radiative flux, both taken in the photosphere. The radiative acceleration is much smaller than the downward-directed solar surface gravity acceleration $g=274$ m s$^{-2}$.

2.4 Energy exchange between radiation and matter

In absence of an absorbing medium, the monochromatic radiative flux divergence $\nabla \cdot {\mathbf{F}} _\nu$ is zero. If a medium (in the solar case a gas or plasma) is present, then the rate per volume with which the material gains energy from the radiation field is given by

$$\begin{aligned} Q_{{\rm rad}}= - \int _0^\infty \nabla \cdot {\mathbf{F}} _\nu \, {{\rm d}}\nu = -\nabla \cdot {\mathbf{F}} , \end{aligned}$$

(8)

with ${\mathbf{F}}$ the total radiative flux. The total radiative flux is an integral over frequency, and as far as the internal energy of the gas is concerned, the exact distribution of the flux over frequency is not important. Typically, radiation-MHD simulations of the solar atmosphere aim to reproduce the correct detailed behaviour of the gas only. The computation of the radiation field only has to reproduce the correct heating and cooling, but does not need to reproduce the correct spectral energy distribution. This allows for large simplifications in treating the radiation without sacrificing too much accuracy in the value of the total flux divergence.

2.5 Explicit expression of the radiative flux divergence

The monochromatic flux is defined in terms of the intensity $I_\nu$ as

$$\begin{aligned} {\mathbf{F}} _\nu = \oint {\mathbf{n}} I_\nu \,{{\rm d}}\varOmega , \end{aligned}$$

(9)

with ${\mathbf{n}}$ the unit vector pointing in the direction of $\varOmega$. Substitution of this equation into the integral over all directions of Eq. (14) yields an expression of the total radiative flux divergence in terms of emissivity or source function $S_\nu =j_\nu /(\kappa _\nu \rho )$, intensity and opacity:

$$\begin{aligned} \nabla \cdot {\mathbf{F}}= & {} \int _0^\infty \oint \left( j_\nu - \kappa _\nu \rho I_\nu \right) \, {{\rm d}}\varOmega \,{{\rm d}}\nu \nonumber \\= & {} \int _0^\infty \oint \kappa _\nu \rho \left( S_\nu - I_\nu \right) \, {{\rm d}}\varOmega \,{{\rm d}}\nu , \end{aligned}$$

(10)

with $\kappa _\nu$ the extinction coefficient per unit mass. If one assumes that both the emissivity and extinction coefficient do not depend on direction then Eq. (10) reduces to

$$\begin{aligned} \nabla \cdot {\mathbf{F}} = \int _0^\infty 4 \pi \kappa _\nu \rho \left( S_\nu - J_\nu \right) \, {{\rm d}}\nu , \end{aligned}$$

(11)

where

$$\begin{aligned} J_\nu = \frac{1}{4 \pi } \oint I_\nu \,{{\rm d}}\varOmega \end{aligned}$$

(12)

is the angle-averaged intensity.

This assumption is reasonable for bound-free transitions and other continuous processes. For bound-bound transitions it is not valid when flow velocities are of the same order or larger than the thermal speed $\sqrt{2kT/m+v_{{{\rm turb}}}^2}$, with m the mass of the atom or ion that is involved, and $v_{{\rm turb}}$ the microturbulent velocity. Nevertheless, the assumption is almost always taken to be valid because it allows large simplifications.

2.6 Light travel time and hydrodynamical timescales

Typical bulk flow velocities in the solar atmosphere range from 1 km s$^{-1}$ in convective upflows to $>300$ km s$^{-1}$ in chromospheric evaporation following solar flares (Graham and Cauzzi 2015). Alfvén velocities up to $2.2\times 10^3$ km s$^{-1}$ have been measured in the solar corona by studying properties of coronal loop oscillations (e.g., Aschwanden et al. 1999; Pascoe et al. 2018). These values are much lower than the speed of light.

The typical hydrodynamic timescale in the photosphere is $\sim 5$ min, it is $\sim 1$ min in the chromosphere and a few minutes in the corona. The light crossing time is of the order of a second even for the largest coronal structures.

Because of the low flow speeds and long hydrodynamic timescales compared to light crossing times, it is typically assumed that light travel times can be ignored in solar radiation hydrodynamics. That means that the computation of the radiative flux divergence at a given time t only depends on the state of the atmosphere at time t. The time-dependent transfer equation for the intensity in direction ${\mathbf{n}}$ then simplifies from

$$\begin{aligned} \frac{1}{c} \frac{\partial I_\nu }{\partial t} + {\mathbf{n}} \cdot \nabla I_\nu = j_\nu - \alpha _\nu I_\nu \end{aligned}$$

(13)

to

$$\begin{aligned} {\mathbf{n}} \cdot \nabla I_\nu = j_\nu - \alpha _\nu I_\nu , \end{aligned}$$

(14)

with $j_\nu$ the emissivity and $\alpha _\nu = \kappa _\nu \rho$ the extinction coefficient per volume. In other words, it is assumed that the radiation field adapts instantaneously to changes in the state of the solar gas.

2.7 Diffusion approximation

At very large optical depth the following approximation holds for $J_\nu$:

$$\begin{aligned} J_\nu \approx S_\nu +\frac{1}{3} \frac{{{\rm d}}^2 S_\nu }{{{\rm d}}\tau _\nu ^2}, \end{aligned}$$

(15)

where the second derivative should be taken along the direction of the gradient of $S_\nu$.¹ For the flux the following approximation holds:

$$\begin{aligned} {\mathbf{F}} _\nu \approx - \frac{4 \pi }{3 \kappa _\nu \rho } \nabla S_\nu . \end{aligned}$$

(16)

At depths where the optical depth at all frequencies is much larger than unity, one can derive the diffusion approximation for the total radiative flux from Eq. (16) if one also assumes that the source function is equal to the Planck function $B_\nu$:

$$\begin{aligned} {\mathbf{F}} \approx - \frac{16}{3} \frac{\sigma T^3}{\kappa _{{{\rm R}}}\rho } \nabla T, \end{aligned}$$

(17)

with the Rosseland opacity defined as

$$\begin{aligned} \frac{1}{\kappa _{{{\rm R}}}} = \left( \int _0^\infty \frac{1}{\kappa _\nu } \frac{{{\rm d}}B_\nu }{{{\rm d}}T} \, {{\rm d}}\nu \right) {\Big /}\left( \int _0^\infty \frac{{{\rm d}}B_\nu }{{{\rm d}}T} \, {{\rm d}}\nu \right) . \end{aligned}$$

(18)

Methods to compute the total flux divergence that rely on the solution of the radiative transfer equation must converge to the diffusion approximation at large depth. Methods that require a numerical approximation of the source function must use one of at least second order in order to recover Eq. (15).

3 Radiative transfer in the photosphere: multi-group radiative transfer

Multi-group radiative transfer was first introduced in Nordlund (1982). It is a method that approximates the frequency integral in Eq. (11) by a sum over a limited number N so called radiation groups (or radiation bins):

$$\begin{aligned} \int _0^\infty 4 \pi \kappa \rho \left( S_\nu - J_\nu \right) \, {{\rm d}}\nu \approx \sum _{i=1}^N 4 \pi \kappa _i \rho \left( S_i - J_i \right) , \end{aligned}$$

(19)

with $\kappa _i$, $S_i$, and $J_i$ the opacity, source function and radiation field in each group, as defined below.

The rationale of the method is the realisation that the $\varLambda$-operator that produces the angle-averaged radiation field from the source function depends on the opacity only, and is a linear operator. While it is customary to write $\varLambda _\nu [\ldots ]$ because the opacity changes with frequency, a better notation would be $\varLambda _{\kappa }[\ldots ]$. For two frequencies $\nu _1$ and $\nu _2$ with identical opacities everywhere in the atmosphere $\kappa _1 = \kappa _2 = \kappa$ one can thus write

$$\begin{aligned} J_{\nu _1} + J_{\nu _2} = \varLambda _{\nu _1} [S_{\nu _1}] + \varLambda _{\nu _2} [S_{\nu _2}] = \varLambda _\kappa [S_{\nu _1}+S_{\nu _2}] \end{aligned}$$

(20)

In order to derive the approximation in Eq. (19) one first approximates the frequency integral with a sum over discrete frequencies $\nu _j$ with summation weights $w_j$. These frequencies are then grouped into N bins, each labeled i, that have similar opacities. Finally, Eq. (20) is used to define a bin-integrated source function:

$$\begin{aligned} \int _0^\infty \kappa \left( S_\nu - J_\nu \right) \, {{\rm d}}\nu \simeq \sum _j w_j \kappa _j \left( S_j - J_j \right) \end{aligned}$$

(21)

$$\begin{aligned} = \sum _{i=1}^N\sum _{j(i)} w_j \kappa _j \left( S_j -J_j \right) = \sum _{i=1}^N \sum _{j(i)} w_j \kappa _j \left( S_j - \varLambda _j[S_j] \right) \end{aligned}$$

(22)

$$\begin{aligned} \approx \sum _{i=1}^N \kappa _i \left( \sum _{j(i)} w_j S_j - \varLambda _i\left[\sum _{j(i)} w_j S_j\right] \right) \equiv \sum _{i=1}^N \kappa _i \left( S_i - \varLambda _i[S_i] \right) \end{aligned}$$

(23)

$$\begin{aligned}&\equiv \sum _{i=1}^N \kappa _i \left( S_i -J_i \right) . \end{aligned}$$

(24)

Here j(i) is the set of all frequencies with similar opacities that are placed in bin i. The bin-integrated source function is constructed from the frequency-dependent source function using

$$\begin{aligned} S_i = \sum _{j(i)} w_j S_j \end{aligned}$$

(25)

What is left is now is to define how to group frequency points in the various bins, how to define the representative opacity for each bin $\kappa _i$, and how to define the frequency-dependent source function $S_j$. The opacity and source function are in general functions of at least frequency, temperature and density. Additional dependencies are on the velocity field in the atmosphere, the ray direction, possible location-dependent abundances, and in the case of non-LTE radiative transfer, on the radiation field.

Section 3.1 discusses how to sort frequencies into opacity groups. In Sect. 3.2 it is assumed that the source function is given by the Planck function. This assumption is relaxed to allow for coherent scattering in Sect. 3.3.

3.1 Sorting frequencies into groups

Grouping frequencies and defining a bin opacity are necessarily somewhat crude. The density and thus opacities drop roughly exponentially with height, and opacities at a given point in the atmosphere vary over many orders of magnitude with frequency. The aim of the multi-group opacities is to approximate the full radiative transfer with sensitivity spanning from the low-opacity continua that form in the photosphere to strong lines that form in the mid and ideally even the upper chromosphere.

The de-facto standard method for grouping opacities is called $\tau$-sorting: first a 1D reference atmosphere is chosen, for which the opacities and vertical optical depth scales at all frequency points j are computed. The atmosphere is then divided into various height bins. Frequencies at which optical depth unity is reached in height bin i are assigned to radiation bin i. This height sorting can in principle be done directly on a geometrical height scale, but most radiation-MHD codes follow Nordlund (1982). They set the borders between the height bins in terms of the Rosseland optical depth scale $\tau _{{{\rm R}}}$: let the borders be defined at a number of Rosseland optical depths $\tau _{{{\rm R}}}^k$. A frequency $\nu _j$ then belongs to bin i if

$$\begin{aligned} \tau _{{{\rm R}}}^{k-1}< \tau _{{{\rm R}}}(\tau _\nu =1) \le \tau _{{{\rm R}}}^k. \end{aligned}$$

(26)

A common choice is to put the border values $\tau _{{{\rm R}}}^k$ equidistantly spaced in $\log {\tau _{{\rm R}}}$. The method is illustrated in Fig. 1.

Fig. 1

Image reproduced with permission from Ludwig (1992)

Illustration of the principle of $\tau$-sorting. The horizontal axis is frequency, while the vertical axis shows the Rosseland optical depth in the 1D reference atmosphere. The solid line is the Rosseland optical depth where the monochromatic optical depth unity is reached. The frequencies are divided into four bins using the three border values $\tau _{{{\rm R}}}^1$–$\tau _{{{\rm R}}}^3$.

Opacities are generally well-approximated by their LTE values in the photosphere. LTE opacities in a static atmosphere are functions of temperature, density and elemental composition only, but even then it requires a large effort to accurately compute them. For historical reasons this is commonly done using an intermediate step called opacity distribution functions (ODFs). ODFs were developed originally to speed up computations used in modeling of 1D LTE radiative equilibrium stellar atmospheres (e.g., Gustafsson et al. 1975; Kurucz 1979).

The method used to compute an appropriate mean opacity in a bin i from the opacities $\kappa _j$ depend on whether one assumes an LTE or non-LTE source function and are described in Sects. 3.2 and 3.3. An illustrative solution assuming an LTE source function is given in Fig. 2.

Fig. 2

Image reproduced with permission from Vögler et al. (2004), copyright by ESO

Illustration of the concept of group mean opacity. Left-hand panel: monochromatic extinction coefficients based on opacity distribution functions for $T=4470\,{{\rm K}}$ and $P = 1.3\times 10^4\,{{\rm Pa}}$, corresponding roughly to the upper photosphere. Right-hand panel: The opacity that is “assigned” to each wavelength using a 5-group scheme. The large continuous variation in the left-hand panel is replaced by only five discrete opacities. Note that the high monochromatic opacities in the UV are replaced by a much lower group opacity.

3.2 Multi group radiative transfer with LTE source function

The source function in the solar atmosphere is in general not equal to the Planck function because at some height in the atmosphere densities become too low to set up Saha–Boltzmann populations through collisions. However, detailed non-LTE computations in 1D models show that in the deep photosphere (roughly defined here as $-100\,{{\rm km}}<z< 200\,{{\rm km}}$, with the $z=0$ defined as the location where $\tau _{500\,{{\rm nm}}}=1$), the source function is almost exactly equal to the Planck function (see Fig. 36 of Vernazza et al. 1981), and the opacities are close to their LTE values. At larger heights this is no longer the case for lines and continua of neutral atoms with a low ionisation potential because they tend to be over-ionized by the radiation from below. Nevertheless, one can expect that assuming LTE for both the source function and the opacity is a good approximation for computing the flux divergence in the photosphere. The source function in group i is then given by

$$\begin{aligned} S_i = \sum _{j(i)} w_j B_j. \end{aligned}$$

(27)

This expression only depends on temperature and can thus easily be precomputed and stored in a 1D lookup table.

Once the opacities $\kappa _j$ are grouped into the N bins, one still needs to compute an appropriate bin opacity $\kappa _i$. Choosing the numerical equivalent of the Rosseland opacity in each bin ensures that the diffusion approximation is recovered deep in the atmosphere:

$$\begin{aligned} \frac{1}{\kappa _i^{{\rm R}}} = \left( \sum _{j(i)} w_j \frac{1}{\kappa _j} \frac{{{\rm d}}B_j}{{{\rm d}}T} \right) {\Big /} \left( \sum _{j(i)} w_j \frac{{{\rm d}}B_j}{{{\rm d}}T} \right) . \end{aligned}$$

(28)

This choice of bin opacity is however not appropriate at low optical depths, where photons are mainly in the free streaming regime. Following approximations valid at small optical depth put forward in Mihalas (1970), and further developed in Ludwig (1992), it turns out that the Planck-mean opacity $\kappa _B$ is a good choice for the outermost layers of the atmosphere:

$$\begin{aligned} \kappa _{i}^B = \frac{\sum _{j(i)} w_j \kappa _j B_j}{ \sum _{j(i)} w_j B_j}. \end{aligned}$$

(29)

Somewhere in the atmosphere one should make a transition from one opacity definition to the other. This can be achieved through defining the group opacity as

$$\begin{aligned} \kappa _i = W \kappa _i^B + (1-W) \kappa _i^{{\rm R}}, \end{aligned}$$

(30)

$$\begin{aligned} W = { {\rm e}}^{-\tau _i/\tau _0}. \end{aligned}$$

(31)

where $\tau _0$ is an adjustable parameter of order unity and $\tau _i$ the vertical optical depth. Computing $\tau _i$ in the MHD simulation is somewhat computationally expensive and instead it is estimated, for example through using the relation between mass density and optical depth in the 1D reference atmosphere used for the sorting of frequencies into bins (Ludwig 1992; Vögler et al. 2004).

Another option is to base the transition on the local mean free path:

$$\begin{aligned} W = {{\rm e}}^{-l_i/(\kappa _i \rho )}, \end{aligned}$$

(32)

where $l_i$ is a typical length scale over which the bin-integrated source function changes (Skartlien 2000). The advantage of this method is that it does not depend on the properties of a 1D reference model.

For a fixed elemental composition, $\kappa _i$ depends only on a combination of any two thermodynamic parameters (for example, e and $\rho$). Like the group-integrated source function, they are commonly precomputed and put in a 2D lookup table.

Extensive discussions of multi-group radiative transfer assuming an LTE source function can be found in the PhD-thesis of Ludwig (1992, in German) and in Vögler et al. (2004).

3.3 Multi group radiative transfer with non-LTE source function

The assumption of LTE for the source function is accurate in the photosphere, but breaks down in the chromosphere and transition region. Here the energy exchange is dominated by the resonance lines of H i, Ca ii, Mg ii, and He ii (see, e.g., Fig. 49 of Vernazza et al. 1981). These lines can have photon destruction probabilities

$$\begin{aligned} \epsilon = \frac{C_{ul}}{A_{ul} + C_{ul} + B_\nu B_{ul}} < 10^{-4}, \end{aligned}$$

(33)

with $C_{ul}$ the downward collisional rate coefficient and $A_{ul}$ and $B_{ul}$ Einstein coefficients for spontaneous and induced deexcitation. Scattering should therefore be taken into account.

Scattering has a strong damping effect on the amplitude of $\nabla \cdot {\mathbf{F}} _\nu /\rho$. In LTE it is given by

$$\begin{aligned} \frac{\nabla \cdot {\mathbf{F}} _\nu ^{{\rm LTE}}}{\rho } = 4 \pi \kappa _\nu \left( B_\nu - J_\nu \right) . \end{aligned}$$

(34)

If one assumes the presence of a coherently scattering line in a two-level atom at frequency $\nu$, then the source function becomes

$$\begin{aligned} S_\nu = (1-\epsilon ) j_\nu + \epsilon B_\nu . \end{aligned}$$

(35)

In that case the flux divergence per mass unit is

$$\begin{aligned} \frac{\nabla \cdot {\mathbf{F}} _\nu ^{{\rm NLTE}}}{\rho } = 4 \pi \epsilon \kappa _\nu \left( B_\nu - J_\nu \right) , \end{aligned}$$

(36)

so that for a given difference between $B_\nu$ and $J_\nu$, the amplitude of the radiative energy exchange in non-LTE can be orders of magnitude smaller than in LTE.

Skartlien (2000) extended the method presented in Sect. 3.2 to include coherent two-level scattering, as an approximation to the much more complicated full non-LTE radiative transfer problem. He assumed that the monochromatic extinction coefficient can still be computed in LTE. This assumption is rather accurate for the resonance lines of H i, Ca ii, Mg ii, whose lower levels are the ground state of a dominant ionisation state. In addition he assumed that the scattering coefficient in a spectral line is given using the approximation by van Regemorter (1962), so that it is independent of the actual line, and only depends on frequency, temperature and electron density. Starting from the monochromatic two level source function $S_\nu = (1-\epsilon _\nu ) J_\nu + \epsilon _\nu B_\nu$ and following a reasoning similar to the one given in Sect. 3.2, he arrives at an expression for the group source function:

$$\begin{aligned} S_i = \epsilon _i J_i + t_i. \end{aligned}$$

(37)

Here $\epsilon _i$ represents a group-averaged scattering probability, and $t_i$ the group-integrated thermal production of photons. He also derives an expression for the group extinction coefficient $\kappa _i$. Similar to the LTE case, $\epsilon _i$, $t_i$, and $\kappa _i$ have different expressions for the diffusion regime and in the free streaming regime in the outer atmosphere. An important difference is that the streaming quantities now depend on the monochromatic mean intensity $J_\nu$ in a 1D reference atmosphere. This means that the quantities must be recalculated for each different target of simulations (e.g., sunspots or quiet Sun). Simulations containing a variety of structures might suffer from inaccuracies because a 1D atmosphere cannot be representative of all atmospheric structures.

Another difference compared to the method employing an LTE source function is that the computation of $Q_{{\rm rad}}$ now requires solution of Eq. (37) together with the transfer equation because $J_i= \varLambda [S_i]$. This is typically done through accelerated $\varLambda$-iteration. Because of the multidimensional geometry, a local approximate $\varLambda$-operator (which is equivalent to Jacobi-iteration, see Olson et al. 1986) is efficient and easily coded, such as in the Oslo Stagger Code (Skartlien 2000). Gauss–Seidel iteration (Trujillo Bueno and Fabiani Bendicho 1995) offers superior convergence speed and has been implemented in the Bifrost code (Hayek et al. 2010; Gudiksen et al. 2011).

3.4 Solving the transfer equation

Most modern radiation-MHD codes are parallelized to make use of large supercomputers with distributed memory. Typically, the simulated domain is represented on a 3D Cartesian grid. The domain is split into smaller subdomains and each CPU handles the required computations on its own subdomain, communicating information to other subdomains as needed. This architecture scales well for the MHD equations: they are local and require communication with neighbouring subdomains only.

Radiation is however intrinsically non-local. An emitted photon can traverse many subdomains before undergoing another interaction with the solar gas, so communication between subdomains is not necessarily local. This problem is shared between radiation-MHD codes that aim to compute a reasonable approximation of the flux divergence, and non-LTE radiative transfer codes such as Multi3D (Leenaarts and Carlsson 2009) and Porta (Štěpán and Trujillo Bueno 2013), which aim to accurately compute the emergent spectrum from a given model atmosphere. Consequently, there is a large amount of literature addressing efficient solutions of the transfer equation in multidimensional geometries and/or decomposed domains (e.g., Kunasz and Auer 1988; Auer 2003; Vögler et al. 2005; Heinemann et al. 2006; Hayek et al. 2010; Ibgui et al. 2013; Štěpán and Trujillo Bueno 2013). I will only briefly touch upon some aspects.

Solving for the flux ${\mathbf{F}}$ or for the angle-averaged radiation field J requires computation of the intensity for a number of different directions at all points on the numerical grid. In practice there are two methods that are used in radiation-MHD codes: short characteristics (SCs) and long characteristics (LCs). Both are illustrated in Fig. 3.

Fig. 3

Illustration of short and long characteristics used for solving the transfer equation in 2D decomposed domains. Grey dots indicate grid points where the intensity should be computed, with grey lines connecting those grid points. Blue lines indicate subdomain boundaries, where it is assumed that the horizontal boundaries are periodic. Red arrows indicate two examples of long characteristics, while the orange arrows indicate short characteristics. Note that the long characteristics cross multiple subdomain boundaries and wrap around the periodic horizontal boundary each subdomain contains a piece of the long characteristic of variable length

Short characteristics solve the transfer equation along short line segments (orange in Fig. 3), starting at a grid cell boundary and ending at a grid point where the intensity is desired. Along the line one typically computes a numerical approximation of the formal solution:

$$\begin{aligned} I(\tau ) = I(\tau =0) \, {{\rm e}}^{-\tau } + \int _0^\tau S(t) \, {{\rm e}}^{t-\tau } {{\rm d}}t, \end{aligned}$$

(38)

where the optical thickness scale has its zero point at the start of the SC and $\tau$ is the optical thickness along the entire SC. Intensities $I(\tau )$ are computed at the grid points (grey circles in Fig. 3). Computation of $I(\tau =0)$, which is needed at the start of the orange arrow, thus requires interpolation from the grid points. SC methods are therefore somewhat diffusive and coherent beams of photons disperse. High-order interpolation schemes can alleviate, but not eliminate, this diffusion. In practice this diffusion is typically not a problem, given that the photosphere is emitting photons everywhere and both the source function and the resulting radiation field are rather smooth.

The transfer equation is solved along all SCs in a sequential order, starting from a known boundary condition (the diffusion approximation at the bottom of the atmosphere for rays going up, and typically zero for SCs going down from the top of the atmosphere). The method was introduced for 2D cartesian geometry by Kunasz and Auer (1988). An in-depth description of the method in 3D Cartesian geometry is given in Ibgui et al. (2013). The ordered fashion in which SCs must be computed leads to complications with spatial domain decomposition. An example method of how to achieve reasonable parallelism despite this ordering can be found in Štěpán and Trujillo Bueno (2013). Short characteristics can be easily computed along any angle. Typical ray quadratures [i.e., the set of angles chosen to numerically compute Eq. (9) or (12)] that are in use are the angle sets from Carlson (1963), or equidistant in azimuth and using a Gaussian quadrature in inclination.

Bruls et al. (1999) present a method to compute SCs on unstructured grids. It is not currently in use in the common radiation-MHD codes that use fixed Cartesian grids, but might be of great use for codes that use unstructured or adaptive grids.

The long characteristic method traces rays from the lower boundary of the domain to the upper boundary (red line segments in Fig. 3). An LC does generally intersect only a few grid points. Interpolation of the source function and the opacity from the grid points to the LC and interpolation of the intensity along the LC back to the grid points is thus necessary. The transfer equation along an LC can be solved using the formal solution [Eq. (38)], or by solving the second-order form of the transfer equation (Feautrier 1964).

Long characteristics allow photons to travel in a straight line and are thus not diffusive. Efficient parallel algorithms exist for solving along LCs in decomposed domains (Heinemann et al. 2006). However, this parallel algorithm is only easily implemented when LCs cross cell boundaries exactly through grid points, which limits the application to grids with fixed spacing and a maximum of 26 directions: both directions along three axes, six face diagonals, and 4 space diagonals in a cubic grid (e.g., Popovas et al. 2019). Arbitrary ray quadratures can be implemented at the expense of code simplicity.

Both the SC and LC method can handle exactly horizontal rays, but such rays require implicit solution methods in case of periodic boundary conditions in the horizontal direction, which is the default for solar atmosphere radiation-MHD simulations. This adds additional coding complications in the usual case that parallelism is implemented through spatial domain composition. Horizonal ray directions are therefore usually avoided.

3.5 Computation of the heating rate from the intensity and source function

Analytically, the two expressions for the heating rate in a bin $Q_i=-\nabla \cdot {\mathbf{F}} _i = 4\pi \kappa _i \rho _i (J_i-S_i)$ are equivalent. In case of actual numerical computation this is no longer the case.

In the deep atmosphere the diffusion approximation holds. Equation (15) shows that while $S_i$ and $J_i$ increase with depth because of the increase in temperature, their difference becomes smaller, and at some point roundoff errors become noticeable. These errors are then amplified by the exponential increase of the density with depth. Using the source function and radiation field to compute $Q_i$ is thus unstable in the deep layers. Instead, computing the flux from the intensity using the numerical equivalent of Eq. (9) and then taking the divergence is stable in the deep layers. Because the radiative flux is small compared to the energy density of the gas, an error in computation of its divergence will not lead to large errors in the internal energy.

In the upper layers the situation is reversed: the radiative flux is large compared to the energy density of the gas [see Table 1 and Eq. (4)], and a small error in computation of the flux divergence from the intensity will lead to a large error in $Q_i$ and e. Using $4\pi \kappa _i \rho _i (J_i-S_i)$ is stable however, because of the generally large split between $S_i$ and $J_i$.

Following the suggestion by Bruls et al. (1999), most 3D codes that use short characteristics compute $Q_i$ using the flux divergence at large optical depth, and switch to using the source function and radiation field around an optical depth of 0.1.

An alternative to the above switching scheme is to rewrite the transfer equation in terms of the quantity $K^I=S-I$ (dropping bin indices i here for brevity). The quantity $K^I$ is proportional to the cooling rate in a specific ray direction. The transfer equation then transforms into

$$\begin{aligned} \frac{{{\rm d}}K^I}{{{\rm d}}\tau } = \frac{{{\rm d}}S}{{{\rm d}}\tau }- K^I, \end{aligned}$$

(39)

which in its integral form is given by

$$\begin{aligned} K^I(\tau ) = K^I(\tau _0) {{\rm e}}^{\tau -\tau _0} + \int _{\tau _0}^\tau {{\rm e}}^{\tau '-\tau } \frac{{{\rm d}}S}{{{\rm d}}\tau '}\, {{\rm d}}\tau '. \end{aligned}$$

(40)

The total heating rate is then computed from

$$\begin{aligned} Q= \oint K^I \,{{\rm d}}\varOmega \end{aligned}$$

(41)

This equation does not suffer from the numerical precision problems caused by the cancellation of the subtraction S and I. Equation (40) has the same form as the formal solution of the normal transfer equation and can be solved efficiently using a variety of methods. If the source function is known (such as when assuming LTE) then solving straight for K is possible without having to solve for I first. Heinemann et al. (2006) describe an elegant method for solving for K in decomposed domains using long characteristics and a direct solution of the transfer equation. In case of a non-LTE source function, such as in Sect. 3.3, then computing S requires solving the transfer equation to obtain I and J. Computing Q from $K^I$ afterwards then offers little benefit over computing Q straight from $\nabla \cdot {\mathbf{F}}$ or $S-J$.

Nordlund (1982) implemented a similar method as Heinemann et al. (2006), but based on the second-order form of the transfer equation:

$$\begin{aligned} \frac{{{\rm d}}^2 P}{{{\rm d}}\tau ^2} = P-S. \end{aligned}$$

(42)

where

$$\begin{aligned} P = \frac{1}{2} \left( I(\varOmega ) + I(-\varOmega ) \right) , \end{aligned}$$

(43)

the average of an ingoing and an outgoing ray (e.g., Hubeny and Mihalas 2014, p. 387). Defining $K_P=P-S$ one arrives at

$$\begin{aligned} \frac{{{\rm d}}^2 K_P}{{{\rm d}}\tau ^2} = K_P - \frac{{{\rm d}}^2 S}{{{\rm d}}\tau ^2}. \end{aligned}$$

(44)

This equation has the same form as the second order transfer equation and can be solved efficiently along a long characteristic using the method of Feautrier (1964).

3.6 Summary and examples of photospheric radiative transfer

Figure 4 shows a summary table of the three major binary choices in multi-group radiative transfer from computing radiative losses and gains in the photosphere: LTE or non-LTE source function, short characteristics or long characteristics, and solving for K or solving for I in order to compute the flux divergence. Each of the resulting 8 options is numbered.

Fig. 4

Summary table of the three major binary choices in multi-group radiative transfer from computing radiative losses and gains in the photosphere

The simulation by Nordlund (1982) is of Type 5. MURaM (Vögler 2004; Rempel 2017), RAMENS (Iijima and Yokoyama 2015; Iijima 2016), and MANCHA3D (Khomenko et al. 2017, 2018) are Type 2 codes. Stellarbox (Wray et al. 2015) is Type 6. CO5BOLD (Freytag et al. 2012) has both Type 2 and Type 6 options, but for stellar atmosphere simulations only Type 2 is supported.

The Oslo version of the Stagger code (Type 8) and the Bifrost code (Gudiksen et al. 2011) (Type 4) are as of October 2019 the only codes using a non-LTE source function.

Figures 5 and 6 demonstrate the result of a 4-group computation in a radiation-MHD simulation. The details of this particular simulation can be found in Carlsson et al. (2016).

Fig. 5

Example of the vertically emergent intensity computed with a 4-group non-LTE scheme. The flux divergence along the horizontal white line is shown in Fig. 6. The model was computed with the Bifrost code (Gudiksen et al. 2011)

Fig. 6

Atmospheric structure and flux divergence per mass unit $Q_i/\rho =-\nabla \cdot {\mathbf{F}} _i/\rho = 4 \pi \kappa _i (J_i-S_i)$ in a vertical slice along the white line in Fig. 5. The top row shows the temperature and density, the two bottom rows show the flux divergence per mass unit in the four opacity groups. Brown is cooling, blue is heating, and the brightness scale for the flux divergence panels is clipped at 20% of the maximum of the absolute value to enhance contrast. The maximum and minimum values are indicated in each panel. The black curves indicate the $\tau _i=1$ height in each bin. Note that the maxima of the heating and cooling per unit mass do not coincide with the $\tau =1$ height. The flux divergence in this simulation is artificially set to zero at the height where the entire atmosphere above has a maximum optical thickness of $10^{-5}$. The total radiative losses above this height are negligible, but the local losses and gains per unit mass are not

Figure 5 shows the vertically emergent intensity in each radiation group. Bins 1 and 2 grouped areas of the spectrum with generally low opacities. The corresponding images are indeed reminiscent of observed optical continuum images of the photosphere. Bins 3 and 4 contain opacities of stronger lines, and resemble images taken in the wings of the Ca ii H&K lines (e.g., Rutten et al. 2004). The images are dominated by reversed granulation; the small bright structures are caused by magnetic field concentrations.

Figure 6 shows the heating rate per mass in each radiation bin. The prominent funnel shape in the mass density panel is caused by the presence of a strong magnetic field concentration that fans out with height. All bins show strong cooling at the top of the granules (Red color just below $z=0.0$ Mm), and bins 2–4 show heating just above the granules. There is strong cooling per mass unit in the chromosphere above $z=1$ Mm. While the optical thickness of the chromosphere above this height is small because of the low density, the heating rate per mass is independent of the mass density, and depends on the value of $\kappa _i$ and the size of $(S_i-J_i)$ only.

Figure 7 demonstrates the accuracy of approximating the spectrum by only a few frequency groups in the non-LTE scheme of Skartlien (2000) on the horizontally averaged values of $Q_i$. It does however not test the assumptions of coherent scattering, LTE equation of state and LTE opacity. The absolute error in $Q = \varSigma Q_i$ is of the order of a few percent, while the error in individual groups i can be as large as 50% (group 3 at $z=0.3$ Mm). Note that Q is dominated by group 1. The absolute value of Q is decreasing with height because of its dependence on the mass density.

Fig. 7

Image reproduced with permission from Skartlien (2000), copyright by AAS

Average radiation heating and cooling per volume (Q) as a function of height in a 3D radiation-hydrodynamics simulations using a 4-group non-LTE radiative transfer scheme. Upper panels: Horizontally averaged radiative heating per volume unit shown in a bilogarithmic plot. The units are arbitrary. Thick lines are the exact group solutions, while the thin lines are the approximate group solutions. The vertical line segments in the uppermost panel show the heights where the average optical depth per group is unity, as measured along the vertical line. Lower panel: Horizontally averaged amplitudes of radiative heating (average absolute value) relative to the amplitude of the exact total heating. Thick curves are the exact amplitude ratios, while the thin curves are ratios from the approximate group solutions.

In Fig. 8 the difference between assuming an LTE or non-LTE source function is demonstrated. The expression for the flux divergence in the multigroup scheme is the same as for the monochromatic case:

$$\begin{aligned} \frac{\nabla \cdot {\mathbf{F}} _i}{\rho } =4 \pi \epsilon _i \kappa _i \left( S_i - J_i \right) , \end{aligned}$$

(45)

where $\epsilon _i = 1$ in LTE, and $\epsilon _i \le 1$ in non-LTE. In the latter case $\epsilon _i$ can be as low as $10^{-4}$. A strict comparison between the scattering and non-scattering cases is not possible, because of the different definition of the group-mean opacities, source functions and thermal emission terms. Nevertheless, the main result from this figure is that the LTE scheme vastly overestimates the radiative cooling in the chromosphere, because $\epsilon _i = 1$ in LTE.

Fig. 8

Image reproduced with permission from Skartlien (2000), copyright by AAS

Horizontal averages of radiative heating per mass unit ($Q/\rho$ for a 4-group scheme assuming LTE and non-LTE. The units are arbitrary, and the scaling is the same for all figures. The amplitude of the flux divergence for the LTE method (labeled old) is much larger above 0.3 Mm than for the non-LTE method (labeleld new). This effect is caused mainly by the larger differences between the source function and the mean intensity. Around 0.0 Mm, the new method produces larger amplitudes in the sense that hot regions are cooled more and cooler regions are heated more. Lower panel: Sum of all groups. The main difference above 1.0 Mm comes from the contribution in group 4.

The amplitude of the cooling and heating in groups 2–4 between $z=0$ Mm and $z=0.3$ Mm is however much larger in non-LTE than in LTE. Skartlien (2000) speculates that this is caused by the smoothing effect that scattering has on $J_i$, but a thorough investigation of this effect has never been done.

4 Radiative losses in the transition region and corona

The corona is optically thin at all wavelengths except for the radio regime. Its radiative losses are dominated by a myriad of EUV lines from highly ionised stages of many different elements (e.g., Curdt et al. 2004; Woods et al. 2012).

In the transition region, which I loosely define here as that part of the atmosphere where hydrogen is ionised but the temperature is below 100 kK, the emission is dominated by lines of lower ionisation stages (typically twice to four times ionised). For most lines, and for most regions on the sun, the TR is optically thin. In solar flares this is not necessarily the case: Kerr et al. (2019) showed that the TR can have appreciable optical thickness in the Si iv 140 nm lines.

Non-equilibrium ionisation effects play a role in the transition region (Olluri et al. 2013; Golding et al. 2017) and corona (Hansteen 1993; Bradshaw et al. 2004; Dzifčáková et al. 2016). Figure 1 of Hansteen (1993) shows an increase in radiative losses at a given density and temperature in the transition region by a factor two in a 1D hydrodynamic situation when non-equilibrium ionisation is used instead of instantaneous ionisation equilibrium.

A fully general treatment of TR and coronal radiative losses in a radiation-MHD simulation would thus involve solution of the full 3D non-equilibrium-non-LTE radiative transfer problem for most ionisation stages for a wide range of elements. This is currently impossible for 3D simulations because of limits on computation speed.

Neglecting radiative non-local transfer effects through assuming optically thin radiative transfer (i.e., assuming $J_\nu =0$ in the rate equations) alleviates already much of the problems without sacrificing much accuracy in most circumstances. If one furthermore excludes hydrogen and helium, the justification being that these elements are fully ionised at high temperatures and do not contribute to line cooling, then changes in ionisation of the elements do not influence the pressure and temperature.

The problem reduces then to solving the rate equations

$$\begin{aligned} \frac{\partial n_i}{\partial t} + \nabla \cdot (n_i{\mathbf{v}} ) = \sum _{j,j \ne i}^{n_{{{\rm l}}}} n_j P_{ji} - n_i \sum _{j,j \ne i}^{n_{{{\rm l}}}} P_{ij}, \end{aligned}$$

(46)

where i sums over all energy levels and ionisation stages of each element under consideration. The rate coefficients $P_{ij}$ contain collisional (de-)excitation and collisional ionisation recombination terms and spontaneous radiative deexcitation and recombination terms, but no radiative terms that involve absorption of an existing photon. Ignoring absorption therefore only allows for cooling. The radiative loss rate per volume in a bound-bound transition between a lower level i and upper level j is then given by

$$\begin{aligned} Q_{ij} = h \nu _{ij} A_{ji} n_j, \end{aligned}$$

(47)

and a similar expression can be written for bound-free transitions. The total cooling rate per volume can then be computed by summing the contributions of all transitions of all elements and including electron-ion free-free radiation. A more detailed description of this method in a 1D radiation-hydrodynamics code is given in Hansteen (1993). To my knowledge this method has not been implemented in 3D codes.

The option of computing non-equilibrium ionisation was added to the Bifrost code by Olluri et al. (2013), but they did not implement the resulting radiative cooling.

Instead, the default method to compute radiative losses in the corona is to assume statistical equilibrium [i.e., the left-hand-side of Eq. (54)] is assumed to be zero. Together with the assumption of no photon absorption these two assumptions together are often called the “coronal approximation”. The cooling a spectral line of element X in ionisation stage m with upper level j and lower level i can then be written as:

$$\begin{aligned} Q_{ij}= h \nu _{ij} \frac{A_{ji}}{n_{{\rm e}}} \, \frac{n_{j,m}}{n_m} \, \frac{n_m}{n_X} \, \frac{n_X}{n_{{\rm H}}} \, n_{{\rm e}}n_{{\rm H}}, \end{aligned}$$

(48)

$$\begin{aligned} \equiv G(T,n_{{\rm e}}) \, n_{{\rm e}}n_{{\rm H}}. \end{aligned}$$

(49)

Here $n_{j,m}/n_m$ is the fraction of all ions in ionisation stage m in level j, $n_m/n_X$ is the fraction of all atoms of species X in ionisation stage m, and $n_X/n_{{\rm H}}$ is the abundance of element X relative to hydrogen. The function $G(T,n_{{\rm e}})$ is only weakly dependent on the electron density: the upper level population is dominated by collisional excitation from the ground state, so that $n_{{\rm e}}n_m \sim A_{ji} n_{j,m}$, and the rate coefficients setting up the ionisation balance are almost linear in the electron density. A residual electron density dependence remains in $G(T,n_{{\rm e}})$ through collisional deexcitation and density-dependent dielectronic recombination (Summers 1972, 1974).

A total coronal radiative loss function $\varLambda (T,n_{{\rm e}})$ can thus be constructed through summing up $Q_{ij}$ for all levels and ionisation stages of all relevant elements and adding the contribution from continuum processes. The electron density dependence of $\varLambda (T,n_{{\rm e}})$ is weak: Landi and Landini (1999) performed a sensitivity study for electron densities in the range $10^{14}$–$10^{20}\,{{\rm m}}^{-3}$, and found a difference in the value of $\varLambda$ of at most 20% (see Fig 9). It is therefore reasonably accurate to pre-compute $\varLambda$ at a fixed typical coronal electron density (say $n_{{\rm e}}=10^{16}\,{{\rm m}}^{-3}$).

Fig. 9

Image reproduced with permission from Landi and Landini (1999), copyright by ESO

Sensitivity of the coronal loss function $\varLambda (T,n_{{\rm e}})$ to the electron density. Full line: $10^{14}$ versus $10^{16}\,{{\rm m}}^{-3}$; dashed line: $10^{14}$ versus $10^{18}\,{{\rm m}}^{-3}$; dash-dotted line: $10^{14}$ versus $10^{20}\,{{\rm m}}^{-3}$.

A larger uncertainty can be introduced by inaccuracies of the elemental abundances. The coronal loss function is dominated by losses from C, Si, O, and Fe, and the abundances of these elements have a large influence. The abundance of C and O is debated after 3D non-LTE computations (Asplund et al. 2004, 2005) gave a different result than calculations using 1D models (Grevesse and Sauval 1998).

Another complication is that the most accurate abundances are derived from photospheric lines, but coronal abundances are generally different from abundances in the photosphere. Feldman et al. (1992) compared coronal abundances for 15 elements to photospheric abundances, and Fig. 10 show the difference in the loss function when using photospheric or coronal abundances. Unfortunately it appears that no systematic re-investigation of coronal abundances have been made since 1992. Furthermore, the coronal abundances are not constant, but depend on the coronal structure (see, e.g., Sect. 2.1 of Cranmer and Winebarger 2019, and references therein). Abundances might thus well be the dominant source of uncertainty in the coronal loss function.

Fig. 10

Image reproduced with permission from Dere et al. (2009), copyright by ESO

Radiative loss functions $\varLambda$ computed assuming $n_{{\rm e}}=10^{15}\,{{\rm m}}^{-3}$ as function of temperature. Upper curve: coronal abundances from Feldman et al. (1992). Lower curve: photospheric abundances from Grevesse and Sauval (1998).

The CHIANTI atomic database² (Dere et al. 1997, 2019) provides an extensive compilation of critically assessed atomic data. In addition to the data it delivers software packages in both the IDL and Python languages for using this data and easy generation of loss functions using a variety of abundances and other input options.

Computation of Q in the corona is then straightforward:

$$\begin{aligned} Q = - \varLambda (T) n_{{\rm e}}n_{{\rm H}}, \end{aligned}$$

(50)

where $\varLambda (T)$ is computed using a 1D lookup table. The hydrogen density is easily computed from the mass density, and the electron density can be approximated accurately assuming full ionisation of both H and He at high temperatures. In the low temperature regime, where H i, He i, and He ii exist in significant amounts, one can employ precomputed tables of the electron density assuming LTE or coronal equilibrium.

In order to avoid contribution from the convection zone, photosphere, and chromosphere, one needs to multiply Q with a cutoff function that drives Q to zero in the lower atmosphere. Bifrost employs a soft cutoff function: ${{\rm exp}}(-P/P_0)$, with P the gas pressure and $P_0$ a typical pressure at the top of the chromosphere (Gudiksen et al. 2011). MURaM uses a hard cutoff at T $=$ 20,000 K (Rempel 2017).

The radiative losses in the TR and corona tend to have sharp peak in the transition region owing to the quadratic density-dependence. Rempel (2017) noticed that the relatively large grid spacing of the simulations compared to the thickness of the TR can lead to inaccuracies in the computation of Q. He proposed a scheme using subgrid interpolation to improve the accuracy.

A demonstration of Q and its relation to temperature and density is shown in Fig. 11 for a simulation carried out with the MURaM code. The cooling is largest (note the logarithmic scale) in the transition region, just at the border between the TR and the chromosphere, where it typically is in the range of 0.1–$1.0\,{{\rm W}}\,{{\rm m}}^{-3}$, and it quickly drops down to coronal values around $10^{-4}$–$10^{-5}\,{{\rm W}}\,{{\rm m}}^{-3}$. Zooming in reveals that the strongly-emitting layer is often only a few pixels wide.

Fig. 11

Figure computed from a simulation by Danilovic

Radiative losses in the transition region and corona in a radiation-MHD simulation computed with the MURaM code. Top: temperature; middle: mass density; bottom: $-Q$, i.e., a positive value means radiative cooling. Also note the strong corrugation of the chromosphere-TR boundary.

5 Radiative transfer in the chromosphere

The assumptions underpinning the methods for photospheric radiative transfer as presented in Sect. 3 are no longer valid in the chromosphere: the chromosphere has a low opacity except in a few strong spectral lines, in ultraviolet continua below 160 nm, and in (sub-)mm continua above 160 $\upmu$m. The lines that dominate the radiative energy exchange are the H&K and infrared triplet of Ca ii, the h&k lines of Mg ii, and H i Ly$\alpha$ (see Fig. 12). The opacity in these lines is severely underestimated in the construction of a bin-averaged opacity. The source function is no longer described by the Planck function, and the assumptions of coherent scattering or complete redistribution are very inaccurate for these strong lines, which should instead be modelled using partially coherent scattering (PRD). The ionisation balance of hydrogen and helium is out of equilibrium.

Fig. 12

Image reproduced with permission from Vernazza et al. (1981), copyright by AAS

Net radiative cooling in the 1D semi-empirical VAL3C model atmosphere. The cooling between $z=700$ km and $z=2120$ km in this model is dominated by five lines from Ca ii and two lines from Mg ii. At larger heights H i Ly$\alpha$ alone is the dominant radiative cooling agent.

Proper inclusion of the radiative cooling in the chromosphere (even when ignoring non-equilibrium effects) thus involves solving the 3D non-LTE radiative statistical-equilibrium transfer problem including PRD. This is in principle possible using dedicated radiative transfer codes, but it is computationally expensive. A single solution to the problem for a single atom costs a CPU time of $\sim 10$ s per grid cell (Sukhorukov and Leenaarts 2017), which is very large compared to the $\sim 5$ $\upmu$s CPU time per grid cell per timestep for radiation-MHD simulations (e.g., Gudiksen et al. 2011). Simplifications that speed up the computation are thus required.

Carlsson and Leenaarts (2012) developed techniques to do so by describing the net effect of all the radiative transfer as a combination of (1) an optically thin radiative loss function which represents the local energy loss through radiation per atom in the right ionisation stage per electron, (2) the probability that this energy escapes the atmosphere, and (3) the fraction of atoms in the ionisation stage under consideration. These three factors must all be determined empirically because there are no obvious general physics-based approximations.

The method approximates the radiative loss per volume owing to species X in ionisation state m as

$$\begin{aligned} Q_{X_m} = - L_{X_m} E_{X_m}(\tau ) \frac{n_{X_m}}{n_X} A_X \frac{n_{{{\rm H}}}}{\rho } n_{{{\rm e}}} \rho . \end{aligned}$$

(51)

Here, $L_{X_m}$ is the optically thin radiative loss function per electron and per particle of element X in ionisation stage m, $E_{X_m}(\tau )$ is the photon escape probability as function of some depth parameter $\tau$, $\frac{n_{X_m}}{n_X}$ is the fraction of element X in ionisation stage m, and $A_X$ the abundance.

The quantities $L_{X_m}$, $E_{X_m}$, and $n_{X_m}/n_X$ are determined from a 1D radiation-hydrodynamic simulation including non-equilibrium ionisation computed with the Radyn code (Carlsson and Stein 1992, 2002) for H i. For Ca ii and Mg ii they were determined from a 2D radiation-MHD simulation with Bifrost that provided the atmospheric structure and subsequent statistical equilibrium radiative transfer calculations including PRD using Multi3D (Leenaarts and Carlsson 2009).

The loss function $L_{X_m}$ can be computed for each grid cell in the simulation by summing the net downward radiative rates multiplied with the energy difference of the transition, summed over all relevant transitions. The joint probability density function (JPDF) of $L_{X_m}$ and gas temperature for Ca ii is shown in the top panel of Fig. 13. Radiative transfer effects make that $L_{X_m}$ is no longer a unique function of T, but the red curve indicates an approximate fit. The ionisation degree can be computed from the atomic level populations in the calibrating simulations (see the bottom panel of Fig. 13). Again they are no longer clean functions of temperature, but the spread is minor and a sensible fit as function of temperature can be made. Finally, the empirical escape probability $E_{X_m}$ is computed from the total radiative losses in the simulation, the radiative loss function and the ionisation degree. The middle panel of Fig. 13 shows $E_{{\text {Ca}}\,\textsc {ii}}$ with the vertical column mass as depth parameter. Again there is considerable spread in the JPDF.

Fig. 13

Image reprinted with permission from Carlsson and Leenaarts (2012), copyright by ESO

Chromospheric radiative losses in Ca ii using empirically calibrated recipes. Top: JPDF of radiative loss function and temperature in the chromosphere of a radiation-MHD model. The blue curve is the coronal approximation, the red curve the adopted fit. Ca ii actually heats (meaning a negative loss function, not shown here because of the logarithmic axis) at low temperatures, and that is why the red curve appears a bad fit. Middle: escape probability, with the adopted fit in red. Bottom: fraction of all Ca as Ca ii. Red is the adopted fit, blue the ionisation balance under coronal equilibrium conditions and green the balance assuming LTE.

Figure 14 displays a comparison of this simple recipe with a detailed calculation. The recipe does a surprisingly good job in reproducing the average cooling given the simplicity of the method. However, in individual grid cells large errors in Q might occur, caused by the spread of values around the red curves in Fig. 13.

Fig. 14

Image reprinted with permission from Carlsson and Leenaarts (2012), copyright by ESO

Test of the validity of tabulated chromospheric radiative losses. The curves show the horizontally averaged radiative cooling (i.e., $-Q/\rho$ as a function of height in a 2D Bifrost radiation-MHD simulation. The colored lines show the results computed from a 2D non-LTE radiative transfer computation (solid) and the approximate recipe (dashed) are shown. For hydrogen there is no detailed radiative transfer computation because the recipe was computed from a 1D simulation with the Radyn code. The black solid land dashed lines show the total cooling from all three elements combined.

Computing this chromospheric radiative loss function in a radiation-MHD simulation is fast and simple; it only involves computing the values of the fitted quantities from 1D lookup tables. The electron density can for example be determined from a 2D lookup table computed assuming LTE. This is not particularly accurate but simple. A more realistic electron density can be obtained by computing the non-equilibrium ionisation balance of hydrogen and/or helium (see Sect. 6).

Carlsson and Leenaarts (2012) also describe techniques to model the absorption by the chromosphere of radiation emitted in the corona. Half of this radiation is emitted downwards and the majority will be absorbed in the chromosphere. The corona emits mainly in the far UV regime, and the dominant source of extinction at these wavelengths are the H i, He i, and He ii continua. The coronal loss function is a frequency-integrated quantity (see Sect. 4), so in order to model the extinction one has to estimate a representative opacity. Carlsson and Leenaarts (2012) propose to use the opacity at the edge of the He i continuum at 50 nm. The contribution to the flux divergence is then given by

$$\begin{aligned} Q_{{{\rm abs}}} = 4\pi \kappa _{{{\rm He}\,{\rm I}}} \rho J_{{{\rm cor}}}, \end{aligned}$$

(52)

with $\kappa _{{{\rm He}\,{\rm I}}}$ the representative opacity. The angle-averaged radiation field $J_{{\rm cor}}$ is computed from the transfer equation:

$$\begin{aligned} \frac{{{\rm d}}I}{{{\rm d}}s} = \eta - \kappa _{{{\rm He}\,{\rm I}}} \rho I, \end{aligned}$$

(53)

with the emissivity given by $\eta = - Q_{{{\rm cor}}}/4\pi$. The quantity $Q_{{{\rm cor}}}$ is the coronal loss function as defined in Eq. (50). In other words: it is assumed that only the corona contributes to the production of photons, the chromosphere can only absorb these photons, and scattering is ignored.

Because most radiation-MHD codes already include a 3D radiative transfer module used for the photospheric radiation, one can relatively easily implement the absorption of coronal radiation at a very modest increase in computation time. Simpler 1D recipes that ignore the spreading of radiation in the horizontal directions are also possible. They offer little else besides a marginal increase in computation speed, and can lead to spuriously large heating when coronal radiation emitted by a small localized source is forced to be absorbed in a single vertical column.

The population of neutral helium in the ground state needed to compute $\kappa _{{{\rm He}\,{\rm I}}}$ can be computed assuming LTE populations without a too large error. Computing it taking into account non-equilibrium ionisation is more accurate but much slower (see Sect. 6).

Figure 15 gives an example of how the various approximations for radiative cooling act in different atmospheric regimes. The photospheric cooling rate per mass is largest in the photosphere, but has a small contribution up into the chromosphere. Coronal losses are relevant throughout the entire corona, but are largest in the transition region. The chromospheric losses are largest just below the transition region owing to Ly$\alpha$ cooling (see also Fig. 12). Modest cooling owing to Ca ii and Mg ii lines and heating from the absorption of coronal radiation happens somewhat deeper in the chromosphere. The lower-right panel, which combines photospheric, chromospheric, and coronal losses, clearly shows that the largest radiative losses per mass unit occur in the transition region.

Fig. 15

Example of chromospheric radiative losses and gains in comparison to the photospheric and coronal losses and gains in a 3D simulation computed with the Bifrost code. Top row: temperature and density in a vertical slice through the simulation. The panel labeled $Q_{{{\rm photo}}}$ shows the radiative losses per mass unit computed with the method described in Sect. 3; $Q_{{{\rm chromo}}}$ displays the chromospheric losses and absorption of coronal radiation per mass unit using the methods from Carlsson and Leenaarts (2012); $Q_{{{\rm corona}}}$ shows the losses per mass unit computed as described in Sect. 4. The panel labeled $Q_{{{\rm total}}}$ shows the sum of the previous three panels, i.e., all radiative losses and gains in the simulation. Brown color indicates cooling, blue color represents heating

6 The equation of state and non-equilibrium ionisation

Non-equilibrium radiative transfer and non-equilibrium ionisation play a role in the solar atmosphere wherever magneto-hydrodynamic timescales are short compared to timescales on which an atomic system reacts to changes. The combined continuity and rate equation for the level population $n_i$ in level i of an atomic species with N energy levels is given by:

$$\begin{aligned} \frac{\partial n_i}{\partial t} + \nabla \cdot (n_i{\mathbf{v}} ) = \sum _{j=1,j \ne i}^N n_j P_{ji} - n_i \sum _{j,j \ne i}^N P_{ij}, \end{aligned}$$

(54)

where $P_{ij}$ is the rate coefficient for transitions from level i to level j. If one ignores the advection term and assumes that all the $P_{ij}$ are constant in time then the equations for all levels together form a set of coupled first-order differential equations

$$\begin{aligned} \frac{\partial {\mathbf{n}} }{\partial t} = P{\mathbf{n}} , \end{aligned}$$

(55)

with ${\mathbf{n}}$ the vector of level populations and P the rate coefficient matrix. This system has the solution:

$$\begin{aligned} {\mathbf{n}} (t) = \sum _{i=1}^N c_i {\mathbf{a}} _i {{\rm e}}^{\lambda _i t}, \end{aligned}$$

(56)

where $a_i$ are the eigenvectors of the rate matrix with corresponding eigenvalues $\lambda _i$, and $c_i$ are constants that depend on the initial condition. One of the eigenvalues is zero, and the corresponding eigenvector represents the equilibrium solution. All other eigenvalues are negative and have an absolute value smaller than one. If the system starts away from the equilibrium solution, then it will evolve towards equilibrium on a characteristic timescale given by

$$\begin{aligned} \tau = \frac{1}{\min (|\lambda _i|)}, \lambda _i \ne 0. \end{aligned}$$

(57)

A detailed discussion of the time dependence of atomic level populations is given in Judge (2005). One of the results from that paper is that the slowest time scale tend to be associated with processes that have small net rate coefficients. The smallest rate coefficients are often associated with ionisation and recombination processes (as well as transitions involving metastable levels).

An illustrative solution for a two-level atom was derived in Carlsson and Stein (2002). The equilibration timescale is given by

$$\begin{aligned} \tau = 1{\Big /}\left( C_{21}+C_{12} + R_{21} \left[ 1 - \frac{n_1}{n_2} \frac{R_{12}}{R_{21}}\right] \right) . \end{aligned}$$

(58)

with $C_{12}$ and $C_{21}$ the upward and downward collisional rate coefficient, $R_{12}$ and $R_{21}$ the radiative rate coefficients and the term between square brackets the net radiative bracket. The collisional coefficients depend linearly (or quadratically in case of collisional recombination) on the electron density and exponentially on temperature. For hydrogen one thus expects a maximum in the chromosphere where the temperature and density are relatively low and strong lines are close to detailed balance so that the net radiative bracket is small.

Carlsson and Stein (2002) did a detailed calculation of this timescale for hydrogen and found that the timescale can be as long as $10^5$ s in the chromosphere (see Fig. 16), and that this is the timescale on which ionisation equilibrium is established.

Fig. 16

Image reproduced with permission from Carlsson and Stein (2002), copyright by AAS

Relaxation timescale for hydrogen ionisation as function of column mass in a 1D radiation hydrodynamics simulation computed with the Radyn code. The numerically determined relaxation timescale is given as the thick solid line. Timescales determined from eigenvalues of the rate matrix are also shown for several cases: full rate matrix (dotted line), collisions only (dot-dashed line), Lyman transitions in detailed balance (thin solid line), and Ly$\alpha$ treated with a constant net radiative bracket (dashed line). The electron density is given as the thick dashed curve.

A similar calculation for helium was done for helium by Golding et al. (2014), finding timescales of the order $10^2$–$10^3$ s in the chromosphere and transition region, again associated with the ionisation equilibrium.

These timescale are larger than the hydrodynamical timescales in the chromosphere and transition region. This has severe consequences for the treatment of radiation hydrodynamics. Because hydrogen and helium are majority species they do not only contribute to the radiative flux divergence, but their ionisation state also influences the temperature, pressure and electron density. The assumption that the equation of state (EOS) can be computed assuming LTE is no longer accurate.

The relaxation timescale itself depends on the radiation field, so that a proper treatment of the EOS now requires solution of Eq. (54) together with an equation for charge conservation

$$\begin{aligned} n_{{\rm e}} = \sum _{i,j,k} (j-1) n_{ijk}, \end{aligned}$$

(59)

and energy conservation

$$\begin{aligned} e = \frac{3}{2} k_{{{\rm B}}}T \left( \sum _{i,j,k} n_{ijk}+ n_{{{\rm e}}}\right) + \sum _{i,j,k} n_{ijk}E_{ijk} . \end{aligned}$$

(60)

Here $n_{ijk}$ are the atomic level populations of species i in ionisation state j and excitation state k and $E_{ijk}$ is the sum of the dissociation, ionization and excitation energy of a particle in state ijk. The rate coefficients $P_{ij}$ contain the radiation field so that the transport equation [Eq. (14)] must be solved too. As stated before, this is too computationally expensive to be of use in 3D simulations.

Sollum (1999) developed approximations to the chromospheric radiation field in hydrogen transitions based on detailed 1D calculations with the Radyn code. He found that the angle-averaged and profile-function-averaged radiation field in a given transition can be approximated by a constant above a certain height in the chromosphere, and by the local Planck function at larger depth. In between these two limits he specifies a smooth transition. His recipes allow for an extreme simplification because the solution of the coupled set of Eqs. (54), (59), and (60) can then be computed without having to solve the transfer equation. A limitation of his method is that the Lyman-$\alpha$ line was set in detailed radiative balance. This limitation makes the solution less accurate in the very upper chromosphere and transition region.

His method was implemented in the Oslo Stagger Code by Leenaarts et al. (2007), and used to perform a 2D radiation-MHD simulation of the solar atmosphere with a non-equilibrium EOS. Helium was still treated assuming LTE populations. The temperature and hydrogen ionisation degree in this simulation is demonstrated in Fig. 17. As a consequence of the long equilibration timescale, the ionisation degree in the chromosphere does not follow the temperature. The ionisation degree is instead rather constant in time for a given Lagrangian fluid element.

Fig. 17

Image reproduced with permission from Leenaarts et al. (2007), copyright by ESO

Temperature and hydrogen ionisation degree $n_{{\rm {\text {H}}\,\textsc {i}}} /(n_{{\rm {\text {H}}\,\textsc {i}}} + n_{{\rm {\text {H}}\,\textsc {ii}}})$ as function of time in two columns of a 2D radiation-MHD simulation that includes non-equilibrium ionisation of hydrogen. Left-hand panel: a column in a magnetic element. Right-hand panel: a column with weak magnetic field. The left-hand panel shows regular shock waves with 3-min period, while the right-hand panel shows a more irregular temperature structure. In both cases the ionisation degree does not follow the temperature structure. Instead it is rather constant.

Golding et al. (2016) developed approximations for the radiation field in helium transitions. The physics of helium ionisation is somewhat more complicated than for hydrogen because ionisation in the chromosphere is mainly driven by UV radiation produced in the corona. The resulting recipes take this into account, but require the solution of the transfer equation in seven radiation bins in order to approximate the radiation field in the continua of helium as well as the He ii 30.4 nm line.

These authors also implemented a semblance of radiative transfer in the Ly $\alpha$ line. They assigned a single-bin source function based on the net radiative rate and representative opacity to the line, and used those to solve the transfer equation ignoring scattering. The resulting radiation field is then used to add an upward radiative rate coefficient in Eq. (54). In this way, the recipes by Sollum (1999) were extended to the very upper chromosphere, but in a rather crude fashion.

The increased realism of an EOS that includes non-equilibrium ionisation comes at a rather steep price in computational efficiency: a simulation with a non-equilibrium EOS is around three (hydrogen only) to five times (hydrogen and helium) slower than a similar simulation with an LTE EOS.

Non-equilibrium ionisation has strong consequences on the structure of the chromosphere and transition region. Most importantly, ionisation can no longer function as efficiently as an energy buffer when the internal energy density of a gas parcel changes.

In LTE energy must be expended to ionise hydrogen and helium if the internal energy density is increased before the temperature can rise. Likewise, a decrease in internal energy leads to an instantaneous transfer of ionisation energy to thermal energy, slowing down the temperature decrease. This leads to characteristic bands of “preferred temperatures” in joint probability density functions of radiation-MHD simulations assuming LTE (see Fig. 18). The temperature of these bands are associated with the temperatures where H i, He i, and He ii ionise according to the Saha–Boltzmann equations.

Fig. 18

Image reprinted with permission from Golding et al. (2016), copyright by AAS

Joint probability density functions of height and temperature in four radiation-MHD simulations of the solar atmosphere. The simulations differ in their treatment of the equation of state. LTE LTE equation of state, HION hydrogen in non-equilibrium, helium in LTE, hydrogen Lyman transitions in detailed balance. LYA-HION as HION but with radiative transfer in Lyman transitions, HELIUM both hydrogen and helium in non-equilibrium. The three horizontal plateaus (at 6, 10, and 22 kK) in the LTE simulation indicate preferred temperatures when using the LTE equation of state. These temperatures are associated with the LTE ionisation of H i, He i, and He ii. The plateaus vanish when we introduce non-equilibrium hydrogen and helium ionisation.

If the ionisation balance is computed in non-equilibrium, and the ionisation/recombination timescale is long, then an increase in internal energy leads directly to a temperature increase. Temperature decreases are also stronger than in LTE because ionisation energy cannot be released quickly enough to counteract cooling. The bands of preferred temperature disappear, and the gas in the chromosphere behaves somewhat like and ideal gas. A clear example of this effect is the increased the amplitude of the temperature jump in acoustic shocks (see Carlsson and Stein 2002; Leenaarts et al. 2007).

Low-temperature areas in the chromosphere have a higher electron density than predicted by the Saha–Boltzmann equations, and the reverse is true for high-temperature areas. This has an effect on the formation of chromospheric lines because the source function couples to the local temperature mainly through collisions with electrons. Non-equilibrium ionisation is also expected to have an effect on the efficiency of heating through ambipolar diffusion (e.g., Martínez-Sykora et al. 2017; Khomenko et al. 2018).

7 Other developments

7.1 Fast approximate radiative transfer in the photosphere

Abbett and Fisher (2012) propose a method for quick evaluation of the radiative losses in the photosphere as an alternative to the methods described in Sect. 3. They assume that each column in a simulation can be treated as an independent plane-parallel atmosphere. The angle-averaged radiation field at vertical optical depth $\tau _\nu$ in a column is then given by:

$$\begin{aligned} J_\nu (\tau _\nu ) = \frac{1}{2} \int _0^\infty S_\nu (\tau _\nu ') E_1 \left( |\tau _\nu -\tau _\nu ' | \right) {{\rm d}}\tau _\nu '. \end{aligned}$$

(61)

The first exponential integral $E_1(x)$ peaks sharply at $x=0$, so that this expression can be approximated as

$$\begin{aligned} J_\nu (\tau _\nu ) \approx \frac{1}{2} S_\nu (\tau _\nu ) \int _0^\infty E_1 \left( |\tau _\nu -\tau _\nu ' | \right) {{\rm d}}\tau _\nu ' \end{aligned}$$

(62)

$$\begin{aligned} \approx S_\nu (\tau _\nu ) \left( 1- \frac{E_2(\tau _\nu )}{2} \right) \end{aligned}$$

(63)

This result is inserted into the equation for the radiative cooling [Eq. (11)] to yield

$$\begin{aligned} Q \approx - 2 \pi \rho \int _0^\infty \kappa _\nu S_\nu E_2(\tau _\nu ) \, {{\rm d}}\nu . \end{aligned}$$

(64)

If one assumes an LTE source function, then the frequency integral can be approximated using a reasoning similar to the one given in Sect. 3.2, and together with some additional assumptions one arrives at the final result:

$$\begin{aligned} Q \approx - 2 \kappa ^B \rho \sigma T^4 E_2(\tau ^B), \end{aligned}$$

(65)

where $\kappa ^B$ is the Planck-averaged opacity (see Sect. 3.2), $\tau ^B$ the optical depth computed from this opacity, and $\sigma$ is the Stefan–Boltzmann constant. With this scheme, computation of the radiative heating only requires a simple 2D table lookup to get $\kappa ^B$, and a column-by-column integration over depth to compute $\tau ^B$.

Figure 19 shows the resulting temperature structure in the photosphere in a simulation where the radiative losses are computed in this simplified fashion. The method is simple and extremely fast and is well suited for problems that do not require very accurate radiative losses, but nevertheless want the radiative losses to drive reasonably realistic-looking convection.

Fig. 19

Image reproduced with permission from Abbett and Fisher (2012), copyright by Springer

Temperature at a constant height in the photosphere of a 3D simulation computed with the RadMHD code. The panel spans a size of $24 \times 24\,{{\rm Mm}}^2$ and has a grid spacing of 21.3 km. The convection pattern strongly resembles the result of simulations that model radiation with higher fidelity, despite the strong simplifications in the computation of the radiative losses.

7.2 Escape probability method

An interesting new method that speeds up non-equilibrium radiative transfer was proposed by Judge (2017). The major time consuming task in radiative transfer is the evaluation of the formal solution for all required frequencies and angles in order to construct the angle-averaged and frequency-averaged radiation field

$$\begin{aligned} {\bar{J}} = \frac{1}{4 \pi } \int _0^\infty \int _\varOmega I_\nu \, {{\rm d}}\varOmega \, {{\rm d}}\nu , \end{aligned}$$

(66)

which appears in the radiative rate coefficients of Eq. (54).

Judge (2017) suggests to compute ${\bar{J}}$ from the equation

$$\begin{aligned} \frac{{{\rm d}}}{{{\rm d}}\tau } \left( S-{\bar{J}}\right) = q^{1/2} \frac{{{\rm d}}}{{{\rm d}}\tau } \left( q^{1/2} S \right) , \end{aligned}$$

(67)

(see Frisch and Frisch 1975; Hummer and Rybicki 1982), where $\tau$ is a vertical optical depth parameter and the function q an escape probability function that only depends on the vertical optical depth. This equation is approximate only: the main approximations are that the source function varies only slowly over an optical path length and that horizontal structure in the atmosphere can be ignored.

The big advantage of using Eq. (67) is that it replaces the repeated evaluations of the transfer equation in order to compute ${\bar{J}}$ in a transition with the solution of a single integral along vertical columns in the atmosphere. Solving the statistical equilibrium non-LTE radiative transfer problem in a MURaM test atmosphere is $\sim 100$ times faster than the full method.

Judge (2017) shows that this method can also be used to solve non-equilibrium problems (i.e, it can be used to simultaneously solve Eqs. (54), (59), and (60)). Radiation-MHD simulations using this method to compute chromospheric radiative energy exchange or the equation state have so far not been reported on. It would be very interesting to see whether the method is fast enough to be used in practice and how it compares to the methods described in Sect. 6.

8 Conclusions and outlook

In this review, I presented the most commonly used methods to approximate the transfer of energy between solar plasma and the radiation field in radiation-MHD simulations. The theory and methods for computing radiative energy exchange in the photosphere are perhaps the most well-developed and well-studied. They are also the most accurate: Pereira et al. (2013) compared a 3D hydrodynamics model of the upper solar convection zone that employed accurate abundances and opacity calculations, together with multi-group LTE radiative transfer with 11 groups. The resulting model reproduces the observed continuum center-to-limb intensity variation, the absolute flux spectrum, the wings of hydrogen lines and the distribution of continuum intensities caused by granulation to a high degree of precision. It seems therefore safe to say that we can model photospheric radiative energy exchange sufficiently accurate for almost all purposes.

Radiative losses in the transition region and corona are much less accurate when using the standard method of an optically loss curve computed from statistical equilibrium, no absorption of radiation and a fixed set of abundances as described in Sect. 4. Inaccuracies caused by ignoring non-equilibrium ionisation effects appear to be largest in the transition region and during solar flares, and might well lead to a factor two error in the value of the loss function $\varLambda (T)$. The choice of abundance has a smaller effect in the transition region, but can lead to differences up to a factor three at higher temperatures in the corona. A systematic test of these effects in 3D radiation-MHD simulations has so far not been done.

Radiative transfer and non-equilibrium ionisation in the chromosphere is probably the least studied and the methods described in Sects. 5 and 6 are likely the most inaccurate of all method used in the photosphere-chromosphere-corona system. A long sequence of approximations is needed to order to make the problem computationally feasible. Testing the accuracy of all approximations combined can only be done in 1D geometry (using Radyn or similar codes), and has so far not been done. A particular worry is the crude way in which the radiative transfer in the Ly-$\alpha$ line is implemented in the non-equilibrium ionisation method,

Interestingly, the way how chromospheric radiative transfer is implemented, and the accuracy of the method used seems however to have only a minor effect on the overall structure of the simulated chromospheres. A MURaM simulation of active regions that only contains single-bin (gray) LTE radiative transfer and a coronal loss function produce chromospheres whose structure resembles the real chromosphere (see, e.g., Bjørgen et al. 2019 and Fig. 20).

Fig. 20

Brightness temperature in the nominal line core of H$\alpha$ in a simulated active region. The radiation-MHD simulation was performed with the MURaM code by Danilovic. The simulation used a single-bin LTE treatment of chromospheric radiative transfer. Nevertheless, the simulation resembles observations rather well. The 3D non-LTE radiative transfer calculation used to obtain the H$\alpha$ intensity was done with the Multi3D code

I argue that this is a consequence of the physics of the chromosphere in the MHD approximation. Thermal conduction in the chromosphere is not efficient, and the only way a Lagrangian fluid element can lose energy is through radiation. In a chromosphere that is in a statistically stable state, the dissipation of non-thermal energy into heating and the radiative cooling must be in balance in a volume and time-integrated sense. An increase in non-thermal energy deposition must be accompanied by an increase in radiative losses. In other words: radiative losses are a reaction to heating, and the chromosphere will adapt its thermodynamic state until heating and cooling are in balance again.

Here is the catch: the radiative losses are very sensitive to temperature, but the mechanisms that dissipate non-thermal energy are not. The radiative losses assuming LTE scale as $T^4$, while the radiative losses in the coronal equilibrium approximation scale exponentially with temperature under chromospheric conditions (see Fig. 10). The non-LTE chromospheric loss function will lie between these two extremes most of the time.

In the MHD approximation there are only two mechanisms that irreversibly convert non-thermal energy into heat: viscosity and electrical resistance. The first one scales as $T^{1/2}$ (assuming a dilute gas of rigid elastic spheres), while the second scales as $\ln T/T^{-3/2}$ (assuming Spitzer resistivity). In practice the viscosity and resistance are replaced by numerical terms that are independent of temperature in almost all numerical simulations.

When the description of the chromospheric radiative losses or the equation of state is not correct, then the simulation will change the temperature compared to the correct solution until balance between energy gains and losses is achieved again. The steep temperature dependence of the radiative losses makes the temperature change modest (say 1000–2000 K, if our methods are not too bad), but the effect on the non-thermal energy dissipation rate is small. The result is a chromosphere with the wrong temperature, but almost correct dissipation, density and velocity structure.

It thus seems that one can study certain aspects of chromospheric physics without a too complicated treatment of the radiative losses and equation of state.³ This does, however, not mean that the work described in Sects. 5 and 6 can be ignored. The only proper way to validate models is to compute the various diagnostics (line profiles and continua) and compare to observations. The emergent intensities in the diagnostics sensitively depend on the temperature and electron density. Comparison of models with observations (such as in Leenaarts et al. 2013; Bjørgen et al. 2018) show that current models are not getting the intensities right, and that they need to be refined. The treatment of chromospheric radiative energy exchange is very likely one of the aspects in need of further improvement.

The reverse comparison is also true: inferred model atmospheres generated through non-LTE inversions of observations (e.g., de la Cruz Rodríguez and van Noort 2017; de la Cruz Rodríguez et al. 2019) can only be used to constrain physical models of the chromosphere if we are sure that the treatment of radiation in the models is done sufficiently accurately.

So what can be improved? I propose the following non-exhaustive list:

• It would be interesting to test, and if necessary improve, the accuracy of Skartlien’s multi-bin radiative transfer with scattering in the mid and upper chromosphere. If successful, this would increase the accuracy compared to the radiative loss tables of Carlsson and Leenaarts (2012), which ignore 3D effects.

• The tables for Ca ii and Mg ii of Carlsson and Leenaarts (2012) should be updated. They were calibrated on a single 2D radiation-MHD simulation with a weak magnetic field using 1.5D radiative transfer. A much larger variety of models is available now, and 3D non-LTE radiative including PRD is now possible (Sukhorukov and Leenaarts 2017).

• The accuracy of the non-equilibrium ionisation methods should be more critically assessed. Comparison against the full physics is only possible in 1D but this will already give more insights in the possible deficiencies of the model, and in particular of the treatment of Ly$\alpha$.

• The escape probability method for approximating the chromospheric radiation field discussed in Sect. 7.2 should be tested and developed further.

• It should be investigated whether non-equilibrium ionisation of elements that are important for radiative losses in the TR and corona can be implemented in a sufficiently computationally efficient fashion. The ionisation balance of helium can be modeled with only one energy level per ionisation stage and parametrized ionisation/recombination rates that include the effects of the excited levels. If a similar thing can be done for other elements, then inclusion of non-equilibrium effects is otherwise straightforward.

Let us hope that improvements along these lines, as well as others, will be presented during the coming years.