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 multidimensional radiationMHD modeling. In this review I discuss the most commonly used approximations for energy exchange between gas and radiation in the photosphere chromosphere and corona.
['The Sun', 'Magnetohydrodynamics', 'Radiative transfer', 'Astrophysics - Solar and Stellar Astrophysics']
Radiation hydrodynamics in simulations of the solar atmosphere
{'No Header': 'Living Reviews in Solar Physics manuscript No.  \n(will be inserted by the editor)', 'Radiation hydrodynamics in simulations of the solar atmosphere': 'Jorrit Leenaarts  \nReceived: date / Accepted: date  \nAbstract 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.  \nKeywords The Sun · Magnetohydrodynamics · Radiative transfer', 'Contents': 'Institute for Solar Physics, Department of Astronomy, Stockholm University  \nAlbaNova University Centre, SE-106 91 Stockholm, Sweden  \nE-mail: jorrit.leenaarts@astro.su.se', '1 Introduction': "The interaction of matter and radiation is an indispensable ingredient of solar and stellar atmospheric modeling. Around 98% of the energy generated in the solar core is transported outwards first by radiative diffusion, then convection, and ultimately escapes into space in the photosphere of the Sun where the overlying solar material becomes transparent. The remainder of the energy escapes the Sun in the form of neutrinos (Turck-Chi'eze 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.  \nModeling 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.  \nWriting general equations that describe the interaction of radiation and matter is not so difficult. Because of limitations in computing time these equations have so far been solved rather completely in one-dimensional geometry only. The Radyn code (Carlsson and Stein 1992, 2002), is perhaps the most well-known example, but other codes are used too (e.g., Flarix, see Kaˇsparov'a et al. 2009).  \nHowever, 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.  \nSince the pioneering simulations of Nordlund (1982), a plethora of codes that aim to model solar and/or stellar atmospheres in 3D have been developed. The ones that I am aware of are: Stagger (e.g., Stein and Nordlund  \n1998, but many different versions of this code exist), MURaM (Vogler 2004; Rempel 2017), Bifrost (Gudiksen et al. 2011), COBOLD (Freytag et al. 2012), Stellarbox (Wray et al. 2015), RAMENS (Iijima and Yokoyama 2015; Iijima 2016), MANCHA3D (Khomenko et al. 2017, 2018), ANTARES (Leitner et al. 2017), and RADMHD (Abbett and Fisher 2012). The latter is the only code that does not use the radiation treatment developed by Nordlund (1982), and instead uses a simpler but much faster method.  \nRadiation-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).  \nIn 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.  \nThe 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.  \nIn the corona the physics of radiative energy losses and gains becomes again somewhat less complex. Most coronal structures are optically thin for all frequencies except in the radio regime. The radio regime lies in the far tail of the Planck function and does not contribute significantly to the radiative losses. For all other frequencies it can be assumed that photon absorption does not take place and that all photons emitted by the gas escape from the corona. It is discussed in Sect. 4.", '2 Fundamentals': 'Excellent books that review the fundamentals of radiation hydrodynamics and modeling of stellar atmospheres are Mihalas and Mihalas (1984) and Hubeny and Mihalas (2014). The first book focusses on fundamental theory, while the second one discusses modeling of 1D static and moving atmospheres. Below I briefly touch upon some general aspects that are relevant for the methods discussed in Sect. 3 - 7.', '2.1 The MHD equations including radiation': "The dominant paradigm for multidimensional modeling of the solar atmosphere has been the magnetohydrodynamics (MHD) approximation. The MHD equations for the density, momentum and internal energy including radiation terms can be written as  \n∂ρ ∂t = -∇ · ( ρ v ) , (1)  \n∂ p ∂t = -∇ · ( v ⊗ p -τ ) -∇ P + J × B + ρ g -∇ P rad , (2)  \n∂e ∂t = -∇ · ( e v ) -P ∇ · v + Q + Q rad , (3)  \nwith ρ 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.  \nThe 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).  \nSeveral radiation-MHD codes have been extended to include some effects beyond MHD, through including the ambipolar diffusion term, the Hall term, and/or Biermann's battery term, to the induction equation (Mart'ınez-Sykora et al. 2012; Cheung and Cameron 2012; Khomenko et al. 2017, 2018) These inclusions retain the single fluid description. This greatly simplifies the treatment of the radiative term Q 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.  \nEfforts 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  \nTable 1 Comparison of energy densities and radiative flux  \ncan 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.", '2.2 Energy density of radiation and matter': 'It is useful to compare the energy density and flux of the radiation field to the energy density of the gas in the solar atmosphere. To get an estimate we assume that the radiation field is isotropic and given by the Planck function at T rad = 5777 K. Then the energy density is  \nE rad = 4 σ c T 4 rad = 0 . 84 J m -3 , (4)  \nwith σ 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  \nF rad = σT 4 rad = 6 . 3 × 10 7 Wm -2 , (5)  \nTable 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.', '2.3 Radiation pressure and force': 'The force exerted by the radiation pressure tensor is typically ignored under normal solar conditions because it is small compared to other forces. To illustrate this one can compare the radiation pressure of isotropic black body radiation to the gas pressure in the photosphere. This pressure is  \n4 σ 3 c T 4 = 0 . 27 Pa (6)  \nat 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:  \na rad ≈ κF c = 3 × 10 -3 ms -2 , (7)  \nwith κ 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 ms -2 .', '2.4 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  \nQ rad = -∫ ∞ 0 ∇· F ν d ν = -∇· F , (8)  \nwith 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, radiationMHDsimulations of the solar atmosphere aim to reproduce the correct detailed behaviour of the gas only. The computation of the radiation field only has to reproduce the correct heating and cooling, but does not need to reproduce the correct spectral energy distribution. This allows for large simplifications in treating the radiation without sacrificing too much accuracy in the value of the total flux divergence.', '2.5 Explicit expression of the radiative flux divergence': 'The monochromatic flux is defined in terms of the intensity I ν as  \nF ν = ∮ n I ν d Ω, (9)  \nwith 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:  \n∇· F = ∫ ∞ 0 ∮ ( j ν -κ ν ρI ν ) d Ω d ν (10) = ∫ ∞ 0 ∮ κ ν ρ ( S ν -I ν ) d Ω d ν,  \nwith κ ν 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  \nwhere  \n∇· F = ∫ ∞ 0 4 πκ ν ρ ( S ν -J ν ) d ν, (11)  \nJ ν = 1 4 π ∮ I ν d Ω (12)  \nis the angle-averaged intensity.  \nThis 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 √ 2 kT/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.', '2.6 Light travel time and hydrodynamical timescales': "Typical bulk flow velocities in the solar atmosphere range from 1 km s -1 in convective upflows to > 300 kms -1 in chromospheric evaporation following solar flares (Graham and Cauzzi 2015). Alfv'en velocities up to 2 . 2 × 10 3 kms -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.  \nThe 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.  \nBecause 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  \n1 c ∂I ν ∂t + n · ∇ I ν = j ν -α ν I ν (13)  \nto  \nn · ∇ I ν = j ν -α ν I ν , (14)  \nwith 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.", '2.7 Diffusion approximation': 'At very large optical depth the following approximation holds for J ν :  \nJ ν ≈ S ν + 1 3 d 2 S ν d τ 2 ν , (15)  \nwhere the second derivative should be taken along the direction of the gradient of S ν 1 . For the flux the following approximation holds:  \nF ν ≈ -4 π 3 κ ν ρ ∇ S ν . (16)  \nAt 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 ν :  \nF ≈ -16 3 σT 3 κ R ρ ∇ T, (17)  \nwith the Rosseland opacity defined as  \n1 κ R = (∫ ∞ 0 1 κ ν d B ν d T d ν )/(∫ ∞ 0 d B ν d T d ν ) . (18)  \nMethods 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):  \n∫ ∞ 0 4 πκρ ( S ν -J ν ) d ν ≈ N ∑ i =1 4 πκ i ρ ( S i -J i ) , (19)  \nwith κ i , S i , and J i the opacity, source function and radiation field in each group, as defined below.  \nThe 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  \nbe Λ κ [ . . . ]. For two frequencies ν 1 and ν 2 with identical opacities everywhere in the atmosphere κ 1 = κ 2 = κ one can thus write  \nJ ν 1 + J ν 2 = Λ ν 1 [ S ν 1 ] + Λ ν 2 [ S ν 2 ] = Λ κ [ S ν 1 + S ν 2 ] (20)  \nIn 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:  \n∫ ∞ 0 κ ( S ν -J ν ) d ν glyph[similarequal] ∑ j w j κ j ( S j -J j ) (21) = N ∑ i =1 ∑ j ( i ) w j κ j ( S j -J j ) = N ∑ i =1 ∑ j ( i ) w j κ j ( S j -Λ j [ S j ]) (22) ≈ N ∑ i =1 κ i   ∑ j ( i ) w j S j -Λ i [ ∑ j ( i ) w j S j ]   ≡ N ∑ i =1 κ i ( S i -Λ i [ S i ]) (23) ≡ N ∑ i =1 κ i ( S i -J i ) . (24)  \nHere 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  \nS i = ∑ j ( i ) w j S j (25)  \nWhat 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.  \nSection 3.1 discusses how to sort frequencies into opacity groups. In Sect. 3.2 it is assumed that the source function is given by the Planck function. This assumption is relaxed to allow for coherent scattering in Sect. 3.3.', '3.1 Sorting frequencies into groups': "Grouping frequencies and defining a bin opacity are necessarily somewhat crude. The density and thus opacities drop roughly exponentially with height, and opacities at a given point in the atmosphere vary over many orders of  \nFig. 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).  \n<!-- image -->  \nmagnitude 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.  \nThe 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  \nτ k -1 R < τ R ( τ ν = 1) ≤ τ k R . (26)  \nA common choice is to put the border values τ k R equidistantly spaced in log τ R . The method is illustrated in Fig. 1.  \nex  \nj  \n(  \ni  \n,  \nl  \n) runs over those  \nB  \ne  \nff  \nect of the sorting pro-  \nstrated in Figs. 2 and 3.  \ny from the ATLAS9 stel-  \nT  \n=  \n4470 K and log  \np  \n=  \noptical depth at 500 nm  \n). After applying the tau-  \nns the spectrum is a step  \nvalues (Fig. 3).  \npproach the exact opac-  \nsmall ODF intervals, the  \ned to converge to the ex-  \nthe detailed reference  \ng process) in the limit  \n. This is due to the fact  \nbetween bin-levels  \nτ  \nl  \nref  \nwill comprise frequen-  \nwith di  \nff  \nerent height  \n. It should be noted that  \n∆  \nν  \ni  \ndenotes the average  \ni+1  \nFig. 2.  \nMonochromatic opacity as a function of wavelength for  \nT  \n=  \nτ  \nref  \nFig. 1.  \nSchematic illustration of the  \nτ  \n-sorting procedure using ODFs.  \nWithin each ODF interval  \n∆  \nν  \nquency bins  \nΩ  \nl  \n,  \ndepending on the height where  \nτ  \nal.: Non-grey RT in numerical simulations of the solar photosphere  \n∆ν  \nν  \ning procedure using ODFs.  \ns  \n∆  \nν  \nwhere  \nτ  \nij  \nare sorted to fre-  \nij  \n=  \n1 is reached  \ne major requirement for  \n. For greater depths, the  \ne inappropriate, but this  \nn transfer becomes in-  \nay.  \nrealized in a convenient  \nep  \nj  \nwithin ODF inter-  \ngives  \nτ  \nij  \nas a function of  \nFig. 1). The ODF steps  \ngrated quantities are ob-  \nsteps in a given bin, viz.  \n(16)  \n(17)  \nwithin one bin are similar, which is the major requirement for  \nval  \nBl  \nand  \nare used to model the line opacities.  \nFig. 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. κ P , l = 1 Bl ∑ i ∆ ν i B ∆ ν i ∑ j ( i , l ) w j ( i , l ) κ i j ( i , l ) , (17) where, for a given bin index l , the index j ( i , l ) runs over those steps in ∆ ν i ,whichareelementsof Ω l . B ∆ ν i denotes the average of B over the ODF interval ∆ ν i . The e ff ect 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 = Fig. 3. Opacity as a function of wavelength for a τ -sorted multigroup model with five opacity bins, for the same thermodynamical conditions as in Fig. 2. using ODFs as basis for sorting frequencies introduces an additional error since the rearrangement of spectral lines inherent in the ODF concept a ff ects the way frequencies are classified. However, since the ODF solution approximates the exact soluFig. 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 Vogler et al. (2004), copyright by ESO.  \n<!-- image -->  \n4  \n.  \n12 (corresponding to a continuum optical depth at 500 nm  \nlog  \nτ  \n500  \n=  \n2 in the solar atmosphere). After applying the tau-  \nsorting procedure with five opacity bins the spectrum is a step  \n-  \ntion well, the error incurred this way is small compared to the  \nconsequence of the binning procedure (see Sect. 3.5).  \nThe choice of a reference atmosphere is another possible  \nFig. 3. Opacity as a function of wavelength for a τ -sorted multigroup model with five opacity bins, for the same thermodynamical conditions as in Fig. 2. 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 di ff erent parts of the spectrum with di ff erent height profiles of the corresponding opacities. It should be noted that 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 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).  \nusing ODFs as basis for sorting frequencies introduces an additional error since the rearrangement of spectral lines inherent in the ODF concept a ff ects 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 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.  \nconsequence of the binning procedure (see Sect. 3.5).  \nThe choice of a reference atmosphere is another possible  \nsource of error connected with the  \nτ  \n-sorting scheme. It is plau-", 'sible that the opacity binning method with τ -sorting shows good results in calculations of static 1D model atmospheres if 3.2 Multi group radiative transfer with LTE source function': "the 'exact' solution (i.e., the 1D atmosphere resulting from the  \nODF 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 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.  \ni  \n, the ODF steps  \n∆  \nν  \nij  \nare sorted to fre-  \nij  \n=  \n1 is reached  \nThe source function in group i is then given by  \nS i = ∑ j ( i ) w j B j . (27)  \nThis expression only depends on temperature and can thus easily be precomputed and stored in a 1D lookup table.  \nOnce 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:  \nThis 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:  \n1 κ R i =   ∑ j ( i ) w j 1 κ j d B j d T   /   ∑ j ( i ) w j d B j d T   . (28)  \nκ B i = ∑ j ( i ) w j κ j B j ∑ j ( i ) w j B j . (29)  \nSomewhere 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  \nκ i = Wκ B i +(1 -W ) κ R i , (30)  \nW = e -τ i /τ 0 . (31)  \nwhere τ 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; Vogler et al. 2004).  \nAnother option is to base the transition on the local mean free path:  \nW = e -l i / ( κ i ρ ) , (32)  \nwhere 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.  \nFor 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.  \nExtensive 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 Vogler et al. (2004).  \n3.3 Multi group radiative transfer with non-LTE source function  \nThe 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  \nglyph[epsilon1] = C ul A ul + C ul + B ν B ul < 10 -4 , (33)  \nwith 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.  \nScattering has a strong damping effect on the amplitude of ∇· F ν /ρ . In LTE it is given by  \n∇· F LTE ν ρ = 4 πκ ν ( B ν -J ν ) . (34)  \nIf one assumes the presence of a coherently scattering line in a two-level atom at frequency ν , then the source function becomes  \nS ν = (1 -glyph[epsilon1] ) j ν + glyph[epsilon1]B ν . (35)  \nIn that case the flux divergence per mass unit is  \n∇· F NLTE ν ρ = 4 πglyph[epsilon1]κ ν ( B ν -J ν ) , (36)  \nso 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 -glyph[epsilon1] ν ) J ν + glyph[epsilon1] ν B ν and following a reasoning similar to the one given in Sect. 3.2, he arrives at an expression for the group source function:  \nS i = glyph[epsilon1] i J i + t i . (37)  \nHere glyph[epsilon1] 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, glyph[epsilon1] i , t i , and κ i have different expressions for the diffusion regime and in the free streaming  \nFig. 3 Illustration of short and long characteristics used for solving the transfer equation in 2D decomposed domains. Grey dots indicate grid points where the intensity should be computed, with grey lines connecting those grid points. Blue lines indicate subdomain boundaries, where it is assumed that the horizontal boundaries are periodic. Red arrows indicate two examples of long characteristics, while the orange arrows indicate short characteristics. Note that the long characteristics cross multiple subdomain boundaries and wrap around the periodic horizontal boundary Each subdomain contains a piece of the long characteristic of variable length.  \n<!-- image -->  \nregime 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.  \nAnother 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).", '3.4 Solving the transfer equation': "Most modern radiation-MHD codes are parallelized to make use of large supercomputers with distributed memory. Typically, the simulated domain is represented on a 3D Cartesian grid. The domain is split into smaller subdomains and each CPU handles the required computations on its own subdomain, communicating information to other subdomains as needed. This architecture scales well for the MHD equations: they are local and require communication with neighbouring subdomains only.  \nRadiation 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 ( ˇ Stˇep'an 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; Vogler et al. 2005; Heinemann et al. 2006; Hayek et al. 2010; Ibgui et al. 2013; ˇ Stˇep'an and Trujillo Bueno 2013). I will only briefly touch upon some aspects.  \nSolving 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.  \nShort 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:  \nI ( τ ) = I ( τ = 0) e -τ + ∫ τ 0 S ( t ) e t -τ d t, (38)  \nwhere 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.  \nThe 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  \ngeometry 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 ˇ Stˇep'an 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.  \nBruls 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.  \nThe 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).  \nLong 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.  \nBoth 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.  \n3.5 Computation of the heating rate from the intensity and source function  \nAnalytically, the two expressions for the heating rate in a bin Q i = -∇· F i = 4 πκ i ρ i ( J i -S i ) are equivalent. In case of actual numerical computation this is no longer the case.  \nIn 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  \ncompute 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.  \nIn 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 .  \nFollowing 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.  \nAn 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  \nd K I d τ = d S d τ -K I , (39)  \nwhich in its integral form is given by  \nK I ( τ ) = K I ( τ 0 )e τ -τ 0 + ∫ τ τ 0 e τ ' -τ d S d τ ' d τ ' . (40)  \nThe total heating rate is then computed from  \nQ = ∮ K I d Ω (41)  \nThis 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 .  \nNordlund (1982) implemented a similar method as Heinemann et al. (2006), but based on the second-order form of the transfer equation:  \nd 2 P d τ 2 = P -S. (42)  \nFig. 4 Summary table of the three major binary choices in multi-group radiative transfer from computing radiative losses and gains in the photosphere.  \n<!-- image -->  \nwhere  \nP = 1 2 ( I ( Ω ) + I ( -Ω )) , (43)  \nthe average of an ingoing and an outgoing ray (e.g., Hubeny and Mihalas 2014, p. 387). Defining K P = P -S one arrives at  \nd 2 K P d τ 2 = K P -d 2 S d τ 2 . (44)  \nThis 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).  \n3.6 Summary and examples of photospheric radiative transfer  \nFigure 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.  \nThe simulation by Nordlund (1982) is of Type 5. MURaM (Vogler 2004; Rempel 2017), RAMENS (Iijima and Yokoyama 2015; Iijima 2016), and MANCHA3D (Khomenko et al. 2017, 2018) are Type 2 codes. Stellarbox (Wray et al. 2015) is Type 6. COBOLD (Freytag et al. 2012) has both Type 2 and Type 6 options, but for stellar atmosphere simulations only Type 2 is supported.  \nFig. 5 Example of the vertically emergent intensity computed with a 4-group non-LTE scheme. The flux divergence along the horizontal white line is shown in Fig. 6. The model was computed with the Bifrost code (Gudiksen et al. 2011).  \n<!-- image -->  \nThe 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.  \nFigures 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).  \nFigure 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  \nFig. 6 Atmospheric structure and flux divergence per mas unit Q i /ρ = -∇ · F i /ρ = 4 πκ i ( J i -S i ) in a vertical slice along the white line in Fig. 5. The top row shows the temperature and density, the two bottom rows show the flux divergence per mass unit in the four opacity groups. Brown is cooling, blue is heating, and the brightness scale for the flux divergence panels is clipped at 20% of the maximum of the absolute value to enhance contrast. The maximum and minimum values are indicated in each panel. The black curves indicate the τ 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.  \n<!-- image -->  \nimages 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.  \nFigure 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.  \nSKARTLIEN  \nF IG . 3.È 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. L ower panel: Horizontally averaged amplitudes of radiative heating (average absolute value) relative to the amplitude of the exact total heating. Thick lines are the exact amplitude ratios, while the thin lines are ratios from the approximate group solutions. Fig. 7 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.  \n<!-- image -->  \napproximate solution. The horizontal averages (  \nthin lines  \n)  \ncoincide very well with the horizontal averages of the exact solutions ( thick lines ). The upper panels show cooling below 0.0 Mm in all groups, and heating in groups 2, 3, and 4 immediately above 0.0 Mm. The lower panel shows the amplitude of the approximate and exact Ñux divergence (horizontal average of the absoFigure 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  \nlute value) relative to the amplitude of the  \ntotal exact  \nÑux  \ndivergence,  \ni.e.,  \nand  \n*  \no  \nT  \n/  \nS  \no  \n/  \no  \nT  \n/  \nS  \no  \n/  \n(brackets denote horizontal average). We see that the exact  \nS  \no  \n/  \nrelative amplitudes (  \nthick lines  \n) coincides very well with the  \napproximate values (  \nthin lines  \n).  \nA sample of the spatial structure is displayed in Figure 4,  \nwhere we have shown the exact heating per mass unit in all  \ngroups as gray-scale images in a vertical slice. Black con-  \ntours mark locations of zero heating, and lighter shades of  \ngray means positive heating. Gas in layers immediately  \nabove the cooling layer is heated in all groups in expanding  \nÑow above granules. Granules are seen as curved horizontal  \nstructures. As up-Ñowing gas expands and cools, the tem-  \nperature falls below the radiation temperature, and the gas  \nis radiatively heated. Note also heating of the cool region  \ni  \ntot  \no  \nT  \nS  \no  \n/  \ni  \ntot  \no  \nT  \nIG  \nF  \n.  \n4.ÈExact group radiative heating per mass  \nthrough the simulation. Full drawn black contours  \nheating. Heating is found at lighter shades of gray, and  \nshades. Dash-dotted black contours show the zero level  \ngroup heating and coincide well with the zero level of  \nWhite curves are the horizontal averages of the radiative  \nunit (normalized to  \nÐ  \nt the plotting window). Black  \nzero level for these curves, and positive values are to the  \nL ower panel:  \nTemperature in the same vertical slice.  \nthe granular layer at the height 0.0 Mm, and heating  \nall groups, and also heating of the cool region  \nchromosphere.  \nbelow 3000 K in the chromosphere as  \nfrom below is converted to thermal energy via  \nThe dash-dotted black contours show the  \nthe approximate heating and coincide well  \ngroups 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.  \nIn 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:  \n∇· F i ρ = 4 πglyph[epsilon1] i κ i ( S i -J i ) , (45)  \nwhere glyph[epsilon1] i = 1 in LTE, and glyph[epsilon1] i ≤ 1 in non-LTE. In the latter case glyph[epsilon1] 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 glyph[epsilon1] i = 1 in LTE.  \nThe amplitude of the cooling and heating in groups 2 - 4 between z = 0 Mm and z = 0 . 3 Mm is however much larger in non-LTE than in LTE. Skartlien (2000) speculates that this is caused by the smoothing effect that scattering has on J i , but a thorough investigation of this effect has never been done.", '4 Radiative losses in the transition region and corona': "The corona is optically thin at all wavelengths except for the radio regime. Its radiative losses are dominated by a myriad of EUV lines from highly ionised stages of many different elements (e.g., Curdt et al. 2004; Woods et al. 2012).  \nIn 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.  \nNon-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ˇc'akov'a 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.  \nA fully general treatment of TR and coronal radiative losses in a radiationMHD 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.  \nNeglecting radiative non-local transfer effects through assuming optically thin radiative transfer (i.e., assuming J ν = 0 in the rate equations) alleviates  \nJ  \nand source  \nbut the scaling  \nfrom the new  \nprevious LTE  \ncaused mainly by  \nfrom di†ering  \nshow the  \nHorizontal aver-  \ncorrespond to  \nequal to  \nthe scaling by  \nall group  \nline scat-  \nsource func-  \nsnapshot that  \nthe same set  \nresults.  \nVol. 536  \nF IG . 7.ÈHorizontal averages of radiative heating per mass unit. The units are arbitrary, and the scaling is the same for all Ð gures. The amplitude of the Ñux divergence for the previous LTE method (old) is much larger above 0.3 Mm than for the current method. This e†ect is caused mainly by the larger di†erences 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. L ower panel: Sum of all groups. The main di†erence above 1.0 Mm comes from the contribution in group 4. 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 n ew). 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.  \n<!-- image -->  \ntering would imply that  \nThe values of the mean  \nJ  \n*  \n\\  \nS  \n*.  \nintensities are also di†ering, mainly because of di†erences in  \ni  \ni  \nthe source function, and partially because of the di†erent  \nheights for which the optical depths are unity.  \nDiamonds and triangles show the heights for where the  \nhorizontally averaged optical depth  \nis unity, for previous  \nand current methods, respectively. The mean intensities in  \nq  \ni  \noptically thin regions are approximately  \nS  \nJ  \n(  \nq  \ni  \n*  \nT  \n^  \n1  \nS  \nS  \n*  \nT  \nThe group mean opacities that determines the  \n2  \ni  \n\\  \n1).  \noptical depths are seen in the bottom panel of Figure 6. The  \ni  \nnew opacities are higher than the old opacities above 0.3  \nalready 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.  \nThe problem reduces then to solving the rate equations  \nglyph[negationslash]  \nglyph[negationslash]  \n∂n i ∂t + ∇· ( n i v ) = n l ∑ j,j = i n j P ji -n i n l ∑ j,j = i P ij , (46)  \nwhere 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  \nQ ij = hν ij A ji n j , (47)  \nand 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.  \nThe 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.  \nInstead, 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:  \nQ ij = hν ij A ji n e n j,m n m n m n X n X n H n e n H , (48)  \n≡ G ( T, n e ) n e n H . (49)  \nHere 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  \nE. Landi & M. Landini: Radiative losses of optically thin coronal plasmas  \nFig. 5. Percentual differences between total emissivity curves calculated assuming different values of the electron density. Full line : 10 8 vs. 10 10 cm -3 ; Dashed line : 10 8 vs. 10 12 cm -3 ; Dash-dotted line : 10 8 Fig. 6. Fig. 9 Sensitivity of the coronal loss function Λ ( T, n e ) to the electron density. Full line: 10 14 vs. 10 16 m -3 ; dashed line: 10 14 vs. 10 18 m -3 ; dash-dotted line: 10 14 vs. 10 20 m -3 . Adapted with permission from Landi and Landini (1999), copyright by ESO.  \n<!-- image -->  \nvs. 10  \n14  \ncm  \n-  \n3  \n.  \nequilibrium, while for lower densities collisional de-excitation in G ( T, n e ) through collisional deexcitation and density-dependent dielectronic recombination (Summers 1972, 1974).  \nPercentual difference between the total emissi  \ntained with the old and new version of the Arcetri  \nadopted electron density is 10  \n10  \ncm  \n-  \n3  \n.  \nAs big improvements have been done in the  \nsion of the Code versus the older version  \n& Monsignori Fossi 1990, we have performed a  \nbetween the present results and those obtained  \nversion of the Arcetri Code. The adopted element  \nare from Allen 1973. There are three main dif  \nthe two versions of the Arcetri Code:  \n(a)  \nthe old  \nall  \nline intensities using the  \nCoronal  \n(b)  \nthe collision rates were calculated  \n(c)  \nradiative data came from different  \nbecomes negligible compared to radiative decay and the Coronal ModelApproximation (yielding density insensitive line Contribution Functions ) may be adopted. Fig. 5 displays the percentual difference Perc.Diff ( i ) = η ( T,N e = 10 8 ) -η ( T,N i e ) η ( T,N e = 10 8 ) (6) (with N i e = 10 10 , 10 12 and 10 14 cm -3 ) between total emissivity calculated imation , 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 ).  \nthan in the present version of the Code.  \nThus, the present comparison allows to check  \nof different assumptions in level population  \nresulting total plasma emissivity.  \nFig. 6 displays the percentual difference  \ncurves 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 hightors and A larger uncertainty can be introduced by inaccuracies of the elemental abundances. The coronal loss function is dominated by losses from C, Si, O, and Fe, and the abundances of these elements have a large influence. The abundance of C and O is debated after 3D non-LTE computations (Asplund et al. 2004, 2005) gave a different result than calculations using 1D models Grevesse and Sauval (1998).  \nest 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 Perc.diff 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  \nreason there are small differences between computations carried  \nout assuming different density values.  \n3.2. Effect of different datasets and approximations  \nin level population computation  \nLevel populations are strongly sensitive to any change or prob-  \nlem in the atomic parameters, collision strengths and transition  \nprobabilities as well as in the approximation adopted for their  \ncalculations, and this affects line radiation. It is therefore im-  \nportant to check the effects of different transition probabilities  \ndatasets on the resulting total emissivity curve.  \nη  \n-  \nη  \nold  \nbetween the two versions of the Code as a function  \ntemperature. It is possible to see that rather high dif  \nto 60%) are found at transition region temperatures,  \ndiscrepancies occur at coronal temperatures. In the  \ntion of Fig. 6 the older version of the Arcetri  \ntotal emissivity than the more recent version at  \ntemperatures. This is due to the presence of few  \ntransitions from O  \niv  \n, O  \nv  \n, C  \niv  \nwhose emissi  \ndifferent values in the two versions of the Code;  \nis due both to the use of different datasets and to  \napproximations used in level population  \nan overestimation of line emissivity for these  \nold version of the Code. The negative section of  \nis due to the much larger number of lines  \nversion.  \n=  \nη  \nold  \nnew  \n926  \nK. P. Dere et al.: CHIANTI - an atomic database for emission lines. IX.  \nFig. 5. Radiative loss rate for coronal abundances ( upper curve ) and photospheric abundances ( lower curve ). 1 . 6 Fig. 10 Radiative loss functions Λ computed assuming n e = 10 15 m -3 as function of temperature. Upper curve: coronal abundances from Feldman et al. (1992). Lower curve: photospheric abundances from Grevesse and Sauval (1998). Adapted with permission from Dere et al. (2009), copyright by ESO.  \n<!-- image -->  \n∫  \nemission measure (  \nn  \ne  \nN  \nHd  \nV  \n) has been calcu  \nization equilibria discussed in Sect. 4. The lo  \nprocesses of bremsstrahlung, radiative  \nation and two-photon radiation. An electron d  \nwas used and two sets of elemental abundan  \nabundances of Feldman et al. (1992)andthe  \ndances of Grevesse & Sauval (1998). The val  \nloss rate are shown in Fig. 5 and in Table 1.  \nThe most recent calculation of radiati  \nthose of Landi & Landini (1999); Colgan et  \n& Landini (1999) used the Arcetri spectr  \nCHIANTI 2 database. Their Fig. 1 shows cal  \nsity of 10  \n10  \ncm  \n-  \n3  \nfor the coronal abundances  \net al. 1992)andtheionizationequilibriaofAr  \n(1985) and Arnaud & Raymond (1992). A c  \nthat the calculations of Landi & Landini (  \nto our current values. The major di  \nff  \nerences  \nwhere the current loss rate is a factor of  \n×  \n10  \n4  \nK where the current loss rate is a  \nthan the values of Landi & Landini (1999).  \nhigher temperature is primarily due to the f  \ndid not include transitions for the hydrogen-li  \n46%lower in the case of oxygen, with the changes largely resulting from the application of 3D model atmospheres to the interions. The di to di ff of Cranmer and Winebarger 2019, and references therein). Abundances might thus well be the dominant source of uncertainty in the coronal loss function.  \nff  \nerences at the lower  \nerences in the way the Arcetri code and  \nlate the radiative losses due to the continuum.  \nThere are significant di  \nff  \nerences between  \nColgan et al. (2008) and those of both Landi  \nand our current values. Detailed comparisons  \n(2008) are not possible as none of the atomi  \nauthors' radiative loss calculation were publi  \npretation 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 preThe CHIANTI atomic database 3 (Dere et al. 1997, 2019) provides an extensive compilation of critically assessed atomic data. In addition to the data it delivers software packages in both the IDL and Python languages for using this data and easy generation of loss functions using a variety of abundances and other input options.  \nnote that for most ions radiative losses are  \nsmall number of strong fine structure tran  \nconfigurations, and so modeling these trans  \nviously been in excellent agreement with the precise helioseisComputation of Q in the corona is then straightforward:  \nmic measurements of parameters such as the sound speed profile. Oxygen is the dominant contributor to opacity in the region Q = -Λ ( T ) n e n H , (50)  \natomic data is crucial to deriving accurate  \nmates. Colgan et al. (2008)treatedentireconfi  \nlevels in their calculations with no fine struc  \nno attempt was made by the authors to asse  \ntheir approach for the dominant low-lying con  \nCHIANTI fine structure levels are used for  \njust below the convection zone, and the revised abundance significantly a ff ected 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 di ff erent measurements and re-analyses of oxygen lines and photospheric models (e.g., Centeno & Socaswhere Λ ( 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.  \nmost accurate atomic data for transitions be  \nare chosen following assessments of the best a  \nliterature. We are thus confident in the accurac  \nradiative losses for the ion models used in the  \n24. New IDL procedures  \nThe parameters that are necessary to calculat  \nsections and rate coe  \nffi  \ncients, r  \nand ionization equilibria are now included  \ndatabase. Accordingly, new Interactive Data  \nfunctions  \nioniz\\_cross  \n,  \nioniz\\_rate  \n,  \nmake\\_ioneq\\_all  \nare also provided.  \nioniz\\_cross  \nreturns the ionization cross s  \nfied ion and electron energy in eV.  \nioniz\\_rate  \nreturns the ionization rate coe  \nfied ion and temperature in K.  \nrecomb\\_rate  \nreturns the recombination r  \nspecified ion and temperature in K  \nmake\\_ioneq\\_all  \ncreates a new ionization e  \nNavarro 2008; Meléndez & Asplund 2008; Ca ff au 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 di ffi cult for a project such as CHIANTI which seeks to provide a single set of reference abundances for users. Our solution is to retain the 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).  \nthe CHIANTI software, and to provide the Grevesse et al. (2007)  \nabundances as an option that can be selected by the user.  \nAlthough the change in the oxygen abundance is the most  \nsignificant in the Grevesse et al. (2007) tabulation, it is also to  \nbe noted that the neon and argon abundances have changed by  \na similar amount. This is on account of the fact that neither of  \nthese elements can be measured in the solar photosphere, and  \nso they have usually been measured relative to oxygen by other  \nmeans (spectroscopy of the solar upper atmosphere, solar ener-  \ngetic particle measurements).  \nCHIANTI users are encouraged to bear in mind the ongo-  \ning debate over solar photospheric abundances and to assess the  \nimpact on their analyses.  \n23. Radiative loss rates  \nFig. 11 Radiative losses in the transition region and corona in a radiation-MHD simulation computed with the MURaM code. Top: temperature; middle: mass density; bottom: -Q , i.e., a positive value means radiative cooling. Also note the strong corrugation of the chromosphere-TR boundary. Figure computed from a simulation by Danilovic et al. (2019).  \n<!-- image -->  \nThe 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.  \nA 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.  \n>  \n00  \nLO  \nLO  \ni-D  \na  \nCO  \nCTi  \n704  \nVERNAZZA, AVRETT, AND LOESER  \nVol. 45  \nFig. 49. - Net radiative cooling rates for Ca n, Mg n, H , H, and other constituents calculated from model C at 50.4 nm, all calculated from model C. In all the calculations we have carried out for various atoms and ions, namely, H, H ', neutral and ionized He, C, Si, Mg, Ca, and O, and neutral Fe, Al, Na, and B, these are the 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 reFig. 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.  \n<!-- image -->  \nonly line and continuum rates that appear to exceed  \n- 4  \n10  \nergs cm  \n-3  \ns '  \n1  \nin absolute value.  \nThe calculated net radiative cooling rates for the Ca n  \ngion are assumed to be zero in each integral. Athay  \n(1976, Table IX-1) estimates that the contributions to  \ndh have the following order of importance: Balmer  \ncontinuum, Ha, Ca n IR triplet, Mg n, Ca n K and H,  \nFe i lines, and Na i D. Athay states, however, that these", 'and Mg n resonance Unes are quite sensitive to the convergence and internal consistency of our numerical 5 Radiative transfer in the chromosphere': "solution. Some of our preliminary solutions gave maxi-  \nmum rates for these lines 2 to 3 times larger than those  \nestimates are not based on detailed computations of the  \nnet energy loss or gain in the bound-free continua. The  \nshown in Figure 49. order we determine is Usted in Table 29. ' American Astronomical Society ¥ Provided by the NASA Astrophysics Data System 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.  \nProper 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 radiationMHD simulations (e.g., Gudiksen et al. 2011). Simplifications that speed up the computation are thus required.  \nCarlsson 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.  \nThe method approximates the radiative loss per volume owing to species X in ionisation state m as  \nQ X m = -L X m E X m ( τ ) n X m n X A X n H ρ n e ρ. (51)  \nHere, L X m is the optically thin radiative loss function per electron and per particle of element X in ionisation stage m , E X m ( τ ) 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.  \nThe quantities L X m , E X m , and n X m /n X are determined from a 1D radiationhydrodynamic simulation including non-equilibrium ionisation computed with the Radyn code (Carlsson and Stein 1992, 2002) for H i . For Ca ii and Mg ii they were determined from a 2D radiation-MHD simulation with Bifrost that provided the atmospheric structure and subsequent statistical equilibrium radiative transfer calculations including PRD using Multi3d (Leenaarts and Carlsson 2009).  \nThe loss function L X m can be computed for each grid cell in the simulation by summing the net downward radiative rates multiplied with the energy difference of the transition, summed over all relevant transitions. The joint probability density function (JPDF) of L X m and gas temperature for Ca ii is shown in the top panel of Fig. 13. Radiative transfer effects make that L X m is no longer a unique function of T , but the red curve indicates an approximate fit. The ionisation degree can be computed from the atomic level populations in the calibrating simulations (see the bottom panel of Fig. 13). Again they are no longer clean functions of temperature, but the spread is minor and a sensible fit as function of temperature can be made. Finally, the empirical escape probability E X m is computed from the total radiative losses in the simulation, the radiative loss function and the ionisation degree. The middle panel of Fig. 13 shows E Ca ii with the vertical column mass as depth parameter. Again there is considerable spread in the JPDF.  \nuchmorelikely  \nP states than to  \ntransition  \nhe ground state  \nmber of radia-  \nstate through a  \n30 A&A539, A39 (2012)  \nhe original col-  \n(4) still holds.  \nlevels give a  \nmes are  \nsmaller  \nstate. In con-  \nnot caused by  \nof the values  \nlowest temper-  \n(not visible in  \ndip downwards.  \nsses in lines and  \ns lead to emit-  \nhere as function  \ncurves give the  \nhere increased  \norption of pho-  \nexcitation pro-  \nto Eq. (4)  \nsible routes in  \nthese processes  \nr equations and  \nempirical es-  \nling as obtained  \nnd the optically  \neally we would  \nction of optical  \ne total cooling  \nl depth scale to  \nA&A539, A39 (2012)  \nlate the escape  \nationale behind  \nnce lines is pro-  \nion stage under  \ncase for both  \non of the atmo-  \nother advantage  \nhydrogen particle  \nsolid), Lyman-  \nβ  \nontinuum (blue)  \nto calculate  \nan MHD simu-  \nLyman lines  \ns than hydrogen  \nregion where  \ned for calculat-  \ngion stays at a  \nlation  \n2  \n.  \nWe ex-  \nsimulation and  \noid most points  \ne, the hydrogen  \nnction. Instead,  \nescape probability  \nperatures above  \nonal to the neu-  \ned fit (red). Only  \nincluded.  \nlose to the total  \n14  \ner spread, both  \ncolumn density  \n4  \n10  \n-  \ncm  \noptical depth at  \ncape probabil-  \n×  \n)) start to break  \nthe Lyman tran-  \nline center and the  \ng  \n-  \n1  \nFig. 4.  \n<!-- image -->  \nFig. 7.  \nFig. 5.  \nMg  \ndepth at Lyman-  \nLyman-  \nand Mg  \nNote that the adopted relations are not the averages of the  \nto contribute and H-  \nα  \neventually dominates below 7 kK. The  \nhydrogen 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 us2 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 Fig. 10. Probability density function of the fraction of Ca atoms in the form of Ca  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).  \nis high to rather  \nattering proba-  \nwards in the at-  \nhave an escape  \nith low escape  \nway on the  \nradiation field.  \nn mass density,  \nas the adopted  \nto the radiative  \nhydrogen atom  \nlculate the aver-  \ner space for the  \ncooling accord-  \nnormalization  \nmore obvious.  \ning the radiative transfer code Multi3d (Leenaarts & Carlsson  \n2009)assumingstatisticalequilibriumandplane-parallelgeom-  \netry with each column in the snapshot treated as an independent  \n1D atmosphere.  \n4.3. Ionizationfraction  \nThe probability density function is close to the total col-  \nlisional excitation curve for both Ca  \nIn Sects. 4.1-4.2 we derived expressions for the radiative losses  \nper atom in the relevant ionization stage per electron. In order  \nand Mg  \n.  \nIn the case  \nto convert these losses to actual losses per volume, the frac-  \nof  \nCa  \nCa  \n3d  \n2  \nthis  \nmay be a bit surprising at first glance. The  \nD doublet is metastable and we do not include tran-  \nsitions from it to the ground state. This implies Eq. (2) does not  \ntion of atoms in the ionization stage under study (H  \n, Ca  \nand  \nFig. 6.  \nbeen excluded.  \nCa  \nFig. 14.  \nCooling as function of column mass from hydrogen transitions  \nat  \nt  \n=  \n3170 s in the RADYN simulation. Total cooling from the detailed  \nFig. 16. As Fig. 15 but for the low-mid chromosphere. detailed solution with the approximate one and not as pictures that at each location either reduces to Q chr( z ) = 2 πχ ( z ) ∫ 1 -1 I ( z ,µ 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).  \n<!-- image -->  \nemissivity associated with this loss is  \nη  \n=  \n-  \nQ  \ncor  \n4  \nπ  \n·  \nBecause the coronal radiative losses are integrated o  \n31 quency one has to choose a single representative opacity  \npute the absorption. A decent choice is to use the  \nionization edge of the ground state of neutral helium (  \nby  \nα  \n). The opacity  \nχ  \nis then given by  \nχ  \n=  \nα  \nN  \nwhere  \nN  \nHeI is the number density of neutral helium,  \nN  \nnumber density of helium,  \nA  \nHe is the abundance of he  \native to hydrogen,  \nN  \nH is the total number density of h  \nparticles and  \nquantity  \nN  \nHeI  \nN  \nH  \nρ  \n/  \nN  \nHe can be pre-computed from a  \nputation as is done for H, Ca and Mg in Sect. 4.3.  \nThe most accurate method for computing the abso  \ncoronal radiation in the chromosphere is to compute the  \nation field resulting from the coronal emissivity and the  \ntion using the radiative transfer equation:  \n-  \nχ  \nI  \nthe intensity and  \ns  \nthe geometrical path length  \nray. The solution can be computed using standard long  \ncharacteristics techniques in decomposed domains (  \net al. 2006; Hayek et al. 2010). Once the intensity is kn  \nFig. 15. Average cooling as function of height in the Bifrost test atmosphere. Total cooling (black), H  (turquoise), Mg  (blue) and Ca  (red). According to the detailed calculations (solid) and according to the recipes (dashed). d I d s = η with I 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 nonLTE 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).  \nheating rate can be directly computed as  \nχ  \nQ chr = with Ω 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.  \n∫  \nΩ  \n-  \nthe solid angle.  \nThe method described above is accurate, but  \nfaster method can be obtained by treating each colum  \nMHD simulation as a plane-parallel atmosphere and a  \nη  \nor  \nχ  \nis zero. For each  \nquantities then become only a function of height  \nz  \nand  \nµ  \n=  \n4  \nπχ  \nJ  \n)d  \n1D  \n,  \nthe cosine of the angle with respect to the ver  \nIf we furthermore assume the emitting region is lo  \ntop of the absorbing region, we can ignore rays that p  \nward from the interior of the star in the integral of  \nThe intensity that impinges on the top of the absorbin  \nfor inward pointing rays is simply the integral of the e  \nof the mean chromospheric energy balance. 7. 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 with µ A39, page 8 of 10 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  \nQ abs = 4 πκ HeI ρJ cor , (52)  \nI  \nd  \nΩ  \nη  \n=  \n4  \nπχ  \nJ  \n3D  \n-  \nη,  \nis a constant dependent on the abundan  \n=  \nαρ  \nN  \nH  \nρ  \nA  \nN  \nHeI  \nN  \nHe  \n,  \nα  \nHeI  \nHe  \nwith κ HeI the representative opacity. The angle-averaged radiation field J cor is computed from the transfer equation:  \nd I d s = η -κ HeI ρI, (53)  \nwith 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.  \nBecause 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.  \nThe 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).  \nFigure 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.", '6 The equation of state and non-equilibrium ionisation': "Non-equilibrium radiative transfer and non-equilibrium ionisation play a role in the solar atmosphere wherever magneto-hydrodynamic timescales are short compared to timescales on which an atomic system reacts to changes. The combined continuity and rate equation for the level population n i in level i of an atomic species with N energy levels is given by:  \nglyph[negationslash]  \nglyph[negationslash]  \n∂n i ∂t + ∇· ( n i v ) = N ∑ j =1 ,j = i n j P ji -n i N ∑ j,j = i P ij , (54)  \nwhere 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  \nFig. 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); Q corona 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.  \n<!-- image -->  \nthen the equations for all levels together form a set of coupled first-order differential equations  \n∂ n ∂t = P n , (55)  \nwith n the vector of level populations and P the rate coefficient matrix. This system has the solution:  \nn ( t ) = N ∑ i =1 c i a i e λ i t , (56)  \nwhere 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  \nilibrium is  \nþ  \nR  \n21NRB  \nÞ  \n:  \n34 CARLSSON & STEIN  \nð  \n14  \nÞ  \nis small compared to  \ns individually, then the  \nradiative rates will make  \non timescale given by  \ntransitions can either  \n(if NRB  \n>  \n0) the relax-  \non timescale (at each  \nnumerical simulation.  \ny: The value of the hydro-  \na given time step and  \non timescale calcula-  \nnsistent with this state of  \nsolving the equations of  \nthe initial population  \nhisatmospherewasthen  \nrature by 1% throughout.  \nthis perturbed state were  \nistical equilibrium giv-  \n1  \nÞ  \n.  \nThe time evolution  \n/C28  \nrelax  \ncal simulation, from the  \nVol. 572  \n<!-- image -->  \norders of magnitude too short. Eigenvalues calculated from a rate matrix  \nof protons at a given the full rate equations, utionwascastinatwoon timescale was calfit to ¼/C0 t : ð 15 Þ Fig. 6. -Relaxation timescale as function of column mass. The numerically determined relaxation timescale is given as the thick solid line. Timescales determined from eigenvalues of the rate matrix are also shown for several cases: full rate matrix ( dotted line ), collisions only ( dot-dashed line ), Lyman transitions in detailed balance ( thin solid line ), and Ly /C11 treated with aconstantnetradiativebracket( dashed line ). 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 Fig. 16 Relaxation timescale for hydrogen ionisation as function of column mass in a 1D radiation hydrodynamics simulation computed with the RADYN code. The numerically determined relaxation timescale is given as the thick solid line. Timescales determined from eigenvalues of the rate matrix are also shown for several cases: full rate matrix (dotted line), collisions only (dot-dashed line), Lyman transitions in detailed balance (thin solid line), and Ly α treated with a constant net radiative bracket (dashed line). The electron density is given as the thick dashed curve. Adapted with permission from Carlsson and Stein (2002), copyright by AAS.  \nwith all Lyman radiative transitions in detailed balance except for Ly  \n/C11  \n,  \nionization/recomwhich is included using its net radiative bracket, give a timescale that closely matches the numerical result. it will evolve towards equilibrium on a characteristic timescale given by  \nis shown in Figure 6  \ncreases outward from the  \nin the midchromosphere  \nτ = 1 min( | λ i | ) , λ i = 0 . (57)  \nglyph[negationslash]  \nhebaseofthetransition of the rate matrix, Pij transition rate per atom the processes controlcan be represented 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 di ff erent assumptions 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).  \n/C21  \ni t  \n:  \nlevel 6 is the continuum.  \nponding to the  \ni  \nth eigen-  \non the initial conditions.  \neigenvector is the equili-  \nare all negative, since the  \nequilibrium value.  \n;  \nð  \n16 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 An illustrative solution for a two-level atom was derived in Carlsson and Stein (2002). The equilibration timescale is given by  \nÞ  \nð 17 Þ 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 )reproducethenumerical 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  \nmagnitude smaller than found in the numerical solution ( thick solid line ). The timescale obtained from the rate matrix with only collisional rates ( dot-dashed line )hasthe τ = 1 /( C 21 + C 12 + R 21 [ 1 -n 1 n 2 R 12 R 21 ]) . (58)  \ntimescale up to the peak in the midchromosphere. This  \nindicates that in this region the large Lyman radiative  \nrates are in fact changing with time so as to maintain  \nnearly detailed balance as the populations change. The  \ngeneral behavior of the relaxation time is thus controlled  \nby the collisional processes. The timescale increases from  \nthe photosphere to the midchromosphere, where the elec-  \ntron density has a minimum approximately as the inverse  \nof the collisional recombination rate (  \nn  \n/C0  \n2  \ne  \n).  \nFrom the  \nmidchromosphere to the transition region the timescale  \n/  \ntions,  \nexponentially 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.  \nCarlsson 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.  \nA 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  \nThese 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.  \nThe 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  \nn e = ∑ i,j,k ( j -1) n ijk , (59)  \nand energy conservation  \ne = 3 2 k B T   ∑ i,j,k n ijk + n e   + ∑ i,j,k n ijk E ijk . (60)  \nHere 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.  \nSollum (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.  \n630  \nJ. Leenaarts et al.: Non-equilibrium hydrogen ionization in 2D simulations of the solar atmosphere  \nFig. 3. Time slices of the gas temperature ( first row ), the ionization degree of hydrogen ( second row ), number density of hydrogen in the n = 2level ( third row ), and the n = 2columndensity( fourth row )inamagneticelement( left-hand column ) and in the internetwork ( right-hand column ). The upper-left magnetic element panel shows dynamic fibrils pushing the corona upward with 3 min periodicity. The upper-right internetwork panel shows rather unstructured shocks and a slowly varying height of the transition region. The snapshot used in Figs. 1 and 2 is indicated by a black dotted line. 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 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 raFig. 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.  \n<!-- image -->  \ngas particles.  \ndiative cooling in the chromosphere and transition region, where  \ndeviations from equilibrium are largest, is thus inconsistent with  \nFig. 3. Time slices of the gas temperature ( first row ), the ionization degree of hydrogen ( second row ), number density of hydrogen in the n = 2level ( third row ), and the n = 2columndensity( fourth row )inamagneticelement( left-hand column ) and in the internetwork ( right-hand column ). The upper-left magnetic element panel shows dynamic fibrils pushing the corona upward with 3 min periodicity. The upper-right internetwork panel shows rather unstructured shocks and a slowly varying height of the transition region. The snapshot used in Figs. 1 and 2 is indicated by a black dotted line. 4. Discussion and conclusions 4.1. Limitationsofthesimulation 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 simthe 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. 4.2. Discussion 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.  \nOur 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 lation. 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. 4.2. Discussion From the analysis of our simulation we obtain the following picture. The internetwork chromosphere is irregularly pervaded 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.  \nof 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. 4. Discussion and conclusions 4.1. Limitationsofthesimulation 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 simuulation, which sets the radiative cooling and heating, employs LTE ionization. For given internal energy and mass density, the From the analysis of our simulation we obtain the following picture. The internetwork chromosphere is irregularly pervaded 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.  \nThe increased realism of an EOS that includes non-equilibrium ionisation comes at a rather steep price in computational efficiency: A simulation with a  \nnon-equilibrium EOS is around three (hydrogen only) to five times (hydrogen and helium) slower than a similar simulation with an LTE EOS.  \nNon-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.  \nIn 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.  \nIf 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).  \nLow-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  \nAbbett 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:  \nJ ν ( τ ν ) = 1 2 ∫ ∞ 0 S ν ( τ ' ν ) E 1 ( | τ ν -τ ' ν | ) d τ ' ν . (61)  \nThe Astrophysical Journal,  \n817:125  \n(  \n13pp  \n)  \n, 2016 February 1  \nGolding, Leenaarts, & Carlsson  \nFigure 2  \n<!-- image -->  \nFigure 8.  \nJoint probability density function of height and logarithmic temperature. The  \n/uniFB01  \ngure includes data from the time interval 1000  \n-  \n3000 s. The three horizontal  \nplateaus ( at 6, 10, and 22 kK ) in the LTE simulation indicate preferred temperatures when using the LTE equation of state. These temperatures are associated with the LTE ionization of H I, He I, and He II. The plateaus vanish when we introduce non-equilibrium hydrogen and helium ionization. Fig. 18 Joint probability density functions of height and temperature in four radiationMHD simulations of the solar atmosphere. The simulations differ in their treatment of the equation of state. LTE: LTE equation of state; HION: hydrogen in non-equilibrium, helium in LTE, hydrogen Lyman transitions in detailed balance. LYA-HION: as HION but with radiative transfer in Lyman transitions; HELIUM: both hydrogen and helium in nonequilibrium. The three horizontal plateaus (at 6, 10, and 22 kK) in the LTE simulation indicate preferred temperatures when using the LTE equation of state. These temperatures are associated with the LTE ionisation of H I, He I, and He II. The plateaus vanish when we introduce non-equilibrium hydrogen and helium ionisation. Adapted with permission from Golding et al. (2016). 14  \nFigure 9.  \nDifferential  \n1000  \n-  \n3000 s. The HELIUM run DEM does not have a bulge associated with  \nthe preferred temperature at  \nT  \n/uniF0A0  \nvalue than the other three runs in the temperature range 11  \nFigure 1 Left: Temperature at the RADMHD model photosphere. Right: Magnetic-field lines threading the low atmosphere over a small subdomain (the box in the left frame indicates the approximate size of the corresponding subdomain). The gray slice indicates the average height of the visible surface. The domain spans 24 × 24 × 12 Mm 3 at a resolution of 512 × 512 × 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.  \n<!-- image -->  \nTemperature in the  \nRADMHD  \nlow chromosphere showing a reverse granulation pattern. Lighter  \n(darker) colors indicate hotter (cooler) temperatures. In the models, this occurs because the radiative cooling  \ndiminishes with height, and the  \np  \nv  \nwork of converging and diverging flows above the intergranular lanes  \nbegins to dominate. The horizontal slice spans 21  \n×  \n12 Mm  \n2  \nat a resolution of 448  \n×  \n256.  \n∇  \n·  \nOccurrence of helium ion fractions as a function of temperature for  \nHe II in the top panel, and He III in the bottom panel. Each column has been  \nnormalized to increase readability. Median values and CHIANTI values are  \nW.P. Abbett, G.H. Fisher  \nThe first exponential integral E 1 ( x ) peaks sharply at x = 0, so that this expression can be approximated as  \nJ ν ( τ ν ) ≈ 1 2 S ν ( τ ν ) ∫ ∞ 0 E 1 ( | τ ν -τ ' ν | ) d τ ' ν (62)  \n≈ S ν ( τ ν ) ( 1 -E 2 ( τ ν ) 2 ) (63)  \nThis result is inserted into the equation for the radiative cooling (Eq. (11)) to yield  \nIf 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:  \nQ ≈ -2 πρ ∫ ∞ 0 κ ν S ν E 2 ( τ ν ) d ν. (64)  \nQ ≈ -2 κ B ρσT 4 E 2 ( τ B ) , (65)  \nwhere κ 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 .  \nFigure 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.", '7.2 Escape probability method': 'An interesting new method that speeds up non-equilibrium radiative transfer was proposed by Judge (2017). The major time consuming task in radiative transfer is the evaluation of the formal solution for all required frequencies and angles in order to construct the angle-averaged and frequency-averaged radiation field  \n¯ J = 1 4 π ∫ ∞ 0 ∫ Ω I ν d Ω d ν, (66)  \nwhich appears in the radiative rate coefficients of Eq. (54).  \nJudge (2017) suggests to compute ¯ J from the equation  \nd d τ ( S -¯ J ) = q 1 / 2 d d τ ( q 1 / 2 S ) , (67)  \n(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:  \nthe 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.  \nThe big advantage of using Eq. (67) is that it replaces the repeated evaluations of the transfer equation in order to compute ¯ J in a transition with the solution of a single integral along vertical columns in the atmosphere. Solving the statistical equilibrium non-LTE radiative transfer problem problem in a MURaM test atmosphere is ∼ 100 times faster than the full method.  \nJudge (2017) shows that this method can also be used to solve non-equilibrium problems (i.e, it can be used to simultaneously solve Eqs. (54), (59), and (60)). Radiation-MHD simulations using this method to compute chromospheric radiative energy exchange or the equation state have so far not been reported on. It would be very interesting to see whether the method is fast enough to be used in practice and how it compares to the methods described in Sect. 6.', '8 Conclusions and outlook': "In this review I presented the most commonly used methods to approximate the transfer of energy between solar plasma and the radiation field in radiationMHD 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.  \nRadiative 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.  \nRadiative 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-  \nFig. 20 Brightness temperature in the nominal line core of H α in a simulated active region. The radiation-MHD simulation was performed with the MURaM code by Danilovic et al. (2019). The simulation used a single-bin LTE treatment of chromospheric radiative transfer. Nevertheless, the simulation resembles observations rather well. The 3D non-LTE radiative transfer calculation used to obtain the H α intensity was done with the Multi3d code.  \n<!-- image -->  \ntions 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,  \nInterestingly, 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).  \nI 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.  \nHere 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.  \nIn 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.  \nWhen 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.  \nIt 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.  \nThe 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.  \nSo what can be improved? I propose the following non-exhaustive list:  \n- -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.\n- -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).\n- -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 α .\n- -The escape probability method for approximating the chromospheric radiation field discussed in Sect. 7.2 should be tested and developed further.\n- -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.  \nLet us hope that improvements along these lines, as well as others, will be presented during the coming years.  \nAcknowledgements This work was supported by grants from the Knut and Alice Wallenberg foundation (2016.0019) and the Swedish Research Council (2017-04099) and benefitted from discussions within the activities of team 399 'Studying magnetic-field-regulated heating in the solar chromosphere' at the International Space Science Institute (ISSI) in Switzerland.", 'References': "Abbett WP, Fisher GH (2012) Radiative Cooling in MHD Models of the Quiet Sun Convection Zone and Corona. Sol. Phys.277(1):3-20, DOI 10.1007/s11207-011-9817-3, 1102.1035  \nAschwanden MJ, Fletcher L, Schrijver CJ, Alexander D (1999) Coronal Loop Oscillations Observed with the Transition Region and Coronal Explorer. ApJ520:880-894, DOI 10.1086/307502  \n- Asplund M, Grevesse N, Sauval AJ, Allende Prieto C, Kiselman D (2004) Line formation in solar granulation. IV. [O I], O I and OH lines and the photospheric O abundance. A&A417:751-768, DOI 10.1051/0004-6361:20034328, astro-ph/0312290\n- Asplund M, Grevesse N, Sauval AJ, Allende Prieto C, Blomme R (2005) Line formation in solar granulation. VI. [C I], C I, CH and C 2 lines and the photospheric C abundance. A&A431:693-705, DOI 10.1051/0004-6361:20041951, astro-ph/0410681  \nAuer L (2003) Insight into Multi-Dimensional Transfer. In: Hubeny I, Mihalas D, Werner K (eds) Stellar Atmosphere Modeling, Astronomical Society of the Pacific, San Francisco, ASP Conference Series, vol 288, p 405  \nBjørgen JP, Sukhorukov AV, Leenaarts J, Carlsson M, de la Cruz Rodr'ıguez J, Scharmer GB, Hansteen VH (2018) Three-dimensional modeling of the Ca II H and K lines in the solar atmosphere. A&A611:A62, DOI 10.1051/0004-6361/201731926, 1712.01045  \nBjørgen JP, Leenaarts J, Rempel M, Cheung MCM, Danilovic S, de la Cruz Rodr'ıguez J, Sukhorukov AV (2019) Three-dimensional modeling of chromospheric spectral lines in a simulated active region. A&A631:A33, DOI 10.1051/0004-6361/201834919, 1906.01098  \nBradshaw SJ, Del Zanna G, Mason HE (2004) On the consequences of a non-equilibrium ionisation balance for compact flare emission and dynamics. A&A425:287-299, DOI 10.1051/0004-6361:20040521  \nBruls JHMJ, Vollmoller P, Schussler M (1999) Computing radiative heating on unstructured spatial grids. A&A348:233-248  \n- Carlson BG (1963) The numerical theory of neutron transport. In: Alder B, Fernbach S (eds) Methods of Computational Physics, vol. 1, Academic Press, New York, pp 1-42\n- Carlsson M, Leenaarts J (2012) Approximations for radiative cooling and heating in the solar chromosphere. A&A539:A39, DOI 10.1051/0004-6361/201118366, 1202.2996\n- Carlsson M, Stein RF (1992) Non-LTE Radiating Acoustic Shocks and CA II K2V Bright Points. ApJ397:L59, DOI 10.1086/186544\n- Carlsson M, Stein RF (2002) Dynamic Hydrogen Ionization. ApJ572(1):626-635, DOI 10. 1086/340293, astro-ph/0202313\n- Carlsson M, Hansteen VH, Gudiksen BV, Leenaarts J, De Pontieu B (2016) A publicly available simulation of an enhanced network region of the Sun. A&A585:A4, DOI 10. 1051/0004-6361/201527226, 1510.07581\n- Cheung MCM, Cameron RH (2012) Magnetohydrodynamics of the Weakly Ionized Solar Photosphere. ApJ750:6, DOI 10.1088/0004-637X/750/1/6, 1202.1937  \nCheung MCM, Rempel M, Chintzoglou G, Chen F, Testa P, Mart'ınez-Sykora J, Sainz Dalda A, DeRosa ML, Malanushenko A, Hansteen V, De Pontieu B, Carlsson M, Gudiksen B, McIntosh SW (2019) A comprehensive three-dimensional radiative magnetohydrodynamic simulation of a solar flare. Nature Astron 3:160-166, DOI 10.1038/ s41550-018-0629-3  \n- Cranmer SR, Winebarger AR (2019) The Properties of the Solar Corona and Its Connection to the Solar Wind. ARA&A57:157-187, DOI 10.1146/annurev-astro-091918-104416, 1811.00461\n- Curdt W, Landi E, Feldman U (2004) The SUMER spectral atlas of solar coronal features. A&A427:1045-1054, DOI 10.1051/0004-6361:20041278\n- de la Cruz Rodr'ıguez J, van Noort M (2017) Radiative Diagnostics in the Solar Photosphere and Chromosphere. Space Sci. Rev.210(1-4):109-143, DOI 10.1007/s11214-016-0294-8, 1609.08324\n- de la Cruz Rodr'ıguez J, Leenaarts J, Danilovic S, Uitenbroek H (2019) STiC: A multiatom non-LTE PRD inversion code for full-Stokes solar observations. A&A623:A74, DOI 10.1051/0004-6361/201834464, 1810.08441  \nDere KP, Landi E, Mason HE, Monsignori Fossi BC, Young PR (1997) CHIANTI - an atomic database for emission lines. A&AS125:149-173, DOI 10.1051/aas:1997368  \nDere KP, Landi E, Young PR, Del Zanna G, Landini M, Mason HE (2009) CHIANTI an atomic database for emission lines. IX. Ionization rates, recombination rates, ionization equilibria for the elements hydrogen through zinc and updated atomic data. A&A498(3):915-929, DOI 10.1051/0004-6361/200911712  \nDere KP, Del Zanna G, Young PR, Landi E, Sutherland RS (2019) CHIANTI-An Atomic Database for Emission Lines. XV. Version 9, Improvements for the X-Ray Satellite Lines. ApJS241(2):22, DOI 10.3847/1538-4365/ab05cf, 1902.05019  \nDzifˇc'akov'a E, Dud'ık J, Mackovjak ˇ S (2016) Non-equilibrium ionization by a periodic electron beam. I. Synthetic coronal spectra and implications for interpretation of observations. A&A589:A68, DOI 10.1051/0004-6361/201527557, 1603.04256  \nFeautrier P (1964) Sur la r'esolution num'erique de l'equation de transfert. C R Acad Sci 258:3189-3191  \nFeldman U, Mandelbaum P, Seely JF, Doschek GA, Gursky H (1992) The Potential for Plasma Diagnostics from Stellar Extreme-Ultraviolet Observations. ApJS81:387, DOI 10.1086/191698  \nFreytag B, Steffen M, Ludwig HG, Wedemeyer-Bohm S, Schaffenberger W, Steiner O (2012) Simulations of stellar convection with CO5BOLD. J Comput Phys 231:919-959, DOI 10.1016/j.jcp.2011.09.026, 1110.6844 √  \nFrisch U, Frisch H (1975) Non-LTE transfer. glyph[epsilon1] revisited. MNRAS173:167-182, DOI 10. 1093/mnras/173.1.167  \nGolding TP, Carlsson M, Leenaarts J (2014) Detailed and Simplified Nonequilibrium Helium Ionization in the Solar Atmosphere. ApJ784(1):30, DOI 10.1088/0004-637X/784/1/30, 1401.7562  \n- Golding TP, Leenaarts J, Carlsson M (2016) Non-equilibrium Helium Ionization in an MHD Simulation of the Solar Atmosphere. ApJ817(2):125, DOI 10.3847/0004-637X/817/2/ 125, 1512.04738\n- Golding TP, Leenaarts J, Carlsson M (2017) Formation of the helium extreme-UV resonance lines. A&A597:A102, DOI 10.1051/0004-6361/201629462, 1610.00352\n- Graham DR, Cauzzi G (2015) Temporal Evolution of Multiple Evaporating Ribbon Sources in a Solar Flare. ApJ807:L22, DOI 10.1088/2041-8205/807/2/L22, 1506.03465\n- Grevesse N, Sauval AJ (1998) Standard Solar Composition. Space Sci. Rev.85:161-174, DOI 10.1023/A:1005161325181\n- Gudiksen BV, Carlsson M, Hansteen VH, Hayek W, Leenaarts J, Mart'ınez-Sykora J (2011) The stellar atmosphere simulation code Bifrost. Code description and validation. A&A531:A154, DOI 10.1051/0004-6361/201116520, 1105.6306  \nGustafsson B, Bell RA, Eriksson K, Nordlund ˚ A (1975) A grid of model atmospheres for metal-deficient giant stars. I. A&A42:407-432  \nHansteen V (1993) A new interpretation of the redshift observed in optically thin transition region lines. ApJ402:741-755, DOI 10.1086/172174  \nHayek W, Asplund M, Carlsson M, Trampedach R, Collet R, Gudiksen BV, Hansteen VH, Leenaarts J (2010) Radiative transfer with scattering for domain-decomposed 3D MHD simulations of cool stellar atmospheres. Numerical methods and application to the quiet, non-magnetic, surface of a solar-type star. A&A517:A49, DOI 10.1051/0004-6361/201014210, 1007.2760  \nHeinemann T, Dobler W, Nordlund ˚ A, Brandenburg A (2006) Radiative transfer in decomposed domains. A&A448:731-737, DOI 10.1051/0004-6361:20053120, astro-ph/ 0503510  \nHubeny I, Mihalas D (2014) Theory of Stellar Atmospheres. Princeton University Press, Princeton, NJ  \n- Hummer DG, Rybicki GB (1982) Second-order escape probability approximations in radiative transfer. ApJ263:925-934, DOI 10.1086/160560  \nIbgui L, Hubeny I, Lanz T, Stehl'e C (2013) IRIS: a generic three-dimensional radiative transfer code. A&A549:A126, DOI 10.1051/0004-6361/201220468, 1211.4870  \nIijima H (2016) Numerical studies of solar chromospheric jets. PhD thesis, Department of Earth and Planetary Science, School of Science, The University of Tokyo, Japan, DOI 10.5281/zenodo.55411  \nIijima H, Yokoyama T (2015) Effect of Coronal Temperature on the Scale of Solar Chromospheric Jets. ApJ812:L30, DOI 10.1088/2041-8205/812/2/L30, 1509.06677  \nJudge PG (2005) Understanding the time dependence of atomic level populations in evolving plasmas. J. Quant. Spec. Radiat. Transf.92(4):479-510, DOI 10.1016/j.jqsrt.2004.08.009 Judge PG (2017) Efficient Radiative Transfer for Dynamically Evolving Stratified Atmospheres. ApJ851:5, DOI 10.3847/1538-4357/aa96a9  \nKaˇsparov'a J, Varady M, Heinzel P, Karlick'y M, Moravec Z (2009) Response of optical hydrogen lines to beam heating. I. Electron beams. A&A499:923-934, DOI 10.1051/ 0004-6361/200811559, 0904.2084  \n- Kerr GS, Carlsson M, Allred JC, Young PR, Daw AN (2019) Si IV Resonance Line Emission during Solar Flares: Non-LTE, Nonequilibrium, Radiation Transfer Simulations. ApJ871(1):23, DOI 10.3847/1538-4357/aaf46e, 1811.11075\n- Khomenko E, Collados M, D'ıaz A, Vitas N (2014) Fluid description of multi-component solar partially ionized plasma. Phys Plasmas 21(9):092901, DOI 10.1063/1.4894106, 1408. 1871\n- Khomenko E, Vitas N, Collados M, de Vicente A (2017) Numerical simulations of quiet Sun magnetic fields seeded by the Biermann battery. A&A604:A66, DOI 10.1051/0004-6361/ 201630291, 1706.06037\n- Khomenko E, Vitas N, Collados M, de Vicente A (2018) Three-dimensional simulations of solar magneto-convection including effects of partial ionization. A&A618:A87, DOI 10.1051/0004-6361/201833048, 1807.01061\n- Kunasz P, Auer LH (1988) Short characteristic integration of radiative transfer problems: formal solution in two-dimensional slabs. J. Quant. Spec. Radiat. Transf.39:67-79, DOI 10.1016/0022-4073(88)90021-0\n- Kurucz RL (1979) Model atmospheres for G, F, A, B, and O stars. ApJS40:1-340, DOI 10.1086/190589\n- Landi E, Landini M (1999) Radiative losses of optically thin coronal plasmas. A&A347:401408\n- Le Chat G, Issautier K, Meyer-Vernet N (2012) The Solar Wind Energy Flux. Sol. Phys.279(1):197-205, DOI 10.1007/s11207-012-9967-y, 1203.1316\n- Leenaarts J, Carlsson M (2009) MULTI3D: A Domain-Decomposed 3D Radiative Transfer Code. In: Lites B, Cheung M, Magara T, Mariska J, Reeves K (eds) The Second Hinode Science Meeting: Beyond Discovery-Toward Understanding, Astronomical Society of the Pacific, San Francisco, ASP Conference Series, vol 415, p 87\n- Leenaarts J, Carlsson M, Hansteen V, Rutten RJ (2007) Non-equilibrium hydrogen ionization in 2D simulations of the solar atmosphere. A&A473:625-632, DOI 10.1051/ 0004-6361:20078161, 0709.3751\n- Leenaarts J, Pereira TMD, Carlsson M, Uitenbroek H, De Pontieu B (2013) The Formation of IRIS Diagnostics. II. The Formation of the Mg II h&k Lines in the Solar Atmosphere. ApJ772(2):90, DOI 10.1088/0004-637X/772/2/90, 1306.0671\n- Leitner P, Lemmerer B, Hanslmeier A, Zaqarashvili T, Veronig A, Grimm-Strele H, Muthsam HJ (2017) Structure of the solar photosphere studied from the radiation hydrodynamics code ANTARES. Ap&SS362(9):181, DOI 10.1007/s10509-017-3151-7, 1708.01156  \nLudwig HG (1992) Nichtgrauer Strahlungstransport in numerischen Simulationen stellarer Konvektion. PhD thesis, CHristian-Albrechts-Universitat  \n- Mart'ınez-Sykora J, De Pontieu B, Hansteen V (2012) Two-dimensional Radiative Magnetohydrodynamic Simulations of the Importance of Partial Ionization in the Chromosphere. ApJ753:161, DOI 10.1088/0004-637X/753/2/161, 1204.5991  \nMart'ınez-Sykora J, De Pontieu B, Hansteen VH, Rouppe van der Voort L, Carlsson M, Pereira TMD (2017) On the generation of solar spicules and Alfv'enic waves. Science 356(6344):1269-1272, DOI 10.1126/science.aah5412, 1710.07559  \n- Mihalas D (1970) Stellar atmospheres. Series of Books in Astronomy and Astrophysics, WH Freeman, San Francisco  \nMihalas D, Mihalas BW (1984) Foundations of radiation hydrodynamics. Oxford University  \nPress, New York  \n- Nordlund ˚ A (1982) Numerical simulations of the solar granulation. I. Basic equations and methods. A&A107:1-10\n- Olluri K, Gudiksen BV, Hansteen VH (2013) Non-equilibrium Ionization in the Bifrost Stellar Atmosphere Code. AJ145:72, DOI 10.1088/0004-6256/145/3/72\n- Olson GL, Auer LH, Buchler JR (1986) A rapidly convergent iterative solution of the nonLTE radiation transfer problem. J. Quant. Spec. Radiat. Transf.35:431-442, DOI 10. 1016/0022-4073(86)90030-0\n- Pascoe DJ, Anfinogentov SA, Goddard CR, Nakariakov VM (2018) Spatiotemporal Analysis of Coronal Loops Using Seismology of Damped Kink Oscillations and Forward Modeling of EUV Intensity Profiles. ApJ860:31, DOI 10.3847/1538-4357/aac2bc\n- Pereira TMD, Asplund M, Collet R, Thaler I, Trampedach R, Leenaarts J (2013) How realistic are solar model atmospheres? A&A554:A118, DOI 10.1051/0004-6361/201321227,", '1304.4932': "Popovas A, Nordlund ˚ A, Ramsey JP (2019) Pebble dynamics and accretion on to rocky planets - II. Radiative models. MNRAS482(1):L107-L111, DOI 10.1093/mnrasl/sly197, 1810.07048  \nRempel M (2017) Extension of the MURaM Radiative MHD Code for Coronal Simulations. ApJ834:10, DOI 10.3847/1538-4357/834/1/10, 1609.09818  \nRutten RJ, de Wijn AG, Sutterlin P (2004) DOT tomography of the solar atmosphere. II. Reversed granulation in Ca II H. A&A416:333-340, DOI 10.1051/0004-6361:20035636  \nSkartlien  \nR (2000) A Multigroup Method for Radiation with Scattering in Three-  \nDimensional Hydrodynamic Simulations. ApJ536:465-480, DOI 10.1086/308934  \nSollum E (1999) Dynamic hydrogen ionization. Master's thesis, University of Oslo  \nStein RF (2012) Solar Surface Magneto-Convection. Living Rev Sol Phys 9:4, DOI 10.12942/ lrsp-2012-4  \nStein RF, Nordlund ˚ A (1998) Simulations of Solar Granulation. I. General Properties. ApJ499:914, DOI 10.1086/305678  \nSukhorukov AV, Leenaarts J (2017) Partial redistribution in 3D non-LTE radiative transfer in solar-atmosphere models. A&A597:A46, DOI 10.1051/0004-6361/201629086, 1606. 05180  \nSummers HP (1972) The density dependent ionization balance of carbon, oxygen and neon in the solar atmosphere. MNRAS158:255, DOI 10.1093/mnras/158.3.255  \nSummers HP (1974) The ionization equilibrium of hydrogen-like to argon-like ions of elements. MNRAS169:663-680, DOI 10.1093/mnras/169.3.663  \nTrujillo Bueno J, Fabiani Bendicho P (1995) A Novel Iterative Scheme for the Very Fast and Accurate Solution of Non-LTE Radiative Transfer Problems. ApJ455:646, DOI 10.1086/176612  \nTurck-Chi'eze S, Couvidat S (2011) Solar neutrinos, helioseismology and the solar internal dynamics. Rep Progr Phys 74(8):086901, DOI 10.1088/0034-4885/74/8/086901, 1009. 0852  \nvan Regemorter H (1962) Rate of Collisional Excitation in Stellar Atmospheres. ApJ136:906, DOI 10.1086/147445  \nVernazza JE, Avrett EH, Loeser R (1981) Structure of the solar chromosphere. III. Models of the EUV brightness components of the quiet-sun. ApJS45:635-725, DOI 10.1086/190731 Vogler A (2004) Three-dimensional simulations of magneto-convection in the solar photosphere. PhD thesis, Georg-August-Universitat, Gottingen  \nVogler A, Bruls JHMJ, Schussler M (2004) Approximations for non-grey radiative transfer in numerical simulations of the solar photosphere. A&A421:741-754, DOI 10.1051/ 0004-6361:20047043  \nVogler A, Shelyag S, Schussler M, Cattaneo F, Emonet T, Linde T (2005) Simulations of magneto-convection in the solar photosphere. Equations, methods, and results of the MURaM code. A&A429:335-351, DOI 10.1051/0004-6361:20041507  \nˇ Stˇep'an J, Trujillo Bueno J (2013) PORTA: A three-dimensional multilevel radiative transfer code for modeling the intensity and polarization of spectral lines with massively parallel computers. A&A557:A143, DOI 10.1051/0004-6361/201321742, 1307.4217  \nWoods TN, Eparvier FG, Hock R, Jones AR, Woodraska D, Judge D, Didkovsky L, Lean J, Mariska J, Warren H (2012) Extreme Ultraviolet Variability Experiment (EVE) on the Solar Dynamics Observatory (SDO): Overview of Science Objectives, Instrument Design, Data Products, and Model Developments. Sol. Phys.275(1-2):115-143, DOI 10.1007/s11207-009-9487-6  \nWray AA, Bensassi K, Kitiashvili IN, Mansour NN, Kosovichev AG (2015) Simulations of Stellar Magnetoconvection using the Radiative MHD Code 'StellarBox'. arXiv e-prints 1507.07999"}