© ESO 2018

1. Introduction

Solar irradiance varies on short timescales from minutes to hours as well as long time scales of days, years, decades and beyond. The variations at the different time scales originate from different processes in the solar atmosphere. The short time scales – other than reconnection processes such as flares – are mostly determined by convection and the longer ones are largely driven by the changing solar surface magnetic field (see e.g., Domingo et al. 2009; Solanki et al. 2013; Yeo et al. 2017).

The results over the past decades clearly indicate that solar variability has an influence on the Earth’s climate; for an overview see e.g., Matthes et al. (2017) and Shindell et al. (2020). However, the exact quantification of the solar influence on climate – besides other natural forcings and the anthropogenic contribution – is still debated. Egorova et al. (2018a) conclude that the Sun must have varied substantially, in order to attribute the temperature increase in the early 19^th century.

To quantify the solar contribution to climate change more precisely, robust solar irradiance datasets are needed as input to the climate models. For those times when space observations provide a decent temporal and spectral coverage, a number of observational composite datasets are available (see e.g., Haberreiter et al. 2017; Coddington et al. 2019; Marchenko et al. 2019). However, when no satellite observations are available irradiance reconstruction models using proxy datasets to describe the state of solar activity back in time are key for our understanding of its impact on climate.

The extent to which solar variability has changed over long time scales is still an open question. Shapiro et al. (2011) proposed a reconstruction of the Total Solar Irradiance (TSI) and Spectral Solar Irradiance (SSI), which accounts for a variable quiet Sun intensity. This approach leads to a significant change of the radiative forcing between the Maunder Minimum and the space era of about 6 ± 3 W m⁻². A later update of this approach suggests a smaller change in long-term irradiance variability (Egorova et al. 2018b). However, other reconstructions methods give a rather flat long-term trend; see e.g., Solanki et al. (2013). Recently, Rempel (2020) determined a linear dependence between the outgoing radiative energy flux and the mean magnetic field strength for the quiet Sun. While the strong sensitivity of TSI to the quiet Sun field strength implies that potential variations of the magnetic field over longer timescales could make a significant contribution to solar irradiance variations, such magnetic variations are not expected if the quiet Sun magnetic field originates primarily from a small-scale dynamo. Nevertheless, to determine the irradiance variations correctly, it is key to quantify the emergent spectrum for various solar surface magnetic elements with a radiative transfer code. State-of-the-art radiative transfer codes are however inherently different with respect to the atomic input data and different numerical schemes to solve the radiative transfer equation. This introduces some uncertainty for irradiance reconstructions. Criscuoli et al. (2020) recently validated the performance of a set of commonly used codes, i.e., the radiative transfer codes COSI (Haberreiter et al. 2008b; Shapiro et al. 2010; Criscuoli et al. 2020), the Rybicki-Hummer (RH) code (Uitenbroek 2003), and the radiation scheme of the Max Planck University of Chicago Radiative MHD (MURaM) code (Rempel 2014). The authors find a good agreement of the spectral synthesis using 1D solar atmosphere structures. In the present paper we go beyond the study by Criscuoli et al. (2020) and use the vertical temperature and density profiles from three 3D MHD simulations, a hydrodynamic (HD) case with 0 G magnetic field, and a 100 G and 200 G MHD snapshot as input for the spectral synthesis to the above-mentioned radiative transfer codes in order to assess their performance. This is deemed particularly necessary, as radiative transfer codes developed to work in 3D geometry may differ not only for the reasons mentioned above, but also for the numerical schemes employed to resolve discontinuities and strong gradients typically present in atmospheres generated by MHD codes (see e.g., Janett 2019). Similar work has already been carried out by Afram et al. (2011) and Norris et al. (2017) by applying a set of simulations with different magnetic field strengths from the MURaM Code (Vögler et al. 2005) as input for the spectral synthesis. To our knowledge, this is the first time that spectral syntheses from 3D MHD simulations carried out by different radiative transfer codes are compared in a quantitative manner.

The paper is organized as follows. First, in Sect. 2 we introduce the 3D MHD simulations used to calculate the synthetic spectra. Second, in Sect. 3 we describe the radiative transfer codes and in Sect. 4 we discuss the results. Finally, in Sect. 5 we summarize our findings.

2. MHD simulations

For this study we calculate the emergent intensities for different snapshots from the simulation runs with the MURaM code (Rempel 2014) for the non-magnetic, or hydrodynamic (HD) case, and for a magnetic field of 100 G and 200 G, respectively. The magnetic cases were branched from the HD setup by adding a uniform vertical field of 100 and 200 G. The simulations ran for about an hour simulated solar time until a statistically relaxed state was reached. These simulations consist of cubes of 384 × 384 × 96 pixels³, corresponding to an area on the solar surface of 8.8 × 8.8 arcsec². The vertical domain extent is 1536 km, with the average τ = 1 level located at about 684 km beneath the top boundary. For further details we refer to Criscuoli et al. (2020). Figure 1, top panels, show the horizontal distribution of the line-of-sight magnetic field strength at the τ = 1 layer, and the temperature variation at the top of the simulation box. The bottom panels give an example of the horizontal variation of the temperature. The depth stratification of the temperature and density are used as input to the calculations discussed in Sect. 3.

Fig. 1. Top panels: photospheric magnetic field strength Bz, τ = 1 for the snapshots of the pure hydrodynamic case (left panel), and the cases with a magnetic field strength of 100 G (middle panel), and 200 G (right panel); bottom panels: horizontal temperature variation of the snapshots at the top of the simulation box for the pure hydrodynamical case (left panel), and the cases with a magnetic field strength of 100 G (middle panel), and 200 G (right panel).

3. Spectral synthesis codes

For the spectral synthesis based on the 3D MHD simulations we use three different radiative transfer codes. First, the COde for Solar Irradiance (COSI, Haberreiter et al. 2008a,b; Shapiro et al. 2010; Criscuoli et al. 2020) allows us to calculate the atomic level populations and the emergent intensity taking into account non-local thermodynamic equilibrium (non-LTE), which is a key element for the calculation of the UV spectral range. The atomic data used in COSI are explained in detail in Shapiro et al. (2010) and the numerical scheme of the radiative transfer goes back to Hamann & Schmutz (1987) and Hubeny (1981). The scheme was first applied to solar studies by Haberreiter et al. (2008a,b). So far, COSI calculations have been based on vertical atmosphere structures, such as the 1D atmosphere structures by Fontenla et al. (1999). In this work, we use for the first time 3D MHD simulations as input to COSI.

Second, we make use of the RH code (Uitenbroek 2001). This radiative transfer code can compute emergent intensities at different viewing angles in different geometries. It allows the computation of several atomic and molecular transitions in both LTE and non-LTE under complete or partial redistribution. Because of its versatility, RH is widely employed for spectral and spectro-polarimetric syntheses of atomic and molecular lines in solar and stellar atmospheres, and more recently it has been employed for solar irradiance reconstructions (Criscuoli et al. 2018; Criscuoli 2019; Berrilli et al. 2020). A massively parallel version, RH1.5D, allows us to compute emergent radiation on a column-by-column basis (Pereira & Uitenbroek 2015) and has been used for syntheses from 3D MHD simulations (e.g., Pereira et al. 2013; Antolin et al. 2018; Peck et al. 2019). We make use of RH1.5D to allow for the efficient calculation of intensities from the simulations. Because we focus on LTE calculations for the vertically emergent intensity, there is no difference between a 1.5D calculation and a full 3D calculation.

Third, the MURaM code is a radiative MHD code that was originally developed by Vögler et al. (2005). We use here the version of Rempel (2014) that uses a different formulation of numerical diffusivities and has been tweaked for computational performance. The radiative transfer scheme is identical to Vögler et al. (2005), but has been expanded to allow for additional diagnostic radiative transfer as described in Rempel (2020). In this study we use MURaM in two capacities: (1) to compute the MHD snapshots analysed and (2) to use the radiative transfer solver of MURaM in order to compute diagnostic intensities for comparison. To this end we use the approach detailed in Rempel (2020) in combination with a RH based opacity table from Criscuoli et al. (2020). The MURaM RT scheme (Vögler et al. 2005) uses short characteristics integration as described in Kunasz & Auer (1988). Since gradients in opacity, density and source function are steep, MURaM uses a linear interpolation for enhanced stability. This can lead to artificial broadening of inclined rays as discussed by Peck et al. (2017). Since we focus in this paper only on the intensity for vertical rays (i.e., coordinate axis aligned rays), intensity diffusion is not a concern. Finally, as already indicated above, MURaM uses the same opacities as RH. The main purpose of this paper is to benchmark the radiative transfer codes and to validate their performance. For the present study to be consistent all radiative transfer calculations are done in LTE. The basic concept of the radiative transfer equations and its solution is given in the Appendix. We are specifically interested in validating the codes for the continuum wavelength. We have identified the wavelength of 665.01 nm to be free of spectral lines, and therefore adopt it for this study, hereby referring to it as 665 nm.

4. Results

For each column in the MHD simulation box we have calculated the intensity spectrum at 665 nm in LTE. For the comparison of the results it turned out to be important to use an identical grid for all three radiative transfer codes. This is particularly important as the MURaM radiation scheme is inherently set up to use a grid that is shifted from the MHD grid by half a grid point in all three dimensions (i.e., MHD quantities are cell-centred while intensities are computed on cell corners). Data are interpolated onto this grid by using the 8-pt average of the surrounding grid cells. Figure 2 shows the normalized intensities for the 100 G snapshot for MURaM (left panel), RH based on the original MHD grid (middle panel), and RH using the 8-pt averaged grid (right panel). The RH intensities based on the original MHD grid show spurious fringes that disappear when using the 8-pt averaged MHD grid. We note that the decrease of the strong gradients when using interpolated atmospheres mostly results from a better sampling of the regions around the optical depth unity, where most of the radiation originates from. To confirm this point, we repeated the RH syntheses using snapshots interpolated on a vertical grid twice the resolution of the original one while keeping the original horizontal grid, and found that the intensity distribution is very similar to the one obtained on the 8-pt averaged snapshots. Overall, using interpolated grids reduces the peak intensity and increases the width of the intensity distribution, which is qualitatively what one would expect when increasing the spatial resolution.

Fig. 2. Relative intensity calculated with the MURaM (top left), COSI (top right), RH based on the original MHD grid (bottom left), and RH based on the 8-pt interpolated MHD grid (bottom right). Each snaphot was normalized to its respective mean intensity.

Figure 3 shows normalized intensity as calculated with the MURaM radiative transfer scheme (top panels) and the RH (bottom panels). From visual inspection it is not possible to identify any difference. In Fig. 4, top panels, we show the corresponding intensity distributions for all three snapshots. The purple lines show the intensity distribution of the RH spectral synthesis using the original MHD grid values, the red dashed lines the calculations with the 8-point interpolation of the grid cells, and the blue dotted line the MURaM intensities. While there are some small deviations, overall, the distributions show consistent results. Since the MHD cubes where computed with the MURaM RT scheme in the first place, and as the MHD grid interpolation removes spurious intensity fringes – while still producing consistent intensity distributions – we apply the 8-pt grid interpolation scheme for the RH and COSI calculations presented here.

Fig. 3. Relative intensity for the HD, 100 G and 200 G snapshots calculated from the MURaM radiation scheme (top row), the RH code (middle row) and COSI (bottom row). The intensities of all snapshots were normalized to the mean of the MURaM HD snapshot, which is 2.71 107 erg s−1 cm−2 nm−1 sr−1. For better visibility of the low-medium contrast features, the color scale only covers the range between 0.5 and above 1.5 and is kept constant outside that range.

To validate the code performance in more detail, in Fig. 4, bottom panels, we compare the distribution of the intensity calculated with the 8-pt interpolated snapshots for COSI (black line), RH (red dashed) and MURaM (blue dotted) obtained for the HD (left panel)), 100 G (middle panel), and 200 G (right panel). The comparison shows that the distributions cover the same range in absolute intensity. Furthermore, the overall shape of the distributions agrees very well, while some systematic differences can be identified. A key result of this comparison is that in the low-intensity wing of the distribution RH and MURaM agree very well, while COSI and RH reproduce the same intensities at the high-intensity wing of the distribution. This finding is systematically present in all three snapshots. As such, it can be concluded that RH generally produces slightly wider intensity distribution than MURaM and COSI. Comparing the COSI and MURaM calculations more closely, these two distributions appear systematically shifted, with COSI producing slightly brighter intensities than MURaM. Furthermore, taking into account the shift between COSI and MURaM, the features with more abundant intensity values at the peak of the distributions are consistently reproduced in all three codes. In summary, COSI, MURaM and RH reproduce consistently the same distribution envelope as well as details with small intensity variations for all snapshots.

Fig. 4. Top panel: intensity distribution function for the full 384 × 384 snapshot calculated with the RH code using the original grid in the MHD simulation (purple dot-dashed lines) and the 8-point interpolated MHD grid (red dashed lines) and the MURaM calculation (blue lines). Bottom panel: intensity distribution function for the HD (left panel), 100 G (middle panel) and 200 G (right panel) snapshots using the COSI (black lines), MURaM solver (blue lines), and RH (red lines) codes.

We are further interested in studying where the largest differences of the intensities come from. Figures 5 and 6 show the difference and the ratio, respectively, between emergent intensities (normalized to their average) obtained with MURaM, RH, and COSI. The images indicate that in the HD snapshots differences of roughly the same amplitude between the three codes are found in both dark intergranular lanes and bright upflow regions. Here, the pixel-to-pixel difference between the RH versus MURaM calculations ranges from –0.04 to 0.12, and the respective ratio from 0.96 to 1.10. For COSI and the MURaM the difference ranges from –0.24 to 0.46, and the ratio between 0.77 and 1.58.

Fig. 5. Top panels: relative pixel-to-pixel difference, (RHi – MURaMi)/< MURaM,HD>, for the HD, 100 G and 200 G snapshots. Bottom panels: same as top panels but for (COSIi – MURaMi)/< MURaM,HD>. For better visibility of the low-medium contrast features, the color scale only covers the given ranges between in the top and bottom panels and is kept constant outside of that range.

Fig. 6. Top panels: ratio of the intensities, RHi/MURaMi, with i referring to the HD, 100 G and 200 G snapshots, respectively. Bottom panels: same as top panel but for COSIi/MURaMi. For better visibility of the low-medium contrast features, the color scale only covers the given ranges between in the top and bottom panels and is kept constant outside of that range.

For the 100 G snapshot, the differences between RH and MURaM range from –0.08 to 0.09 and the ratio between 0.92 to 1.10, while for the case of COSI and MURaM the difference ranges from –0.58 to 0.63, and the ratio between 0.60 and 1.73. The difference and ratio for the 200 G simulation range from –0.05 to 0.16 and 0.95 to 1.13, respectively. For COSI and MURaM the difference ranges for the 200 G simulation ranges from –0.88 to 0.56, and the ratio between 0.55 and 1.71. The larger pixel-to-pixel variations of the COSI calculations compared to RH and MURaM stem from higher intensity values, mostly in the intergranular lanes, which will be further investigated. Comparing the averages of the intensities, the ratios of RH and MURaM give 1.012, 1.013, and 1.014 and for COSI and MURaM it is 1.042, 1.041, and 1.040 for the HD, 100 G, and 200 G snapshots, respectively. From this we conclude that while the pixel-to-pixel differences of the codes can vary by several percent. The averaged snapshots agree to 1.4% for the codes using the same opacities and to about 3–4% for the codes with different opacity sets.

Understanding how the intensity scales with the magnetic field strength is especially important in the context of irradiance studies, as some irradiance reconstruction models make use of the magnetic flux as proxy for the radiative intensity (e.g., Krivova et al. 2003; Foukal et al. 2011; Yeo et al. 2017). As such, the relation between brightness and magnetic flux has been the subject of several observational (e.g., Ortiz et al. 2002; Criscuoli et al. 2017) and theoretical studies (e.g., Röhrbein et al. 2011; Peck et al. 2019). Therefore, we are further interested in finding how well the intensities calculated by the different codes agree for different magnetic field strengths. Figure 7, top panels, shows the magnetic field strengths for each snapshot, interpolated to the height where τ = 1. The bottom panels show the respective absolute field strength for both snapshots for the same layer segmented at 500 G intervals. We then investigate how the normalized intensity depends on absolute magnetic field strength at τ = 1 for subsequent 100 G bins in both snapshots. Figure 8 compares the mean intensity for each of the bins for the 100 G and 200 G snapshots calculated with COSI (black), MURaM (blue) and RH (red). We note that the shape of these curves is heavily influenced by realization noise, since they are based on single MHD snapshots (see Rempel 2020 for a detailed discussion on the influence of realization noise on the emergent intensity). Consequently, while we cannot meaningfully compare the differences between the 100 and 200 G snapshots, we can compare the differences of the 3 radiative transfer codes for each of the setups, which is the goal of this work. Nonetheless, the overall shape of the dependence of the intensity as a function of magnetic field strength is in line with the findings by Rempel (2020). Specifically, the intensities in regions with small-medium magnetic field strength (100–900 G) are darker than in the non-magnetic areas. Second, the shape of the intensity variation as a function of magnetic field strength shows systematic differences for the 100 G and 200 G snapshots. In particular, the turning point of the intensity in the 100 G snapshot is found at around 500 G, while for the 200 G snapshot it is at around 350 G. This difference can be explained by the fact that the stronger magnetic field in the 200 G snapshot pushes the weak field, which is naturally confined to the intergranular lanes, towards the edges of granular cells and as such to areas with higher intensity. Third, while all three codes agree very well in their response to the varying plasma properties some detailed performance differences can still be identified. The scatter amongst the codes is systematically lower for the 100 G than for the 200 G snapshots. This goes in line with the differences found in Fig. 5, where the RH and MURaM code show the largest deviations in difference (top panels) and ratio (bottom panels) for the 200 G snapshot, and specifically in the intergranular lanes. Looking at the codes individually, in the 100 G snapshots all three codes agree within almost about 1–2% and the codes do not seem to produce systematic differences for that individual snapshot. This is somewhat different for the 200 G snapshot. Here the codes differ by about 2–3%. Moreover, the MURaM code tends to give systematically brighter intensity for the weak-medium magnetic field strengths. An indication of that can be also found in Fig. 4, top panel, where the MURaM calculations give consistently higher intensities at the low-intensity side of the peak of the distributions than the RH code. Finally, the average intensities for all magnetic field strength larger than 1000 G lead to a normalized intensity above unity.

Fig. 7. Top panel: magnetic field strengths at the height where τ = 1 for the 100 G and 200 G snapshots. Bottom panel: mask of the 100 G and 200 G snapshot to identify which pixels fall into 100 G intervals of the absolute magnetic field strength over the range for each snapshot.

Fig. 8. Normalized intensity as a function of absolute magnetic field strength Bτ = 1 for COSI (black), MURaM (blue), and RH (red). The diamonds give the normalized intensity determined from the 100 G simulation, while the stars the ones from the 200 G run. The last bin contains all magnetic field strength > 1000 G.

5. Conclusions

We have determined the emergent intensity from 3D MHD simulation snapshots for a non-magnetic case, and 100 G and 200 G simulation using the radiative transfer codes COSI, RH and the MURaM radiation scheme. We compared the difference, ratio and distribution of the intensities and find an overall good agreement amongst the codes. We find that while the absolute intensities produced by the codes agree very well, systematic differences on the code performance could be identified from the distribution functions. In particular, the RH code produces slightly wider distributions. The RH and MURaM calculations agree very well in the low-intensity range, while RH and COSI match very well on the higher intensities. While the pixel-to-pixel differences of the codes can vary by several percent, the averaged intensities for RH versus MURaM differ up to 1.4% and COSI versus RH and MURaM to about 3–4%, respectively.

Furthermore, we investigated how well the codes reproduce the intensities for different magnetic field strengths at the τ = 1 layer. Here, we find the differences between the codes to be about 1–2% for the 100 G snapshot and about 2–3% for the 200 G snapshot. The overall shape of the change in intensity as a function of magnetic field strength is in line with previous work. For the COSI and RH calculations we carried out an 8-point averaging of the cell corners of the MHD grid onto the cell centers as it is done in the MURaM code. In addition to ensuring consistency, this technique also removes fringes that appear when using the original MHD grid with RH and COSI. While we focused in this study on computing intensities, it is likely that this technique should also be considered for polarized radiative transfer as well. Tests performed using snapshots interpolated on the vertical grid suggest that MHD simulations with substantially higher resolution in the vertical τ-scale might not cause this issue. A detailed study investigating resolution dependence would need to follow.

Acknowledgments

We kindly acknowledge the support by the International Space Science Institute (ISSI), Bern, Switzerland during the International Team on Modeling Solar Irradiance with 3D MHD simulations, lead by Serena Criscuoli. MH kindly acknowledges support by Daniel Karbacher. TMDP’s research is supported by the Research Council of Norway through its Centres of Excellence scheme, project number 262622. The National Solar Observatory is operated by the Association of Universities for Research in Astronomy, Inc. (AURA) under cooperative agreement with the National Science Foundation. This material is based upon work supported by the National Center for Atmospheric Research, which is a major facility sponsored by the National Science Foundation under Cooperative Agreement No. 1852977.

References

Appendix A: Basic concept of radiative transfer

In the following we describe the basic concept of how the radative transfer is solved in the case of local thermodynamic equilibrium and then discuss the differences in the implementation in the case of MURaM, RH and COSI.

In local thermodynamic equilibrium (LTE) the velocity distribution of the atoms, ions and electrons is Maxwellian. Furthermore, the ionization rates and population numbers are a function of the local temperature and density. The assumption of LTE is only true if the collisional processes dominate over the radiative processes, e.g., for high densities, which is a suitable assumption for our study. The standard radiative transport equation in Eq. A.1

$$\begin{array}{c}\hfill \mu \frac{d{I}_{\nu \mu }}{\mathit{dz}}={\eta }_{\nu }-{\chi }_{\nu }{I}_{\nu \mu }\end{array}$$(A.1)

describes the change of the specific intensity I_ν at frequency ν along the path dz in a plane-parallel atmosphere due to the emission and absorption, where η is the emissivity and χ the extinction coefficient. In our study μ = cos Θ = 1 as we consider the radiative transfer in the vertical direction. Using the optical depth

$$\begin{array}{c}\hfill {\tau }_{\nu }={\displaystyle {\int }_0^{z^{\prime }}{\chi }_{\nu }dz,}\end{array}$$(A.2)

and the source function, defined as

$$\begin{array}{c}\hfill S_{\nu }\equiv \frac{{\eta }_{\nu }}{{\chi }_{\nu }},\end{array}$$(A.3)

the intensity change can be described as a function of the change of the optical depth as

$$\begin{array}{c}\hfill \frac{d{I}_{\nu }}{d\tau }=I(d{\tau }_{\nu })-S(d{\tau }_{\nu }).\end{array}$$(A.4)

Multiplying Eq. A.4 with e^−τ gives

$$\begin{array}{c}\hfill \frac{d{I}_{\nu }{e}^{-{\tau }_{\nu }}}{d{\tau }_{\nu }}=I_{\nu }{e}^{-\tau }-S_{\nu }{e}^{-{\tau }_{\nu }}\end{array}$$(A.5)

If S_ν is known, the integration of Eq. A.5 gives

$$\begin{array}{c}\hfill I_{\nu }({\tau }_1)=I_{\nu }({\tau }_2)e^{-({\tau }_2-{\tau }_1)}+{\displaystyle {\int }_{{\tau }_1}^{{\tau }_2}{S}_{\nu }(t)e^{-(t)}dt,}\end{array}$$(A.6)

It is assumed S_ν is a linear function of τ_ν. Then Eq. A.6 can be written as:

$$\begin{array}{c}\hfill I_{\nu }({\tau }_1)=I_{\nu }({\tau }_2)e^{-({\tau }_2-{\tau }_1)}+{\displaystyle {\int }_{{\tau }_1}^{{\tau }_2}{e}^{-t}[S_1+S_2(t-{\tau }_1)/({\tau }_2-{\tau }_1)]dt}\end{array}$$(A.7)

Clearly, the determination of the opacities and the source function has an impact in the emergent intensity. In the following, we give the basic concept of the numerical scheme of all three codes that is used to determine the emergent intensity in LTE and vertical direction (μ = 1).

Appendix B: MURaM

MURaM interpolates the opacity from a RH opacity table κ = f(ρ, T) using a bi-linear table interpolation. Starting with the values at the positions: z₁ (κ₁, 𝜚₁, S₁) and z₂ (κ₂, 𝜚₂, S₂) we compute the τ scale using Δτ = τ₂ − τ₁ > 0:

$$\begin{array}{c}\hfill \mathrm{\Delta }\tau =-\mathrm{\Delta }z[({\kappa }_1{\varrho }_1+{\kappa }_2{\varrho }_2)/3+({\kappa }_1{\varrho }_2+{\kappa }_2{\varrho }_1)/6]\end{array}$$(B.1)

The outgoing intensity at the top of the domain is computed solely for diagnostics using a scheme that is separate from the radiative transfer that is used to compute the intensity throughout the simulation domain (needed for radiative heating/cooling). The contribution to the outgoing intensity at the top of the domain from the interval [τ₁, τ₂] in regions with Δτ > 10⁻⁵ is given by:

$$\begin{array}{cc}\hfill c& =(1-e^{-\mathrm{\Delta }\tau })/\mathrm{\Delta }\tau \hfill \end{array}$$(B.2)

$$\begin{array}{cc}\hfill \mathrm{\Delta }I& =[S_1(1-c)+S_2(c-e^{-\mathrm{\Delta }\tau })]e^{-{\tau }_1}\hfill \end{array}$$(B.3)

otherwise:

$$\begin{array}{cc}\hfill \mathrm{\Delta }I& =0.5\mathrm{\Delta }\tau [S_1+S_2]e^{-{\tau }_1}\hfill \end{array}$$(B.4)

Appendix C: RH

In RH, the formal solution of Eq. (A.6), is achieved by piecewise integration using the formalism of characteristics (see Kunasz & Auer 1988). In this work we used the solver with linear interpolation of the source function between points τ₁ and τ₂, according to the expression:

$$\begin{array}{c}\hfill I_{\nu }({\tau }_1)=I_{\nu }({\tau }_2)e^{-\mathrm{\Delta }\tau }+w_0{S}_1+w_1\frac{\mathrm{\Delta }S}{\mathrm{\Delta }\tau },\end{array}$$(C.1)

where

$$\begin{array}{cc}\hfill \mathrm{\Delta }\tau & =0.5\mathrm{\Delta }z({\chi }_1+{\chi }_2),\hfill \\ \hfill \mathrm{\Delta }S& =S_2-S_1,\hfill \\ \hfill w_0& =1-e^{-\mathrm{\Delta }\tau },\hfill \\ \hfill w_1& =w_0-\mathrm{\Delta }\tau e^{-\mathrm{\Delta }\tau }.\hfill \end{array}$$

To preserve numerical accuracy, when Δτ > 50, w₀ = w₁ = 1 and when Δτ < 5 ⋅ 10⁻⁴, w₀ and w₁ are approximated by

$$\begin{array}{cc}\hfill w_0& =\mathrm{\Delta }\tau (1-\frac{\mathrm{\Delta }\tau }{2}),\hfill \\ \hfill w_1& =\mathrm{\Delta }{\tau }^2(\frac{1}{2}-\frac{\mathrm{\Delta }\tau }{3}).\hfill \end{array}$$

In RH the extinction coefficients are explicitly calculated, while in MURaM they are interpolated from the RH opacity table (κ ≡ χ/ρ = f(ρ, T)). For consistency, both the explicit RH synthesis and the opacity table employed in MURaM were computed under LTE using the same set of input parameters. In summary, the synthesis included the full list of atomic and molecular bound-free transitions from the Kurucz website, photo-ionization of 12 atomic species (including updated Fe atom cross-sections) and photo-dissociation of 52 diatomic molecules. The computation took also into account Thomson and Rayleigh scattering, although scattering contribution is expected to be negligible in continua along the vertical line-of-sight (e.g., Fabbian & Moreno-Insertis 2015).

Appendix D: COSI

The numerical scheme of COSI consists of two modules. The first, also referred to as hminus (Haberreiter et al. 2008b; Shapiro et al. 2010) calculates the level populations of the atoms, which under the assumption of LTE follow the Boltzmann distribution

$$\begin{array}{c}\hfill \frac{n_{\mathrm{u}}}{n_{\mathrm{l}}}=\frac{g_{\mathrm{u}}}{g_{\mathrm{l}}}{e}^{-({\chi }_{\mathrm{u}}-{\chi }_{\mathrm{l}})/kT}=\frac{g_u}{g_l}{e}^{-(h\nu /kT)}.\end{array}$$(D.1)

The second module of the code, also referred to as fioss, uses the level populations to calculate the emergent intensity, in our case in the vertical direction (μ = 1). In COSI the opacities are determined taking into account all radiative (negligible in the case of LTE) and collisional processes for absorption, emission and ionization. In COSI the solution of Eq. (A.6) is obtained via consideration of the change of the intensity between 3 subsequent depth points, n = 1,2,3 and τ₃ > τ₂ > τ₁:

if τ > 10⁻¹⁰ then

$$\begin{array}{cc}\hfill I_{\nu }& =I_{\nu }{e}^{-\mathrm{\Delta }\tau }+d{S}_{\nu }/d\tau =\hfill \end{array}$$(D.2)

$$\begin{array}{cc}\hfill & =I_{\nu }{e}^{-\mathrm{\Delta }\tau }+S_{\nu }{w}_{\tau }\hfill \end{array}$$(D.3)

with

$$\begin{array}{cc}\hfill w_{\tau }& =w_a+w_b\hfill \end{array}$$(D.4)

$$\begin{array}{cc}\hfill w_a& =(e^{-{\tau }_1}-e^{-{\tau }_2})/\mathrm{\Delta }{\tau }_1\hfill \end{array}$$(D.5)

$$\begin{array}{cc}\hfill w_b& =(e^{-{\tau }_3}-e^{-{\tau }_2})/\mathrm{\Delta }{\tau }_2\hfill \end{array}$$(D.6)

$$\begin{array}{cc}\hfill \mathrm{\Delta }{\tau }_1& =0.5\mathrm{\Delta }{z}_1({\chi }_1+{\chi }_2)\hfill \end{array}$$(D.7)

$$\begin{array}{cc}\hfill \mathrm{\Delta }{\tau }_2& =0.5\mathrm{\Delta }{z}_2({\chi }_2+{\chi }_3).\hfill \end{array}$$(D.8)

We note that the S₂e^−τ₂ term (from equation ) cancels for consecutive depth points with the next term, which is then equal to the new S₁-term.

All Figures

Fig. 1. Top panels: photospheric magnetic field strength Bz, τ = 1 for the snapshots of the pure hydrodynamic case (left panel), and the cases with a magnetic field strength of 100 G (middle panel), and 200 G (right panel); bottom panels: horizontal temperature variation of the snapshots at the top of the simulation box for the pure hydrodynamical case (left panel), and the cases with a magnetic field strength of 100 G (middle panel), and 200 G (right panel).

In the text

	Fig. 2. Relative intensity calculated with the MURaM (top left), COSI (top right), RH based on the original MHD grid (bottom left), and RH based on the 8-pt interpolated MHD grid (bottom right). Each snaphot was normalized to its respective mean intensity.
In the text

Fig. 3. Relative intensity for the HD, 100 G and 200 G snapshots calculated from the MURaM radiation scheme (top row), the RH code (middle row) and COSI (bottom row). The intensities of all snapshots were normalized to the mean of the MURaM HD snapshot, which is 2.71 107 erg s−1 cm−2 nm−1 sr−1. For better visibility of the low-medium contrast features, the color scale only covers the range between 0.5 and above 1.5 and is kept constant outside that range.

In the text

Fig. 4. Top panel: intensity distribution function for the full 384 × 384 snapshot calculated with the RH code using the original grid in the MHD simulation (purple dot-dashed lines) and the 8-point interpolated MHD grid (red dashed lines) and the MURaM calculation (blue lines). Bottom panel: intensity distribution function for the HD (left panel), 100 G (middle panel) and 200 G (right panel) snapshots using the COSI (black lines), MURaM solver (blue lines), and RH (red lines) codes.

In the text

	Fig. 5. Top panels: relative pixel-to-pixel difference, (RHi – MURaMi)/< MURaM,HD>, for the HD, 100 G and 200 G snapshots. Bottom panels: same as top panels but for (COSIi – MURaMi)/< MURaM,HD>. For better visibility of the low-medium contrast features, the color scale only covers the given ranges between in the top and bottom panels and is kept constant outside of that range.
In the text

	Fig. 6. Top panels: ratio of the intensities, RHi/MURaMi, with i referring to the HD, 100 G and 200 G snapshots, respectively. Bottom panels: same as top panel but for COSIi/MURaMi. For better visibility of the low-medium contrast features, the color scale only covers the given ranges between in the top and bottom panels and is kept constant outside of that range.
In the text

	Fig. 7. Top panel: magnetic field strengths at the height where τ = 1 for the 100 G and 200 G snapshots. Bottom panel: mask of the 100 G and 200 G snapshot to identify which pixels fall into 100 G intervals of the absolute magnetic field strength over the range for each snapshot.
In the text

	Fig. 8. Normalized intensity as a function of absolute magnetic field strength Bτ = 1 for COSI (black), MURaM (blue), and RH (red). The diamonds give the normalized intensity determined from the 100 G simulation, while the stars the ones from the 200 G run. The last bin contains all magnetic field strength > 1000 G.
In the text