Radiation hydrodynamics in simulations of the solar atmosphere

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.


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 ∼ 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 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.
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 Sections 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.

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.

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 with ρ the mass density, v the velocity vector, p the momentum density, τ the stress tensor, P the gas pressure, J the current density, B the magnetic field vector, g the acceleration due to gravity, P rad the radiation pressure tensor, e the internal energy, Q 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. 2017Khomenko et al. , 2018 These inclusions retain the single fluid description. This greatly simplifies the treatment of the radiative term Q 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 multifluid models.

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 rad = 5777 K. Then the energy density is with σ 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 F rad = σT 4 rad = 6.3 × 10 7 W m −2 , 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 = nk 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 chromospere, 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.

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 4σ 3c T 4 = 0.27 Pa (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: with κ 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 .

Energy exchange between radiation and matter
In absence of an absorbing medium, the monochromatic radiative flux divergence ∇ · F ν 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 with 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.

Explicit expression of the radiative flux divergence
The monochromatic flux is defined in terms of the intensity I ν as with n the unit vector pointing in the direction of Ω. 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 ν = j ν /(κ ν ρ), intensity and opacity: Radiation hydrodynamics in simulations of the solar atmosphere 7 with κ ν the extinction coefficient per unit mass. If one assumes that both the emissivity and extinction coefficient do not depend on direction then Eqs. (10) reduce to where 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 2kT /m + v 2 turb , with m the mass of the atom or ion that is involved, and v turb the microturbulent velocity. Nevertheless, the assumption is almost always taken to be valid because it allows large simplifications.

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×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 ∼ 5 minutes, it is ∼ 1 minute 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 n then simplifies from with j ν the emissivity and α ν = κ ν ρ 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.

Diffusion approximation
At very large optical depth the following approximation holds for J ν : where the second derivative should be taken along the direction of the gradient of S ν 1 . For the flux the following approximation holds: 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 from Eq. (16) if one also assumes that the source function is equal to the Planck function B ν : with the Rosseland opacity defined as 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): with κ 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 Λ-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 Λ ν [. . .] because the opacity changes with frequency, a better notation would be Λ κ [. . .]. For two frequencies ν 1 and ν 2 with identical opacities everywhere in the atmosphere κ 1 = κ 2 = κ one can thus write In order to derive the approximation in Eq. (19) one first approximates the frequency integral with a sum over discrete frequencies ν 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: 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 frequencydependent source function using 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 κ 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.

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 Fig. 1 Illustration of the principle of τ -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 τ 1 R -τ 3 R . Adapted with permission from Ludwig (1992). 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 τ -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 τ R : Let the borders be defined at a number of Rosseland optical depths τ k R . A frequency ν j then belongs to bin i if A common choice is to put the border values τ k R equidistantly spaced in log τ R . The method is illustrated in Fig. 1 ex j(i, l) runs over those B∆ν i denotes the average effect of the sorting prostrated in Figs. 2 and 3. y from the ATLAS9 stel-T = 4470 K and log p = optical depth at 500 nm ). After applying the tauns the spectrum is a step te values (Fig. 3). pproach the exact opacsmall ODF intervals, the ed to converge to the exn the detailed reference ing process) in the limit s. This is due to the fact s between bin-levels τ l ref n will comprise frequenum with different height s. It should be noted that  using ODFs as basis for sorting frequencies introduces an additional error since the rearrangement of spectral lines inherent in the ODF concept affects the way frequencies are classified. However, since the ODF solution approximates the exact solution well, the error incurred this way is small compared to the consequence of the binning procedure (see Sect. 3.5).
The choice of a reference atmosphere is another possible source of error connected with the τ-sorting scheme. It is plausible that the opacity binning method with τ-sorting shows good results in calculations of static 1D model atmospheres if the "exact" solution (i.e., the 1D atmosphere resulting from the ODF approach) is chosen as the reference atmosphere. In multidimensional time-dependent simulations, however, the physical parameters may deviate considerably from a 1D reference stratification. As a consequence, the assignment of frequencies to bins might lead to an inappropriate representation of the opacities. In order to test the applicability of the τ-sorting procedure in non-planeparallel cases, we performed test τ ref Fig. 1. Schematic illustration of the τ-sorting procedure using ODFs. Within each ODF interval ∆νi, the ODF steps ∆νij are sorted to frequency bins Ωl, depending on the height where τij = 1 is reached (indicated by bold arrows). within one bin are similar, which is the major requirement for the validity of the multigroup approach. For greater depths, the way frequencies are gathered may be inappropriate, but this shortcoming is less severe as radiation transfer becomes increasingly grey in these regions anyway.
The τ-sorting procedure can be realized in a convenient way by using ODFs. For a given step j within ODF interval ∆νi κijϱ is integrated along z, which gives τij as a function of height in the reference atmosphere (see Fig. 1). The ODF steps are then sorted into bins Ωl. Bin-integrated quantities are obtained as weighted sums over the ODF steps in a given bin, viz. and where, for a given bin index l, the index j(i, l) runs over those steps in ∆νi, which are elements of Ωl. B∆ν i denotes the average of B over the ODF interval ∆νi. The effect of the sorting procedure on the opacity spectrum is illustrated in Figs. 2 and 3. Figure 2 shows the ODF-based opacity from the ATLAS9 stellar atmosphere code (Kurucz 1993) for T = 4470 K and log p = 4.12 (corresponding to a continuum optical depth at 500 nm log τ500 = −2 in the solar atmosphere). After applying the tausorting procedure with five opacity bins the spectrum is a step function which can assume five discrete values (Fig. 3).
In contrast to the ODFs, which approach the exact opacity spectrum in the limit of infinitely small ODF intervals, the multigroup solution can not be expected to converge to the exact solution (i.e. the solution based on the detailed reference spectrum which was used in the sorting process) in the limit of a very large number of opacity bins. This is due to the fact that, even for infinitely small intervals between bin-levels τ l ref (or κ l in the case of κ-sorting), each bin will comprise frequencies from different parts of the spectrum with different height profiles of the corresponding opacities. It should be noted that Fig. 2. Monochromatic opacity as a function of wavelength for T = 4470 K and log p = 4.12, corresponding to log τ500 = −2 in the solar atmosphere. ATLAS9 (Kurucz 1993) Opacity Distribution Functions are used to model the line opacities. using ODFs as basis for sorting frequencies introduces an additional error since the rearrangement of spectral lines inherent in the ODF concept affects the way frequencies are classified. However, since the ODF solution approximates the exact solution well, the error incurred this way is small compared to the consequence of the binning procedure (see Sect. 3.5).
The choice of a reference atmosphere is another possible source of error connected with the τ-sorting scheme. It is plausible that the opacity binning method with τ-sorting shows good results in calculations of static 1D model atmospheres if the "exact" solution (i.e., the 1D atmosphere resulting from the ODF approach) is chosen as the reference atmosphere. In multidimensional time-dependent simulations, however, the physical parameters may deviate considerably from a 1D reference stratification. As a consequence, the assignment of frequencies to bins might lead to an inappropriate representation of the opacities. In order to test the applicability of the τ-sorting procedure in non-planeparallel cases, we performed test Fig. 2 Illustration of the concept of group mean opacity. Left-hand panel: Monochromatic extinction coefficients based on Opacity Distribution Functions for T = 4470 K and P = 1.3 × 10 4 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 continous 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. Adapted with permission from Vögler et al. (2004), copyright by ESO.
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 κ 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.

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 km < z < 200 km, with the z = 0 defined as the location where τ 500 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 This expression only depends on temperature and can thus easily be precomputed and stored in a 1D lookup table.
Once the opacities κ j are grouped into the N bins, one still needs to compute an appropriate bin opacity κ i . Choosing the numerical equivalent of the Rosseland opacity in each bin ensures that the diffusion approximation is recovered deep in the atmosphere: 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 κ B is a good choice for the outermost layers of the atmosphere: 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 where τ 0 is an adjustable parameter of order unity and τ i the vertical optical depth. Computing τ 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: 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, κ i depends only on a combination of any two thermodynamic parameters (for example e and ρ). Like the groupintegrated 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).

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 for example Fig. 49 of Vernazza et al. 1981). These lines can have photon destruction probabilities 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 ∇ · F ν /ρ. In LTE it is given by If one assumes the presence of a coherently scattering line in a two-level atom at frequency ν, then the source function becomes In that case the flux divergence per mass unit is so that for a given difference between B ν and J ν , 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 ν = (1 − ν )J ν + ν B ν and following a reasoning similar to the one given in Sect. 3.2, he arrives at an expression for the group source function: Here i represents a group-averaged scattering probability, and t i the groupintegrated thermal production of photons. He also derives an expression for the group extinction coefficient κ i . Similar to the LTE case, i , t i , and κ 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 ν 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 rad now requires solution of Eq. (37) together with the transfer equation because J i = Λ[S i ]. This is typically done through accelerated Λ-iteration. Because of the multidimensional geometry, a local approximate Λ-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).

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), wich 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;Stěpán and Trujillo Bueno 2013). I will only briefly touch upon some aspects.
Solving for the flux 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.
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: where the optical thickness scale has its zero point at the start of the SC and τ is the optical thickness along the entire SC. Intensities I(τ ) are computed at the grid points (grey circles in Fig. 3). Computation of I(τ = 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 In case of actual numerical computation this is no longer the case.
In the deep atmosphere the diffusion approximation holds. Eq. (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πκ i ρ 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 which in its integral form is given by The total heating rate is then computed from This equation does not suffer from the numerical precision problems caused by the cancellation of the subtraction S and I. Eq. (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 ∇ · 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: where 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 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). Figure 4 shows a summary table of the three major binary choices in multigroup 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. The simulation by Nordlund (1982) is of Type 5. MURaM (Vögler 2004;Rempel 2017), RAMENS (Iijima and Yokoyama 2015; Iijima 2016), and MAN-CHA3D (Khomenko et al. 2017(Khomenko et al. , 2018 are Type 2 codes. Stellarbox (Wray et al. 2015) is Type 6. COBOLD (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 2 (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.

Summary and examples of photospheric radiative transfer
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). 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  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 τ i = 1 height in each bin. Note that the maxima of the heating and cooling per unit mass do not coincide with the τ = 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.
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 κ i and the size of (S i − J i ) only. The lower panel shows the amplitude of the approximate and exact Ñux divergence (horizontal average of the absolute value) relative to the amplitude of the total exact Ñux divergence, i.e., and

SKARTLIEN
o T (brackets denote horizontal average). We see that the exact relative amplitudes (thick lines) coincides very well with the approximate values (thin lines).
A sample of the spatial structure is displayed in Figure 4, where we have shown the exact heating per mass unit in all groups as gray-scale images in a vertical slice. Black contours mark locations of zero heating, and lighter shades of gray means positive heating. Gas in layers immediately above the cooling layer is heated in all groups in expanding Ñow above granules. Granules are seen as curved horizontal structures. As up-Ñowing gas expands and cools, the temperature falls below the radiation temperature, and the gas is radiatively heated. Note also heating of the cool region FIG. 4.ÈExact group radiative heating per mass unit through the simulation. Full drawn black contours show heating. Heating is found at lighter shades of gray, and c shades. Dash-dotted black contours show the zero level of group heating and coincide well with the zero level of t White curves are the horizontal averages of the radiative unit (normalized to Ðt the plotting window). Black vertic zero level for these curves, and positive values are to the L ower panel : Temperature in the same vertical slice. No the granular layer at the height 0.0 Mm, and heating imm all groups, and also heating of the cool region below chromosphere. below 3000 K in the chromosphere as rad from below is converted to thermal energy via The dash-dotted black contours show the the approximate heating and coincide well w 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. Adapted with permission from Skartlien (2000), copyright by AAS. 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 = Σ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.
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: where i = 1 in LTE, and i ≤ 1 in non-LTE. In the latter case 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 i = 1 in LTE. 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.

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-equilibriumnon-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 ν = 0 in the rate equations) alleviates  tering would imply that The values of the mean J i * \ S i *. intensities are also di †ering, mainly because of di †erences in the source function, and partially because of the di †erent heights for which the optical depths are unity.
Diamonds and triangles show the heights for where the horizontally averaged optical depth is unity, for previous q i and current methods, respectively. The mean intensities in optically thin regions are approximately SJ i *T^1 2 SS i *T The group mean opacities that determines the (q i \ 1). optical depths are seen in the bottom panel of Figure 6. The new opacities are higher than the old opacities above 0.3 Fig. 8 Horizontal averages of radiative heating per mass unit (Q/ρ 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. Adapted with permission from Skartlien (2000), copyright by AAS.
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 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 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 a spectral line of element X in ionisation stage m with upper level j and lower level i can then be written as: ≡ G(T, n e ) n e n H .
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 H is the abundance of element X relative to hydrogen. The function G(T, n 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 e n m ∼ 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 equilibrium, while for lower densities collisional de-excitation becomes negligible compared to radiative decay and the Coronal Model Approximation (yielding density insensitive line Contribution Functions) may be adopted.
(with N i e = 10 10 , 10 12 and 10 14 cm −3 ) between total emissivity curves calculated at different densities as a function of electron temperature. As expected, the greatest differences are found with the curves at 10 12 − −10 14 cm −3 , which are very similar, because density-dependence affects line emissivity mostly between 10 8 and 10 10 cm −3 . Differences are always smaller than 25% and show a marked temperature dependence, being highest at transition region and coronal temperatures and decreasing down to zero at the edges of the selected temperature range.
The maximum at coronal temperatures is given by the presence of a host of strong density dependent lines formed in quiet corona, mainly from Fe, Mg and Si ions. The high temperature tail is dominated by strong, density insensitive lines and freefree continuum; the low temperaure tail is dominated by density insensitive transition region and chromospheric lines and for this reason there are small differences between computations carried out assuming different density values.

Effect of different datasets and approximations in level population computation
Level populations are strongly sensitive to any change or problem in the atomic parameters, collision strengths and transition probabilities as well as in the approximation adopted for their calculations, and this affects line radiation. It is therefore important to check the effects of different transition probabilities datasets on the resulting total emissivity curve.  6. Percentual difference between the total emissivit tained with the old and new version of the Arcetri Spectra adopted electron density is 10 10 cm −3 .
As big improvements have been done in the p sion of the Code versus the older version described & Monsignori Fossi 1990, we have performed a c between the present results and those obtained usin version of the Arcetri Code. The adopted element a are from Allen 1973. There are three main differenc the two versions of the Arcetri Code: (a) the old calculated all line intensities using the Coronal Mod imation, (b) the collision rates were calculated using tors and (c) radiative data came from different literat than in the present version of the Code.
Thus, the present comparison allows to check also of different assumptions in level population calculat resulting total plasma emissivity. Fig. 6 displays the percentual difference between the two versions of the Code as a function temperature. It is possible to see that rather high diff to 60%) are found at transition region temperatures, a discrepancies occur at coronal temperatures. In the p tion of Fig. 6 the older version of the Arcetri Code total emissivity than the more recent version at transi temperatures. This is due to the presence of few transitions from O iv, O v, C iv whose emissivities different values in the two versions of the Code; thei is due both to the use of different datasets and to th approximations used in level population calculation an overestimation of line emissivity for these transi old version of the Code. The negative section of t is due to the much larger number of lines included version. in G(T, n e ) through collisional deexcitation and density-dependent dielectronic recombination (Summers 1972(Summers , 1974. A total coronal radiative loss function Λ(T, n 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 Λ(T, n e ) is weak: Landi and Landini (1999) performed a sensitivity study for electron densities in the range 10 14 -10 20 m −3 , and found a difference in the value of Λ of at most 20% (see Fig 9). It is therefore reasonably accurate to pre-compute Λ at a fixed typical coronal electron density (say n e = 10 16 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(Asplund et al. , 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 for example Sec 2.1 46% lower in the case of oxygen, with the changes largely resulting from the application of 3D model atmospheres to the interpretation of photospheric absorption lines instead of 1D models. Since solar photospheric abundances are a vital reference point for many fields in astronomy the new abundances have caused a re-assessment of many previous works. The most dramatic has been in the study of the solar interior, where models had previously been in excellent agreement with the precise helioseismic measurements of parameters such as the sound speed profile. Oxygen is the dominant contributor to opacity in the region just below the convection zone, and the revised abundance significantly affected the solar structure in this region, breaking the excellent agreement with measurements. The crucial role of oxygen in the solar interior models has led to a number of different measurements and re-analyses of oxygen lines and photospheric models (e.g., Centeno & Socas-Navarro 2008; Meléndez & Asplund 2008;Caffau et al. 2008, and references therein) leading to oxygen abundances ranging between the value given in Grevesse & Sauval (1998) and that given in Grevesse et al. (2007). This makes it difficult for a project such as CHIANTI which seeks to provide a single set of reference abundances for users. Our solution is to retain the original Grevesse & Sauval (1998) abundances as the default for the CHIANTI software, and to provide the Grevesse et al. (2007) abundances as an option that can be selected by the user.
Although the change in the oxygen abundance is the most significant in the Grevesse et al. (2007) tabulation, it is also to be noted that the neon and argon abundances have changed by a similar amount. This is on account of the fact that neither of these elements can be measured in the solar photosphere, and so they have usually been measured relative to oxygen by other means (spectroscopy of the solar upper atmosphere, solar energetic particle measurements).
CHIANTI users are encouraged to bear in mind the ongoing debate over solar photospheric abundances and to assess the impact on their analyses.

Radiative loss rates
emission measure ( n e N H dV) has been calcul ization equilibria discussed in Sect. 4. The loss processes of bremsstrahlung, radiative recomb ation and two-photon radiation. An electron de was used and two sets of elemental abundanc abundances of Feldman et al. (1992) and the ph dances of Grevesse & Sauval (1998). The valu loss rate are shown in Fig. 5 and in Table 1.
The most recent calculation of radiativ those of Landi & Landini (1999); Colgan et & Landini (1999) used the Arcetri spectra CHIANTI 2 database. Their Fig. 1 shows calc sity of 10 10 cm −3 for the coronal abundances us et al. 1992) and the ionization equilibria of Arn (1985) and Arnaud & Raymond (1992). A co that the calculations of Landi & Landini (1999 to our current values. The major differences where the current loss rate is a factor of 1 1.6 × 10 4 K where the current loss rate is a than the values of Landi & Landini (1999). The higher temperature is primarily due to the fact did not include transitions for the hydrogen-lik ions. The differences at the lower temperature to differences in the way the Arcetri code and late the radiative losses due to the continuum.
There are significant differences between t Colgan et al. (2008) and those of both Landi & and our current values. Detailed comparisons w (2008) are not possible as none of the atomic authors' radiative loss calculation were publis note that for most ions radiative losses are small number of strong fine structure transiti configurations, and so modeling these transitio atomic data is crucial to deriving accurate ra mates. Colgan et al. (2008) treated entire config levels in their calculations with no fine structu no attempt was made by the authors to asses their approach for the dominant low-lying confi CHIANTI fine structure levels are used for e most accurate atomic data for transitions betw are chosen following assessments of the best av literature. We are thus confident in the accuracy radiative losses for the ion models used in the d

New IDL procedures
The parameters that are necessary to calculate sections and rate coefficients, recombination and ionization equilibria are now included i database. Accordingly, new Interactive Data functions ioniz_cross, ioniz_rate, re make_ioneq_all are also provided.
ioniz_cross returns the ionization cross se fied ion and electron energy in eV. ioniz_rate returns the ionization rate coeffi fied ion and temperature in K. recomb_rate returns the recombination rate specified ion and temperature in K make_ioneq_all creates a new ionization eq The CHIANTI atomic database 3 (Dere et al. 1997(Dere et al. , 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: where Λ(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: 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 Figure 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 W m −3 , and it quickly drops down to coronal values around 10 −4 -10 −5 W m −3 . Zooming in reveals that the strongly-emitting layer is often only a few pixels wide. The calculated net radiative cooling rates for the Ca n and Mg n resonance Unes are quite sensitive to the convergence and internal consistency of our numerical solution. Some of our preliminary solutions gave maximum rates for these lines 2 to 3 times larger than those shown in Figure 49.
In Table 29 we list the integrated net cooling rates fïïdh for the Ca n and Mg n lines, La, and H", i.e., the area under each $ as a function oî h in Figure 49. Negative values of 0 in the temperature-minimum region are assumed to be zero in each integral. Athay (1976, Table IX-1) estimates that the contributions to dh have the following order of importance: Balmer continuum, Ha, Ca n IR triplet, Mg n, Ca n K and H, Fe i lines, and Na i D. Athay states, however, that these estimates are not based on detailed computations of the net energy loss or gain in the bound-free continua. The order we determine is Usted in Table 29. © American Astronomical Society • Provided by the NASA Astrophysics Data System Fig. 12 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α alone is the dominant radiative cooling agent. Adapted with permission TODO:ASK PERMISSION from Vernazza et al. (1981), copyright by AAS.

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 µ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α (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 equilbrium.
Proper inclusion of the radiative cooling in the chromosphere (even when ignoring non-equilibrium effects) thus involves solving the 3D non-LTE radiative statistical-equilbrium 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 ∼ 10 s per grid cell (Sukhorukov and Leenaarts 2017), which is very large compared to the ∼ 5 µ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 Here, L Xm is the optically thin radiative loss function per electron and per particle of element X in ionisation stage m, E Xm (τ ) is the photon escape probability as function of some depth parameter τ , n Xm n X is the fraction of element X in ionisation stage m, and A X the abundance.
The quantities L Xm , E Xm , and n Xm /n X are determined from a 1D radiationhydrodynamic simulation including non-equilibrium ionisation computed with the Radyn code 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 Xm 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 Xm and gas temperature for Ca ii is shown in the top panel of Fig. 13. Radiative transfer effects make that L Xm 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 Xm 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 Ca ii with the vertical column mass as depth parameter. Again there is considerable spread in the JPDF.
A&A 539, A39 (2012) sses in lines and phere as function e curves give the cording to Eq. (4) hydrogen particle solid), Lyman-β continuum (blue) ts than hydrogen lation 2 . We exsimulation and oid most points peratures above lose to the total er spread, both scape probabil-)) start to break the Lyman tranwards in the athave an escape ith low escape as the adopted to the radiative hydrogen atom normalization radiative losses. Fig. 4. Probability density function of the radiative losses in lines of Ca II as function of temperature in the Bifrost test atmosphere for points above 1.3 Mm height. The curves give the adopted fit (red) and the total collisional excitation according to Eq. (4) (blue).

Fig. 5. Probability density function of the radiative losses in lines of
Mg II as function of temperature in the Bifrost test atmosphere for points above 1.3 Mm height. The curves give the adopted fit (red) and the total collisional excitation according to Eq. (4) (blue).
Lyman-β and the Lyman-continuum are the most important ones). At the low temperature end, the Lyman-continuum starts to contribute and H-α eventually dominates below 7 kK. The hydrogen optically thin radiation dominates over contributions from other elements below 32 kK. Figures 4, 5 show the probability density function of the total radiative losses per singly ionized calcium and magnesium atom, respectively, as function of temperature for the Bifrost simulation. This cooling was computed in detail using the radiative transfer code Multi3d (Leenaarts & Carlsson 2009) assuming statistical equilibrium and plane-parallel geometry with each column in the snapshot treated as an independent 1D atmosphere.
The probability density function is close to the total collisional excitation curve for both Ca II and Mg II. In the case of Ca II this may be a bit surprising at first glance. The Ca II 3d 2 D doublet is metastable and we do not include transitions from it to the ground state. This implies Eq. (2) . 6. Probability density function of the empirical escape probability as function of approximate optical depth at Lyman-α line center and the adopted fit (red). Only points with a cooling above 5 × 10 9 erg s −1 g −1 are included and the startup phase of the simulation (first 1000 s) has been excluded. Fig. 7. Probability density function of the empirical escape probability as function of column mass for Ca II and the adopted fit (red). Only points with a a cooling above 5 × 10 8 erg s −1 g −1 are included.
as function of the chosen depth-scale. We adopt this ratio as the empirical escape probability as function of depth. Figure 6 shows the probability density function of the empirical escape probability as function of approximate optical depth at Lyman-α line center together with the adopted relation. Figures 7 and 8 show the probability density function of the empirical escape probability as function of column mass for Ca II and Mg II together with the adopted relation.
Note that the adopted relations are not the averages of the points shown in Figs. 6−8 but the ratio between the average of the actual cooling and the average of the optically thin cooling such that points with large cooling get larger weight. The adopted fit may therefore look like a poor representation of the PDFs of Figs. 6-8. If these figures had only included the points with the largest cooling, the correspondence would have been more obvious.

Ionization fraction
In Sects. 4.1-4.2 we derived expressions for the radiative losses per atom in the relevant ionization stage per electron. In order to convert these losses to actual losses per volume, the fraction of atoms in the ionization stage under study (H I, Ca II and Mg II) must be known. In general this quantity is a function of A&A 539, A39 (2012) escape probability ted fit (red). Only included. Fig. 10. Probability density function of the fraction of Ca atoms in the form of Ca II as function of temperature in the Bifrost test atmosphere between heights of 0.5 Mm and 2.0 Mm. Curves show the adopted fit to the PDF (red), the coronal approximation (green) and the coronal approximation with a two-level atom (blue). Fig. 13 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. Adapted with permission from Carlsson and Leenaarts (2012). Fig. 14. Cooling as function of column mass from hydrogen transitions at t = 3170 s in the RADYN simulation. Total cooling from the detailed calculation (black solid), from the recipes (black dashed), Lyman-α (blue), H-α (red) and the Lyman continuum (green).  detailed solution with the approximate one and not as pictures of the mean chromospheric energy balance.

Heating from incident radiation field
Part of the optically thin radiative losses from the corona escapes outwards, but an equal amount of energy is emitted towards the that can be pre-computed and tabulated from atomic da emissivity associated with this loss is Because the coronal radiative losses are integrated ov quency one has to choose a single representative opacity pute the absorption. A decent choice is to use the opacit ionization edge of the ground state of neutral helium (d by α). The opacity χ is then given by where N He I is the number density of neutral helium, N H number density of helium, A He is the abundance of heli ative to hydrogen, N H is the total number density of hy particles and N H ρ is a constant dependent on the abundanc quantity N He I /N He can be pre-computed from a calibratio putation as is done for H, Ca and Mg in Sect. 4.3.
The most accurate method for computing the absorp coronal radiation in the chromosphere is to compute the 3 ation field resulting from the coronal emissivity and the tion using the radiative transfer equation: with I the intensity and s the geometrical path length ray. The solution can be computed using standard long characteristics techniques in decomposed domains (Hein et al. 2006;Hayek et al. 2010). Once the intensity is kno heating rate can be directly computed as with Ω the solid angle. The method described above is accurate, but s faster method can be obtained by treating each column MHD simulation as a plane-parallel atmosphere and as that at each location either η or χ is zero. For each col quantities then become only a function of height z and E reduces to with µ the cosine of the angle with respect to the vertical If we furthermore assume the emitting region is loc top of the absorbing region, we can ignore rays that po ward from the interior of the star in the integral of E The intensity that impinges on the top of the absorbing for inward pointing rays is simply the integral of the em A39, page 8 of 10 Fig. 14 Test of the validity of tabuled chromospheric radiative losses. The curves show the horizontally averaged radiative cooling (i.e., −Q/ρ 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. Adapted with permission from Carlsson and Leenaarts (2012). 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.
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 with κ HeI the representative opacity. The angle-averaged radiation field J cor is computed from the transfer equation: with the emissivity given by η = −Q cor /4π. The quantity Q 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 κ HeI 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α 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.

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: 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 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 photo shows the radiative losses per mass unit computed with the method described in Sect. 3; Q chromo displays the chromospheric losses and absorption of coronal radiation per mass unit using the methods from Carlsson and Leenaarts (2012); Qcorona shows the losses per mass unit computed as described in Sect. 4. The panel labeled Q 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.
then the equations for all levels together form a set of coupled first-order differential equations ∂n ∂t = P n, with n the vector of level populations and P the rate coefficient matrix. This system has the solution: where a i are the eigenvectors of the rate matrix with corresponding eigenvalues λ 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 uilibrium is te is small compared to s individually, then the radiative rates will make ation timescale given by ve transitions can either (if NRB > 0) the relaxation timescale (at each ur numerical simulation. y: The value of the hydrom a given time step and xation timescale calculansistent with this state of solving the equations of ed the initial population his atmosphere was then rature by 1% throughout. this perturbed state were tatistical equilibrium giv-1Þ. The time evolution er of protons at a given the full rate equations, ution was cast in a twoxation timescale was calfit to rogen ionization/recomrical simulation, from the ion, is shown in Figure 6 creases outward from the in the midchromosphere he base of the transition rs of the rate matrix, P ij transition rate per atom ify the processes controlrocess can be represented pulations, : level 6 is the continuum. sponding to the ith eigenon the initial conditions. eigenvector is the equiliare all negative, since the eir equilibrium value. The numerically determined relaxation timescale is compared with the timescale from the eigenvalue calculation, which is the inverse of the smallest (in absolute value) nonzero eigenvalue, in Figure 6. We have calculated the eigenvalues using several different assumptions about the rate matrix. Note, first of all, that when all the processes are included in the rate matrix the timescale obtained from the eigenvalues (dotted line) is orders of magnitude smaller than found in the numerical solution (thick solid line). The timescale obtained from the rate matrix with only collisional rates (dot-dashed line) has the same general pattern as the numerical timescale but is generally slightly larger. When radiative rates are added, but with all the Lyman transitions assumed to be in detailed balance, the timescale from the eigenvalues of the rate matrix (thin solid line) reproduce the numerical timescale up to the peak in the midchromosphere. This indicates that in this region the large Lyman radiative rates are in fact changing with time so as to maintain nearly detailed balance as the populations change. The general behavior of the relaxation time is thus controlled by the collisional processes. The timescale increases from the photosphere to the midchromosphere, where the electron density has a minimum approximately as the inverse of the collisional recombination rate (/n À2 e ). From the midchromosphere to the transition region the timescale The electron density is given as the thick dashed curve. The ionization/recombination timescale becomes very long in the chromosphere. Eigenvalues calculated from a rate matrix with all the radiative and collisional rates give a timescale several orders of magnitude too short. Eigenvalues calculated from a rate matrix with all Lyman radiative transitions in detailed balance except for Ly, which is included using its net radiative bracket, give a timescale that closely matches the numerical result.

CARLSSON & STEIN
Vol. 572 it will evolve towards equilibrium on a characteristic timescale given by 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 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.
A similar calculation for helium was done for helium by Golding et al. (2014), finding timescales of the order 10 2 s to 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 and energy conservation 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 solved without having to solve the transfer equation. A limitation of his method is that the Lyman-α line was set in detailed radiative balance. This limitation makes the solution less accurate in the very upper chromosphere and transition region. of about 2500 K. These low values result from the reverse process: over-ionization compared with LTE. More energy remains stored as ionization energy, leaving less kinetic energy for the gas particles.

Limitations of the simulation
Our implementation of non-equilibrium hydrogen ionization has various limitations.
First, the assumption that all Lyman transitions are in detailed balance is justified up to the transition region (Sollum 1999). However, the transition region is optically thin in most Lyman features, requiring detailed radiative transfer modeling to evaluate their influence on the hydrogen populations.
Second, the multi-group radiative transfer within the simulation, which sets the radiative cooling and heating, employs LTE ionization. For given internal energy and mass density, the radiative transfer uses the group-mean opacity, scattering probability and Planck function based on the corresponding LTE (or coronal equilibrium) temperature and electron density. The radiative cooling in the chromosphere and transition region, where deviations from equilibrium are largest, is thus inconsistent with the non-equilibrium temperature and electron density as computed in the simulation.
Third, the cool parts of the simulation chromosphere often reach the limiting temperature of 2400 K allowed in the simulation. It is not clear how low the actual chromospheric minima may reach because radiative heating in the hydrogen continua and other strong spectral features is not taken into account in the simulation, only their radiative cooling.

Discussion
From the analysis of our simulation we obtain the following picture. The internetwork chromosphere is irregularly pervaded of about 2500 K. These low values result from the reverse process: over-ionization compared with LTE. More energy remains stored as ionization energy, leaving less kinetic energy for the gas particles.

Limitations of the simulation
Our implementation of non-equilibrium hydrogen ionization has various limitations.
First, the assumption that all Lyman transitions are in detailed balance is justified up to the transition region (Sollum 1999). However, the transition region is optically thin in most Lyman features, requiring detailed radiative transfer modeling to evaluate their influence on the hydrogen populations.
Second, the multi-group radiative transfer within the simulation, which sets the radiative cooling and heating, employs LTE ionization. For given internal energy and mass density, the radiative transfer uses the group-mean opacity, scattering probability and Planck function based on the corresponding LTE (or coronal equilibrium) temperature and electron density. The radiative cooling in the chromosphere and transition region, where deviations from equilibrium are largest, is thus inconsistent with the non-equilibrium temperature and electron density as computed in the simulation.
Third, the cool parts of the simulation chromosphere often reach the limiting temperature of 2400 K allowed in the simulation. It is not clear how low the actual chromospheric minima may reach because radiative heating in the hydrogen continua and other strong spectral features is not taken into account in the simulation, only their radiative cooling.

Discussion
From the analysis of our simulation we obtain the following picture. The internetwork chromosphere is irregularly pervaded Fig. 17 Temperature and hydrogen ionisation degree n H i /(n H i + n H 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-minute 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. Adapted with permission from Leenaarts et al. (2007), copyright by ESO.
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. 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 α 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.
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 en 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 approxmiate 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 τ ν in a column is then given by:

Figure 2
Temperature in the RADMHD low chromosphere showing a reverse granulation pattern. Lighter (darker) colors indicate hotter (cooler) temperatures. In the models, this occurs because the radiative cooling diminishes with height, and the p∇ · v work of converging and diverging flows above the intergranular lanes begins to dominate. The horizontal slice spans 21 × 12 Mm 2 at a resolution of 448 × 256.

Fig. 19
Temperature at a constant height in the photosphere of a 3D simulation computed with the RADMHD code. The panel spans a size of 24 × 24 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. Adapted with permission from Abbett and Fisher (2012), copyright by Springer.
The first exponential integral E 1 (x) peaks sharply at x = 0, so that this expression can be approximated as This result is inserted into the equation for the radiative cooling (Eq. (11)) to yield 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: where κ B is the Planck-averaged opacity (see Sect. 3.2), τ B the optical depth computed from this opacity, and σ is the Stefan-Boltzmann constant. With this scheme, computation of the radiative heating only requires a simple 2D table lookup to get κ B , and a column-by-column integration over depth to compute τ 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.

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 fieldJ which appears in the radiative rate coefficients of Eq. (54). Judge (2017) suggests to computeJ from the equation (see Frisch and Frisch 1975;Hummer and Rybicki 1982), where τ 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 computeJ 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 problem in a MURaM test atmosphere is ∼ 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.

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 wellstudied. 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 Λ(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 Sections 5 and 6 are likely the most inaccurate of all method used in the photosphere-chromospherecorona system. A long sequence of approximations is needed to order to make the problem computationally feasible. Testing the accuracy of all approxima- tions 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-α line is implemented in the nonequilibrium 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 for example Bjørgen et al. (2019) and Fig. 20).
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 K-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 4 . This does however not mean that the work described in Sections 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α. -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.