© The Authors 2025

1. Introduction

The new generation of solar telescopes, such as the Daniel K. Inouye Solar Telescope (DKIST, Rimmele et al. 2020) in Hawaii or the upcoming European Solar Telescope (EST, Quintero Noda et al. 2022) on the Canary Islands, help us observe the dynamics in the solar atmosphere at unprecedented spatial, spectral, and temporal scales. In particular, these telescopes can in principle resolve structures well below the horizontal photon mean free path of 50–100 km in the solar photosphere under the assumption that the solar atmosphere is plane-parallel (Hinode Review Team et al. 2019). However, observations taken at the diffraction limits of these telescopes may capture the imprints of horizontal radiative transfer (RT) or the three-dimensional (3D) RT effects.

Atmospheric structures (in density and temperature) at spatial scales well below the typical photon mean free path, λ, may be produced by magneto-hydrodynamic (MHD) processes, in particular by turbulence and magnetoturbulence. However, whether these structures can be observed has been a matter of debate. A series of investigations with two-dimensional (2D) idealized arrangements of the atmospheric structures have been performed (see, e.g., Stenholm & Stenflo 1977, 1978; Kneer 1981; Bruls & von der Lühe 2001; Holzreuter & Solanki 2012, and references therein). While these papers all considered the Fe I lines at 525 nm, Kneer (1981) looked into the horizontal transfer effects in the solar chromosphere. These investigations have shown that the horizontal RT in combination with strong line scattering may lead to a spatial smearing of the atmospheric structures, and therefore to a reduced contrast in the observations. The author also found that the spatial smearing due to horizontal RT is insignificant in the photosphere, and can be more dominant at chromospheric heights where lines are strongly scattering.

More recently, Judge et al. (2015) reconciled the different results by finding that the amount of scattering in a specific atomic line is the main responsible ingredient determining the degree of spatial smearing due to NLTE effects, especially in the chromosphere. The authors predicted that this type of smearing results in a solar ‘fog’ whose effects are comparable to that of the point spread function (PSF) of a telescope. The chromospheric spectral lines are well known to be affected by the 3D non-local thermodynamic equilibrium (NLTE) effects (for e.g., Leenaarts & Carlsson 2012; Bjørgen et al. 2018, 2019; Judge et al. 2020). However, to our knowledge no study has been done that tests how clearly the new telescopes will be able to see structures at scales down to 20–30 km in the solar photosphere using high-resolution state-of-the-art 3D MHD simulations. The question is particularly acute for the magnetic field, which is typically measured in lines of Fe I that are affected by departures from local thermodynamic equilibrium (LTE) (e.g. Solanki & Steenbock 1988; Smitha et al. 2020, 2023).

In our previous work (Holzreuter & Solanki 2012, 2013, 2015, papers I, II, and, III hereafter) on photospheric Fe I lines, we found that many atmospheric structures smaller than λ remained observable, confirming earlier investigations (e.g. Bruls & von der Lühe 2001) on that topic. In papers II and III, in which we used radiation-MHD simulations to model the solar photosphere, we also found that the results of the idealized models of flux tubes (FTs) and flux sheets (FSs) do not correspond well to the ones obtained using the much more realistic MHD simulations. The main shortcoming of these earlier studies was that the MHD simulations employed had a grid spacing of 20 km in the horizontal direction, which was too large to address the kind of resolution that is expected to be achieved by the new telescopes. Therefore, in this work we extend these investigations using a 3D atmospheric input with a much finer grid spacing to check whether the principal results of papers II and III continue to hold.

2. Model ingredients

2.1. Model atmosphere

The input model atmosphere for our calculations was taken from a 3D radiation MHD simulation calculated with the MURaM code (Vögler et al. 2005). The cube employed was the same as the one used in Smitha et al. (2020, 2021). The cube was composed of 1024  ×  1024  ×  256 voxels spanning a geometrical range of approximately 6 Mm  ×  6 Mm  ×  1.5 Mm. The initial seed field was homogeneous and vertical, with mixed polarity in a 2  ×  2 checkerboard pattern; that is, the upper left and lower right quadrants of the whole domain were filled with a negative homogeneous vertical field of 200 G each, and the other two with a positive one. This simulation, with periodic horizontal boundaries and a vertical field at the top boundary, with inflows and outflows allowed but no passage of magnetic flux at the bottom boundary, was allowed to run for nearly 30 minutes of solar time after the introduction of the initial field before a snapshot was chosen for analysis. The magnetic flux was thus conserved during the simulations.

The 2  ×  2 checkerboard pattern is still visible in the final snapshot from the map of its vertical magnetic field at log(τ_c⁵⁰⁰) = 0.0 shown in Figure 1. When the map is divided into quadrants, the top left and bottom right parts carry mostly negative polarity fields, while in the remaining two parts those with positive polarity prevail; hence, the 2  ×  2 checkerboard pattern. For the RT calculations, the optically thick part at the bottom was cut as in paper III, ending up with a 1024 × 1024  ×  150 data cube spanning a height range from approximately 430 km below the average τ_c⁵⁰⁰ = 1 (hereafter referred to as τ_c) level to 740 km above it. The smallest τ_c value at the bottom of the cube was ≈38, while the average value at that depth amounts to almost 1000.

Fig. 1. Map of the vertical component of the magnetic field at the surface (log(τc500) = 0.0).

Figure 2 displays the continuum intensity, I_c, (close in wavelength to 525 nm) from the snapshot chosen for analysis. As the available computing resources were not large enough to make a true 3D NLTE calculation with the full atmospheric model, we had to choose between lowering the resolution for the 3D RT (as has typically been done in the literature) and restricting the spectral synthesis to only a part of the data cube. As the high resolution of the simulation is central to the aims of this paper, we selected a representative sub-domain – containing FT- and FS-like structures – for which we performed a full 3D NLTE calculation. Figure 2 indicates the location of these two sub-sections of the full box in which 3D RT was carried out (yellow boxes).

Fig. 2. Continuum intensity at a wavelength close to the 525 nm lines. The squares indicate the two selected sub-domains that were calculated in full 3D NLTE. The bigger dashed squares indicate the domain used for calculations, while the smaller solid squares within the dashed squares represent the areas that were used for our investigation.

Some important atmospheric parameters for the selected sub-domain at the same geometric heights as in papers II and III – z = 0 km (average height of τ_c = 1) and z = 240 km (approximate average formation height of the cores of the considered lines) – are given in Figure 3. The average height of τ_c = 1 in the selected sub-domain is equal to that for the whole atmosphere. The selected sub-domain is located at the boundary of three granules. At z = 0, each of the three intergranular lanes contains strong magnetic FSs of up to 1.8 kG. The density increases towards the boundary of the granules, forming thick and relatively cool downflowing walls close to the flux elements. In the centre of the magnetic elements the downflows are often reduced or even inverted. The lower right intergranular lane even shows strong upflows in the centre of the FS. At z = 240 km, the situation is much more complex when the temperature and velocities are considered. Especially in the intergranular lanes, which are much broader at that height, the atmospheric parameters vary strongly on a very short scale of a few tens of kilometres.

Fig. 3. Main atmospheric input parameters of the model sub-domain used in this work at heights z = 0 km (i.e. at the height of average τc = 1; panels on the right) and z = 240 km (the approximate line core formation region of the investigated lines; panels on the left). Plotted, from top to bottom, are the vertical component of the magnetic field, Bz, the hydrogen number density, nH, the vertical component of the velocity, vz, and the temperature, T.

2.2. RT calculations

As in the previous papers of this series, we used the code RH (Uitenbroek 2000, 2001) with the adaptions mentioned in papers I and II. The same atomic model (23 levels, 33 lines, and ≈1300 wavelength points) was used as in paper III. To investigate the influence of horizontal RT and NLTE effects separately, we again calculated the line spectra in LTE, 1D NLTE, and true 3D NLTE. For details, we refer the reader to these papers.

Cutting a domain from the full atmospheric box violates the assumption of periodic boundary conditions in the 3D version of the RH code; in other words, at the boundaries of the sub-models the atmospheric quantities are not continuous anymore. At a point close to the boundary, the local state of the atoms can then erroneously be influenced by a region close to the opposite boundary; that is, from a part of the atmosphere reasonably far away, definitely further than the horizontal photon mean free path. To investigate the influence of the discontinuity, we selected two sub-domains, each spanning 256  ×  256 grid points but one shifted by 64 grid points along both the x and y co-ordinate against the other. For both sub-domains, a full 3D NLTE calculation was performed. The two sub-domains selected are marked by the dashed squares in Figure 2.

We chose the equivalent width (EW) of the 630.15 nm line to check the influence of the non-periodic boundary conditions on the results of the 3D RT at the boundaries. The EW calculated in any two 3D runs should be equal if the effects at the non-periodic boundaries are negligible. By visual inspection of the profiles at the boundaries of each sub-domain, we made sure that they have a similar shape in both runs, ruling out the possibility that the EW values are not equal by chance due to differently shaped profiles.

Figure 4 presents the relative differences of the EW values for the 630.15 nm line in pixels common to both sub-domains, each calculated with one of the 3D runs (δE_{3D − 3D}) where

Fig. 4. Spatial distribution of the relative differences between the EW values of the Fe I 630.15 nm spectral line calculated in the two spatially shifted sub-domains in 3D NLTE. The colours denote the deviation in the EW between the two runs. The rectangle indicates the area chosen for our investigation, which lies completely inside the region with less than 1% difference.

$$\begin{array}{c}\hfill \delta E_{M_i-M_j}=\frac{{\mathrm{EW}}^{M_i}-{\mathrm{EW}}^{M_j}}{{\mathrm{EW}}^{M_j}},\end{array}$$(1)

with M_i and M_j each referring to one of the three RT calculation methods (3D, 1D, LTE), or, in the above case, the two spatially shifted 3D runs.

As was expected, the effects are largest close to the boundary and decrease inwards. At the boundary itself, at some selected locations, the difference between the two EWs is as high as 30%. It decreases very rapidly towards the inner part of the common spatial area. By cutting 32 points, corresponding to slightly less than 200 km at each boundary, we have an agreement in EW that is clearly better than one percent in the whole area. In the following, all of the results were derived in the area left after the removal of these 32 outer spatial points.

3. Results

In paper III, the horizontal size of a voxel in the MHD input atmosphere was roughly 21 km. The cube used in this investigation has a horizontal voxel size of slightly below 6 km, surpassing all 3D NLTE RT calculations performed so far (to our knowledge). This enables us to investigate effects of horizontal RT at a spatial scale clearly below the resolution limit of the new generation of telescopes. In the following, we investigate the strength of horizontal RT effects on spatial scales smaller than the horizontal photon path of roughly 100 km in the photosphere.

3.1. Horizontal radiation transport in a highly resolved flux element

As in papers I and III of this series, we first selected a cut through a flux element, along which the influence of 1D NLTE and horizontal RT (3D NLTE) was studied. Figure 5 presents the spatial distribution of the continuum intensity, I_c, in the part of the atmosphere where the 3D NLTE calculation was performed (see also the upper right solid square in Figures 2, and 3). The yellow line indicates the position of the cut, which extends from one granular element to another, cutting a relatively small (≈80 km diameter) but typical FS containing a strong vertical magnetic field of B_z(τ_c = 1) > 1.5 kG.

Fig. 5. Enlarged section (≈1 Mm  ×  1 Mm) of the continuum intensity image (Figure 2) showing the part of the atmosphere where the 3D NLTE calculation was performed. The yellow line and the five yellow dots (a–e) indicate selected positions used in the text for further analysis. The yellow dots were drawn larger (2  ×  2) than in the actual image for better visibility. The black contours indicate the boundary of the magnetic element corresponding roughly to Bz(z = 0) = 1200 G.

Figure 6 shows different quantities along the cut at the height of z = 0 km (the average τ_c = 1 level). In the top panel, the temperature and density profiles along the cut are shown. While the density increases slowly as we approach the FS and abruptly falls within the magnetic part of the FS, the lowering of the temperature is more continuous and much broader than the flux element itself. On the basis of the temperature curve alone, the spatial location of the FS could not even be identified. This is in good agreement with the scenario found in paper III for MURaM simulations with a lower resolution, but does not conform to the ones of paper I and Stenholm & Stenflo (1977, 1978), Bruls & von der Lühe (2001), in which idealized models with a strong temperature gradient at the FS-FT boundary were used. This demonstrates that magnetic diagnostics are much better suited to identifying structures than the intensity.

Fig. 6. Selected quantities along a cut through a FS indicated in Figure 5. Top panel: Temperature (black) and density (blue). Middle panel: Vertical magnetic field (black) and vertical velocity (blue, negative values denote downflows). Bottom panel: δE1D − 3D (black) and continuum intensity, Ic (blue). All the quantities were taken at a geometrical height of z = 0 km. The vertical dotted lines indicate the boundary region of the flux element over which the density decreases and the magnetic field increases.

The middle panel presents the variation in the vertical magnetic field, B_z, and the vertical velocity, v_z (positive values denote upflows). In the FS, B_z reaches approximately 1.5 kG, whereas in the granules it typically lies below 100 G (although it reaches nearly 300 G at one location). In the small region with intense downflows just outside the FS, one observes a field enhancement of up to 500 G but with opposite polarity. The opposite polarity flux is associated with the kilo-Gauss FT. It can be brought to the surface through flux emergence due to turbulent motions (Chitta et al. 2019). Another possibility is that the opposite polarity flux is brought down from the canopy by the downflows (Pietarila et al. 2011). Both of these are due to small-scale intergranular lane turbulence distorting the field lines at the edge of the tubes. This corresponds to the opposite polarity fields deduced by Bühler et al. (2015) from Hinode data. As was expected, the vertical velocity points upwards in the granules and strongly downwards close to the flux element. In the inner parts of the flux element, v_z is much smaller. The strong downflows are located in areas of weak or vanishing B_z. In the inner parts of the selected FS, the velocity is not homogeneous. At some x locations of the FS, weak to medium-strength downflows occur, while at others, upflows may be observed (see also Figure 3).

In the lowest panel of Figure 6, the EW ratio, δE_{1D − 3D}, is plotted. It serves as a measure of line weakening (δE > 0) or strengthening (δE < 0) because of horizontal RT. The line weakening observed in the 3D NLTE calculation as described for the first time by Stenholm & Stenflo (1977) is only present in the boundary area outside the magnetic element. Thanks to the higher resolution of the simulations used here, it is even more obvious than in paper III that there is no line weakening in the inner parts of the FS. Contrary to the findings from idealized FT models (e.g. Stenholm & Stenflo 1977; Bruls & von der Lühe 2001; paper I), a distinct line strengthening is found in the centre of the flux element due to the cooler surroundings (see paper III and Smitha et al. 2020). Although along the cut, the temperature increases away from the centre of the FS (top panel in Figure 6), the surroundings of the FS are on average cooler compared to the centre, as is seen in Figure 3. Due to horizontal transfer, these cooler surroundings result in an overall weakening of the line profile. The blue curve in the lowest panel of Figure 6 represents the continuum intensity, I_c, at ≈525 nm. The variation in I_c is the opposite of that in δE_{1D − 3D}, confirming the findings of papers II and III, in which we found that in areas with higher I_c, the horizontal effects tend to strengthen the spectral lines, whereas in areas with lower I_c, horizontal transfer weakens the spectral lines.

3.2. Resolution of walls: still an open question

The thickness and the temperature profile of the wall of a magnetic element are central factors determining the line weakening or strengthening of Fe I lines calculated in 3D NLTE (see also paper I, in which we investigated the influence of wall thickness). An open question in paper III was whether a better spatial resolution would give rise to a different picture of this situation. There, the boundary of a flux element – defined by the drop in density from 90% to 10% – extended over two to three voxels, corresponding to a spatial scale of approximately 50 km. With our current high-resolution atmosphere, the thickness of the wall still extends approximately over three voxels, corresponding now to less than 20 km. Whether the wall is fully resolved or not still remains unanswered. Presumably, the thickness of three voxels is limited by the MHD modelling rather than by any other influencing factors. It is quite possible that with an even higher resolution in the MHD calculation, the wall would be even thinner. One could be tempted to believe that a thinner wall such as this could then decrease the separation of the cooler inner parts of the FS from the outer hotter surroundings, as is the case in FT cartoon models, and therefore increase the UV irradiation into the flux element. However, if we look at Figure 6, we find that the temperature gradient is not determined by either the density drop or the field strength enhancement. It even spreads far beyond the width of the downflow area outside the wall. We conclude that, even when the wall is not resolved, our results for the temperature profile, and therefore the observed horizontal RT effects, will not change qualitatively if the resolution is increased further. However, quantitative changes are of course possible.

The downflow area is much broader than three voxels, and therefore likely to be fully resolved. In paper III, the downflows spread over the whole flux element and beyond, whereas they are now clearly concentrated in the immediate surroundings of the magnetic element. Thus, we can conclude that the fundamental picture is not likely to change with a further increase in the spatial resolution of the simulation box. Owing to the better localization of the wall and especially the downflows, one could speculate that the strength of horizontal RT effects should increase with increasing resolution, at least to some extent. However, if we compare the strength of the line weakening or strengthening in the present atmosphere with that of paper III, the effects are even slightly smaller in the present, high-resolution case.

3.3. Variation in transfer effects on smaller scales

We now investigate the spatial scales on which 1DNLTE and horizontal RT effects may vary. Whether the photon mean free path, λ, could serve as a lower limit for the scales on which the 3D NLTE effects can be observed has been debated for a long time. Based on 3D-NLTE calculations of the chromospheric Ca I 422.7 nm line, Judge et al. (2015) argued that the amount of scattering in an observed line is the main factor controlling the extent of smearing due to horizontal RT effects. Since the photospheric lines are in general less scattering than the chromospheric lines, it should in principle be possible to find variations at the photosphere on spatial scales below λ because the temperature, the density, and the magnetic field, respectively, vary on scales much smaller than λ (e.g. Schüssler 1986, see also the width of the FS walls in Fig. 6).

We show in this section that 3D NLTE effects in the iron lines considered here may vary on scales much smaller than λ. We selected five locations in our atmosphere represented by the yellow dots in Figure 5. From the bottom to the top: (a) is located in the middle of a granule, far away from any magnetic elements, (b) lies in a downflow region with low I_c adjacent to a FS, (c) is situated at the boundary of the FS, between bright and dark I_c regions, (d) lies in the middle of the FS (high I_c), and (e), again, is located at the boundary but on the other side of the FS. We note that locations (b) to (e) lie close together, each less than 25 km from its neighbours.

Figure 7 shows (typical) line profiles of the 630.15 nm and the 525.02 nm lines at the selected positions. In the granule (a), the profiles for 1D NLTE and 3D NLTE virtually coincide. No major effects due to horizontal RT occur, as was expected from the relatively homogeneous atmospheric conditions. The LTE profiles differ slightly from their NLTE counterparts, the LTE line core intensity in the 630.15 nm line being slightly higher, the one in the 525.02 nm line slightly lower. At position (b) in the downflow area close to the FS, the 3D NLTE lines are weakened relative to their 1D NLTE counterparts owing to strong UV irradiation from the hot area in the granule at the bottom. We note that this location corresponds best to the situation given in the cartoon models of Stenholm & Stenflo (1977, 1978), Bruls & von der Lühe (2001) and paper I. At position (c), however, the differences between the 3D and 1D NLTE line intensities disappear. The place marks the reversal point of the effects between locations (b) and (d), respectively. In the middle of the FS, at location (d) a considerable strengthening of the 3-d NLTE lines is observed (as also found in paper III). Position (e) at the other boundary is much more complex as it lies close to the FS boundary where the FS folds on itself. Hence, the profiles at (e) are influenced by a complex mix of both hotter and cooler surroundings. The line core intensities in both the 630.15 nm line and the 525.02 nm line are weakened by the horizontal transfer effects. We note that the 630.15 nm line is affected by both scattering as well as UV overionization effects, while for the 525.02 nm line, UV overionization is the dominant NLTE mechanism (see paper II, III, Smitha et al. 2020, 2021, 2023). Hence, at points (a)–(e), the 1D NLTE profile of the 630.15 nm line can be either stronger or weaker than the LTE profile, whereas for the 525.02 nm line the 1D NLTE profile is always weaker than the LTE profile.

Fig. 7. Intensity profiles of the 630.15 nm (upper row) and the 525.02 nm (lower row) lines at five selected spatial positions indicated in Figure 5. Ordering: (a) location with lowest y co-ordinate (granule) to (e) location with largest y co-ordinate (at the upper boundary of the FS). Legends explaining the colours of the curves are given in panel a.

These sample profiles confirm the findings of paper III. Furthermore, they show that variations due to horizontal RT may be observed on rather small spatial scales. The 3D NLTE line strengthening at position (d) and the line weakening at position (e) lie only 4 voxels apart from each other; that is, the spatial distance between (d) and (e) is less than 25 km. Similarly, locations (b) and (d) lie 45 km apart. These are typical distances at which line strengthening reverts to line weakening, as can be seen also from the δE^{1D − 3D} curve in the bottom panel of Figure 6, which displays the effect of horizontal RT along the cut indicated in Fig. 5. There, δE^{1D − 3D} also varies strongly within approximately 50 km. Over a given region, the larger the local horizontal gradients in temperature and intensity, the stronger the variations in profiles due to horizontal RT (Carlsson et al. 2004).

Figure 8 shows two sample profiles of Stokes V/I_c at locations (d) and (e). The same effect as in the I profiles (see Figure 7) can be found. The strength of the NLTE and horizontal RT effects is of a similar magnitude as in Stokes I, a finding that has been stated already in Stenholm & Stenflo (1978) based on a simple FT cartoon model. In some cases, we observe that the V/I_c maxima are slightly shifted in wavelength relative to the Stokes I profile. This is most probably owing to the height dependence of the vertical velocity, as the Stokes V/I_c profiles are formed deeper in the atmosphere than the Stokes I line core. Therefore, the Stokes V/I_c profiles probe a slightly different height regime.

Fig. 8. Sample profiles of Stokes V/Ic of the 525.02 nm line at the two locations, (d) and (e). See Figures 5 and 7 for an explanation of the locations.

3.4. Contrasts and spatial smearing

The horizontal RT can affect the quality of the observed images in two ways: one, by reducing the image contrasts (Holzreuter & Solanki 2013, 2015), and two, by smearing out the spatial structures (Judge et al. 2020). In the sections below, we discuss these two effects, one after the other.

3.4.1. Effect on contrasts

In order to understand effects of horizontal transfer on the brightness contrast in the line cores, we investigated how strongly the RMS contrast changes depending on which of the three different computation methods is used (LTE, 1D NLTE, and 3D NLTE). The RMS contrast is defined as the standard deviation divided by the mean, expressed as a percentage (Danilovic et al. 2008). In every case, the RMS contrasts were calculated at the rest wavelengths of the lines.

A large part of the spatial region chosen for the 3D NLTE computations (Figure 5) is covered by granules that appear dark in the line core and that have a low RMS contrast. As the contrast depends on the area selected, in other words the structures contained and their sizes, we investigated the influence of the computation methods on different geometrical scales. To better quantify the effect of 3D NLTE on contrasts, we chose a smaller region approximately corresponding in size to the mean free photon path, λ. The region contains a part of a magnetic element harbouring strong density variations. This region is shown in Figure 9 for all lines and computation methods, with the RMS contrast in each case indicated in the respective panel of Figure 9.

Fig. 9. Comparison of RMS contrasts over a small region of the magnetic element for the three computation methods: LTE (first column), 1D NLTE (second column), and 3D NLTE (third column). In each row, we show the intensity images at the rest wavelength of the respective line, indicated in blue. The colour scale is the same as in Figure 10.

Between the 525 nm line pair and the 630 nm line pair, the contrast is largest for the 524.71 nm line because it has a lower excitation potential and is formed deeper in the atmosphere compared to the lines in the 630 nm pair. For the 525 nm pair, the change in contrast going from LTE to 1D NLTE is quite small, while the 630 nm line pair, interestingly, shows an increase in contrast. Figures 5 and 6 of paper II demonstrate that the contrast for the lines in 630 nm pair can either decrease or increase going from LTE to 1D NLTE, depending on whether the contrast is computed at the minimum of the line core intensity or at the rest wavelength of the line. On the other hand, for the 525 nm pair, the contrast decreases in both cases. Since the contrasts in Figure 9 are computed at the rest wavelength of the lines, we see an increase in contrast going from LTE to 1D NLTE for the 630 nm pair. On the other hand, the 3D NLTE computations in all cases show reduced contrasts relative to 1D NLTE, with the 630.1 nm line showing the largest reduction. This decrease in contrast in the 3D NLTE computations can be attributed to the horizontal transfer effects.

A similar analysis was carried out in paper II for the purely hydrodynamic case. There, the difference in contrast was larger between LTE and 1D NLTE than between 3D NLTE and 1D NLTE. This could be due not only to differences in the simulations but also the actual values of the contrasts being dependent on the area considered and the type of features it contains. The is supported by the histogram plotted in Figure 10, in which we divided the entire region into small squares of 5 pixels × 5 pixels. In each square, the contrast was determined for the intensities computed in LTE, 1D NLTE, and 3D NLTE. Then, a histogram was created for regions that show a contrast ≥10%. Only the 630.15 nm line was chosen for this analysis, since the contrast in this line showed the strongest response to 3D effects compared to other lines, as is seen from Figure 9. This is because, of all the four lines considered, the 630.15 nm is most affected by scattering (see paper III, Smitha et al. 2020), and the stronger the scattering effects in a given line, the greater the reduction in contrast due to horizontal transfer effects, according to Judge et al. (2015).

Fig. 10. Histogram of RMS contrast of small areas: we divided the full area of the intensity image into a grid of squares, each of which had a size of 5 pixels × 5 pixels. The grid is indicated as orange dots on the intensity image shown in the left panel. A histogram of the RMS contrasts computed at every square in the grid is shown in the right panel, for the three computation methods: LTE (black), 1D NLTE (red), and 3D NLTE (blue). The x axis has been clipped to show only bins with contrasts ≥10%.

From Figure 10, we can learn how the computation method influences the contrast on very small scales (below 25 km). The general trend of a contrast reduction by horizontal RT, in other words by 3D NLTE computations, is clearly seen from the histogram. In some regions, the contrast from 1D NLTE reaches values as high as 50%, but for 3D NLTE the contrast values do not exceed 37%. The number of areas showing large contrasts is also strongly reduced for 3D NLTE. For example, the number of areas showing contrast values of ≥25% has more than halved, from approximately 40 (1D NLTE) to less than 20 (3D NLTE). Overall, on very small scales (here, below 25 km), horizontal transfer effects may degrade the contrast, but on a more global scale the influence is much smaller, depending of course on the type of feature considered. As magnetic structures show very strong intensity variations on even smaller scales, it is still possible to resolve the structures – given the necessary telescope resolution – though with reduced contrast. This can also be seen from Figure 9.

3.4.2. Spatial smearing

In addition to reduction in image contrast, horizontal transfer effects spatially smear structures in the solar atmosphere. Signatures of spatial smearing due to 3D effects are evident by comparing, for example, the 1D NLTE and 3D NLTE maps for the 630.1 nm line in Figure 9, in which we present a blow-up of a small region around the FT. In observations, spatial smearing is also introduced by the diffraction limit of the telescope. To compare the two effects, we spatially degraded our computations. Once again, we choose the 630.1 nm line, since it is more affected by the 3D NLTE effects than the other lines considered in this study.

The intensity images from both 1D NLTE and 3D NLTE computations were spatially degraded to match the specifications of the Visible Spectro-polarimeter (ViSP; de Wijn et al. 2022) at DKIST. The degradation was done in two steps: first, the intensity was convolved with a Gaussian, mimicking a PSF that has a FWHM of 0.07 arcsec corresponding to the theoretical spatial resolution of ViSP at 630 nm¹. The convolved profiles were then spatially re-binned to the detector pixel size. No spectral degradation was applied.

In Figure 11, we present the intensity images for the 630.1 nm line at its rest wavelength, from both 1D NLTE and 3D NLTE computations before and after the two-step spatial degradation. The contrast decreases from 1D NLTE to 3D NLTE computations irrespective of whether instrumental degradation is applied or not (see Figure 11). Surprisingly, the decrease in the global RMS contrast due to the PSF of the telescope (5–6%) is comparable to that caused from horizontal transfer in 3D NLTE calculations (5–6%). Convolution with the PSF of DKIST causes a much stronger spatial smearing than that caused by 3D NLTE effects. The finer structures within the magnetic element are still visible after carrying out 3D NLTE computations, but are no longer visible after applying the DKIST PSF. Re-binning to the detector pixel size, which is nearly five times larger than the grid spacing in the MHD simulations, further degrades the image quality by reducing the smoothness of the intensity image, but has little effect on the contrast or, seemingly, on the resolution.

Fig. 11. Comparison of intensity maps at the rest wavelength of the 630.1 nm line from 1D NLTE and 3D NLTE computations. The top row shows images at the full resolution of the MHD cube, while for the images in the lower rows a spatial degradation has been applied, first by convolving with a PSF (middle row) and then by re-binning to the detector pixel resolution (bottom row).

4. Discussion and conclusions

In the present paper, we have investigated how strongly the horizontal transfer of radiation influences the visibility of structures at very small spatial scales in widely used photospheric spectral lines. To this end, we carried out 3D RT in a high-resolution 3D radiation MHD simulation box of a plage region with a horizontal grid spacing of 6 km, which is three times finer than in our previous studies. This allowed us to test also the effects of horizontal transfer on the finest spatial scales that are accessible with the largest current solar telescope, the DKIST on Maui.

First, we tested if the effects introduced by horizontal transfer and found in our paper III are reproduced at higher resolution. They are, and they are even clearer at higher resolution, which leads to sharper boundaries of magnetic elements, and thus to a cleaner separation of the gas properties inside and outside of the elements.

After that, we considered intensity contrasts and found that in the line cores the contrast generally decreases from 1D NLTE to 3D NLTE. On very small scales below 25 km, the decrease may be considerable, but on larger scales it remains small. Of the four lines considered here, the 630.1 nm line shows the largest reduction in contrast. This reduction is moderate, between 4–6%, at both full resolution and after spatial degradation.

Purely from the horizontal transfer effects in the photospheric layers, we do not see any evidence of the solar fog predicted by Judge et al. (2015). The solar fog is mainly caused by scattering, and in the photospheric lines analysed by us the effects of scattering are very small. The 525 nm line pair exhibits much less scattering than the 630 nm pair (see also Table 1 of Judge et al. 2015). For the 525 nm pair, UV overionization is the dominant NLTE effect and its source function closely follows the Planck function. The 630 nm pair is affected by both UV overionization and scattering. This explains why the 630.1 nm line, which is the stronger of the pair, shows a bigger reduction in contrast compared to all the other lines. However, the scattering coefficient (1 − ϵ), where ϵ is the collisional coupling parameter, is much smaller for the 630 nm line pair than for chromospheric lines. Even on very small scales at which we find a considerable reduction in contrast, it is still possible to see the large intensity variations produced by magnetic elements. Hence, the spatial smearing and reduction in contrast due to horizontal transfer effects resulting in a solar fog is of concern only for lines formed in the chromosphere.

Our investigations of the formation of photospheric iron lines in a high-resolution 3D MHD simulation conform with the findings from the 2D idealized FS models of Stenholm & Stenflo (1977), Bruls & von der Lühe (2001) and paper I, in the sense that magnetic elements display intensity variations on scales well below the photospheric photon mean free path. We conclude that the spatial resolution achievable in the studied lines is only slightly degraded by horizontal RT, even at scales well below the horizontal photon mean free path in a plane-parallel atmosphere. This degradation is clearly smaller than that introduced by the PSF of the telescope and the detector for the new generation of large solar telescopes (DKIST and EST). Therefore, the degradation due to horizontal RT should not hamper the capabilities of these telescopes to capture fine-scale structures smaller than the photon mean free path, at least in the photosphere.

Acknowledgments

The authors thank R. Cameron for running the high resolution MURaM simulations. HNS thanks L. P. Chitta and D. Przybylski for helpful discussions. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 101097844 – project WINSUN).

References

All Figures

	Fig. 1. Map of the vertical component of the magnetic field at the surface (log(τc500) = 0.0).
In the text

	Fig. 2. Continuum intensity at a wavelength close to the 525 nm lines. The squares indicate the two selected sub-domains that were calculated in full 3D NLTE. The bigger dashed squares indicate the domain used for calculations, while the smaller solid squares within the dashed squares represent the areas that were used for our investigation.
In the text

Fig. 3. Main atmospheric input parameters of the model sub-domain used in this work at heights z = 0 km (i.e. at the height of average τc = 1; panels on the right) and z = 240 km (the approximate line core formation region of the investigated lines; panels on the left). Plotted, from top to bottom, are the vertical component of the magnetic field, Bz, the hydrogen number density, nH, the vertical component of the velocity, vz, and the temperature, T.

In the text

	Fig. 4. Spatial distribution of the relative differences between the EW values of the Fe I 630.15 nm spectral line calculated in the two spatially shifted sub-domains in 3D NLTE. The colours denote the deviation in the EW between the two runs. The rectangle indicates the area chosen for our investigation, which lies completely inside the region with less than 1% difference.
In the text

Fig. 5. Enlarged section (≈1 Mm  ×  1 Mm) of the continuum intensity image (Figure 2) showing the part of the atmosphere where the 3D NLTE calculation was performed. The yellow line and the five yellow dots (a–e) indicate selected positions used in the text for further analysis. The yellow dots were drawn larger (2  ×  2) than in the actual image for better visibility. The black contours indicate the boundary of the magnetic element corresponding roughly to Bz(z = 0) = 1200 G.

In the text

Fig. 6. Selected quantities along a cut through a FS indicated in Figure 5. Top panel: Temperature (black) and density (blue). Middle panel: Vertical magnetic field (black) and vertical velocity (blue, negative values denote downflows). Bottom panel: δE1D − 3D (black) and continuum intensity, Ic (blue). All the quantities were taken at a geometrical height of z = 0 km. The vertical dotted lines indicate the boundary region of the flux element over which the density decreases and the magnetic field increases.

In the text

	Fig. 7. Intensity profiles of the 630.15 nm (upper row) and the 525.02 nm (lower row) lines at five selected spatial positions indicated in Figure 5. Ordering: (a) location with lowest y co-ordinate (granule) to (e) location with largest y co-ordinate (at the upper boundary of the FS). Legends explaining the colours of the curves are given in panel a.
In the text

	Fig. 8. Sample profiles of Stokes V/Ic of the 525.02 nm line at the two locations, (d) and (e). See Figures 5 and 7 for an explanation of the locations.
In the text

	Fig. 9. Comparison of RMS contrasts over a small region of the magnetic element for the three computation methods: LTE (first column), 1D NLTE (second column), and 3D NLTE (third column). In each row, we show the intensity images at the rest wavelength of the respective line, indicated in blue. The colour scale is the same as in Figure 10.
In the text

Fig. 10. Histogram of RMS contrast of small areas: we divided the full area of the intensity image into a grid of squares, each of which had a size of 5 pixels × 5 pixels. The grid is indicated as orange dots on the intensity image shown in the left panel. A histogram of the RMS contrasts computed at every square in the grid is shown in the right panel, for the three computation methods: LTE (black), 1D NLTE (red), and 3D NLTE (blue). The x axis has been clipped to show only bins with contrasts ≥10%.

In the text

	Fig. 11. Comparison of intensity maps at the rest wavelength of the 630.1 nm line from 1D NLTE and 3D NLTE computations. The top row shows images at the full resolution of the MHD cube, while for the images in the lower rows a spatial degradation has been applied, first by convolving with a PSF (middle row) and then by re-binning to the detector pixel resolution (bottom row).
In the text