© ESO, 2012

1. Introduction

The magnetic fields in both the solar network and active region plages is concentrated into elements that are generally described by magnetic flux tubes and sheets (Stenflo 1973; Spruit 1976; Solanki 1993). These magnetic elements form the footpoints of either coronal funnels and loops and chromospheric canopies. As such, they channel the energy transported from the solar interior to the upper solar atmosphere (Narain & Ulmschneider 1996) and play an important role in modulating the Sun’s total and spectral irradiance (Domingo et al. 2009).

Magnetic elements are thought to be formed by the magneto-convective processes of flux expulsion and convective collapse. Flux expulsion separates the magnetic flux from the convection (Parker 1963a,b; Weiss 1966; Galloway & Weiss 1981; Hurlburt & Toomre 1988), while convective collapse brings field strengths up to kG levels (e.g. Parker 1978; Spruit 1979; Grossmann-Doerth et al. 1998; Lites et al. 2008; Danilovic et al. 2010). The magnetic elements have a range of sizes down to as small as 100 km, close to the highest currently achievable resolution (Lagg et al. 2010), and probably lower than that too (Riethmüller et al., in prep.).

Radiation plays a dominant role in determining the properties of magnetic elements, both during the formation – i.e. during the convective collapse, where radiation leads to the cooling of the gas in the magnetic features, cf. Venkatakrishnan (1986) – and during the normal lifetime of a magnetic element. Since the kG magnetic field suppresses the convection within the magnetic elements, their photospheric layers are heated mainly by radiation. Owing to the evacuation produced by the excess magnetic pressure within the magnetic elements, radiation flows in through the walls and heats the interiors of the magnetic elements (Spruit 1976). Multidimensional radiative transfer (RT) plays a central role in setting this temperature and is a standard part of modern MHD simulations of magnetoconvection in the Sun’s surface layers (Nordlund 1982, 1983, 1985; Deinzer et al. 1984a,b; Vögler et al. 2004, 2005).

Diagnostic RT, however, is generally carried out in 1D or at most in so-called 1.5D, i.e. along many parallel rays passing through a multidimensional model (e.g. Bünte et al. 1993; Solanki 1989). This difference in the treatment of the energetics and the diagnostics is a basic inconsistency that is worth it to be investigated. We henceforth also use the term 1D for 1.5D RT. Although there have been a number of 3D RT studies of convective features (Carlsson 2009, and references therein), the multi-dimensional RT diagnostics of magnetic elements are to our knowledge restricted to only a few papers (Stenholm & Stenflo 1977, 1978; Bruls & von der Lühe 2001).

Owing to the very restricted computer capacities at the time, Stenholm and Stenflo in their pioneering papers had to make some compromises, restricting themselves to a strictly cylindrical flux tube and a two-level atom. In the present work, we use a more sophisticated (but still simple) flux-tube model where its expansion with height is computed consistently using the thin-tube (thin flux-sheet) approximation (Leighton 1963; Parker 1974a,b; Spruit & Roberts 1983, and references therein), a multi-level atom, and true 3D NLTE RT. A comparison with 3D radiation MHD simulations has shown that the thin-tube (or thin flux-sheet) approximation provides a reasonable description of the magnetic flux concentrations in the simulations (Yelles Chaouche et al. 2009).

The thin FS geometry was studied using 2D RT by Bruls & von der Lühe (2001) and part of our work resembles theirs. However, while Bruls & von der Lühe (2001) aimed mainly at establishing that features narrower than the horizontal photon mean free path can still be resolved, the aim of the present paper is broader. We aim to quantify the influence of multi-dimensional NLTE diagnostic RT relative to LTE and 1D NLTE RT. To this end, we carry out a parameter study probing the influence of magnetic field strength, flux-tube radius, flux-tube atmospheric model, and UV opacity on the Stokes profiles of four commonly used spectral lines. Although MHD simulations have reached a high level of realism, we prefer to be able to change parameters at will in a simpler model. This allows us to get a clearer understanding of the effects of multi-dimensional RT in comparison with 1D NLTE and 1D LTE. 3D diagnostic RT in hydrodynamic and MHD simulations will be treated in forthcoming papers.

The structure of our work is as follows: in Sect. 2, we provide a summary of our model ingredients. In Sect. 3, we describe the results of our parameter study in detail. In Sect. 4, we discuss the influence of the parameter variations on widely used quantities, such as contrasts and equivalent widths, focussing on the differences between three calculation methods (LTE, 1D NLTE, and 3D NLTE). Our conclusions are presented in Sect. 5.

2. Model ingredients

2.1. Radiative transfer calculations

All our calculations were done with the RT code RH adapted to our needs. This versatile code was developed by Uitenbroek (2000) and is based on the iteration scheme of Rybicki & Hummer (1991, 1992, 1994). It uses an improved version (Socas-Navarro et al. 2000) of the short characteristics (SC) method of Olson & Kunasz (1987) and Kunasz & Auer (1988) for Stokes I and Rees et al. (1989) for the full Stokes vector. The code is able to do LTE or NLTE calculations in 1D, 2D, and 3D. If not otherwise stated, the calculations presented here, were always done with the 3D version of the code (Uitenbroek 2006). In addition to several adaptions for convenience or to enhance the performance, we made a principal change to the formal solver that has a considerable influence on our results. In the Appendix, we summarize our improved method for the formal solution used for both Stokes I and full Stokes vector calculations, as well as some pitfalls we encountered with both the old and new methods.

2.2. Atomic model and iron abundance

We used an Fe i atom including 22 Fe i levels and the ground state of Fe ii, which allowed us to perform a realistic NLTE calculation of selected multiplets containing transitions of interest. The model, kindly provided by J. H. M. J. Bruls, employs energy level data compiled by Corliss & Sugar (1982) and line data from Fuhr et al. (1988). It was enhanced with oscillator strength data from Thevenin (1989, 1990).

We selected the well-studied and often observed 524.7 nm, 525.0 nm, 630.1 nm, and 630.2 nm lines for this work, which are often analyzed either individually or as line pairs (Stenflo 1973; Elmore et al. 1992; Tsuneta et al. 2008; Solanki et al. 2010; Martínez Pillet et al. 2011). Altogether, the atom contained 23 continuum and 33 line transitions with a grand total of some 1300 wavelength points (we tried to keep this number small to reduce the computing resources needed). In this model atom, the oscillator strengths were as displayed in Table 1. The f-values of the 524.7 nm and 525.0 nm lines were enhanced by roughly a factor of from two to three against those found in the literature (Blackwell et al. 1979), following the correction procedure used in Bruls & von der Lühe (2001). An iron abundance of 7.44 was employed following Asplund et al. (2000). We note that we were mainly interested in differences between the line profiles obtained in LTE, 1D NLTE, and 3D NLTE under different conditions and not in the absolute line profiles. In particular, we did not compare our line profiles with observations, so that details such as the exact elemental abundance or f-values are not of central consequence and do not influence our conclusions.

The well-known problem of UV overionization (Bruls et al. 1992), which leads to wrong population numbers because of missing iron line opacity in the UV, was eliminated by enhancing the opacities in the relevant wavelength range (opacity fudging) as proposed by Bruls et al. (1992). The resulting differences were calculated and are discussed in Sect. 3.2.

2.3. Model atmospheres and computational setup

We studied both magnetic flux tubes (displaying axial symmetry) and flux sheets (displaying translational symmetry). In both cases, we employed the thin-tube approximation to describe their internal structure. For most computations, however, we used magnetic fields only to determine the Wilson depression and the expanding geometry of the flux tube, to better isolate the influence of multidimensional RT, as in Stenholm & Stenflo (1977). In Sect. 3.9, however, we also consider results computed by including the formal solution of the Unno-Rachkowski equation, i.e. where we computed Stokes profiles using the so-called field-free method (Rees 1969). Thus, we included the influence of the Zeeman effect on the line profiles (splitting, polarization), although the results presented here neglect the influence of the Zeeman effect on the population numbers.

The above-mentioned flux sheet could in principle be calculated with a 2D code as long as one considers radiation emerging only along rays lying in a plane perpendicular to the plane of symmetry. However, to achieve an unbiased comparison with flux tubes, we used the 3D code version. The qualitative effects in a 2D flux sheet (FS) are the same as in a true 3D flux tube (3D FT), but are achieved at much lower computing costs. We note that an equidistant grid (in the xy-direction enforced by the code) is not very well-suited to mapping an expanding tube or sheet, therefore the number of grid points has to be quite high. To achieve reliable results, we used models with 288  ×  4  ×  256 grid points in the flux sheet geometry, where the z-axis denotes the height in the atmosphere (see also Fig. 1 for a cut across the flux sheet). In Sect. 3.7, we present the results for a true 3D flux tube geometry having 144  ×  144  ×  128 grid points in the cartesian geometry of the RH code. Although only half the number of pixels are introduced along the x and z axes when compared to the flux sheet case, the total number of grid points is an order of magnitude larger. To keep the computing time from growing accordingly, the number of rays was reduced (standard set A6 of RH code instead of 54 vertical directions per octant, see Sect. 2.3.2).

To compare our 3D with our 1D calculations, we constructed plane-parallel atmospheres by vertically cutting the 3D cube at selected x-positions. The 1D NLTE RH code was then run with the so-constructed 1D atmospheres as input (we stick to the term 1D NLTE for this method, which is often called 1.5D NLTE). In addition, for comparison purposes 1.5D computations in LTE were carried out for each flux sheet or tube (henceforth simply called LTE computations).

2.3.1. Standard model

We varied a range of parameters such as radius, Wilson depression, thickness of the walls, filling factor, and flux sheet atmosphere (see Sect. 3 for details). As a reference, we defined a standard model with a fixed parameter set:

Atmosphere: if not otherwise stated, the inner and outer atmospheres were analogous to those in Stenholm & Stenflo (1977) and therefore identical. This choice allowed us to more clearly segregate the effects of horizontal RT. We used the Harvard-Smithonian reference atmosphere (HSRA) described by Gingerich et al. (1971). In the field-free part, the model was extended downwards by combining it with the convection zone model of Spruit (1974). Figure 1 shows the geometry and the temperature distribution in the xz plane for our standard model. The thick solid line follows τ_c = 1 and outlines the smooth but quite steep transition between the flux sheet interior and exterior. Inside the nominal boundary (indicated by the vertical lines as in many subsequent figures), the hot walls of the outer atmosphere, positioned below the continuum level of the quiet Sun, lie close to the local τ_c = 1 level, so that their presence strongly influences the radiation field sensed by an atom within the magnetic feature. A constant microturbulent velocity of 0.8 km s^-1 was added to every grid point of the model in direct correspondence with the original model. In the following, for simplicity, we refer to the outer, field-free atmosphere as either the quiet Sun (QS) or QS part. We note that the outer atmosphere ceases to exist at the height where the canopy fills the whole width of the atmosphere except for one single pixel at the edge of the atmosphere that was chosen to be pure HSRA without any magnetic field.

Field strength: the magnetic field strength of our standard sheet was B = 1600 G at ${\mathit{\tau }}_{\mathrm{c}}^{\mathrm{QS}}\mathrm{=}\mathrm{1}$ and B = 0 G outside, resulting in a Wilson depression of 191 km. This B value cannot be directly compared with measured intrinsic field strengths, which depend on the height of formation of the employed spectral line(s), diagnostic method, etc. The field strength dropped with height and consequently the magnetic element expanded with height. We note that a flux sheet exhibits a more rapid expansion with height than a round flux tube as a consequence of flux conservation if one assumes a fixed run of field strength with height because B × R = const. for a sheet in contrast to B × R² = const. for a tube, where R is the radius of the FT or the horizontal location of the FS boundary.

Wall thickness: the sheet boundaries were smoothed in the x-direction with a Fermi function which drops from 90 to 10 percent within 20 km, that is slightly above the minimal value of the wall thickness estimated by Schüssler (1986).

Radius, filling factor, and grid points: our standard model had a radius of R = 100 km and a filling factor α = 0.2. The grid points were equidistant in each dimension, where dx = 3.47 km, dy = 100 km (irrelevant for a 2D geometry), and dz = 5 km, starting at 460 km below τ_c = 1 of the field-free part of the model. Such a low starting point was needed to ensure that the model is optically thick at the bottom, even in the presence of a strongly evacuated flux sheet.

2.3.2. Angular discretization

Multi-dimensional problems basically need higher angular resolution to achieve the same accuracy of angular quadrature as in a 1D problem (van Noort et al. 2002). The very specific geometry of FS/FTs compelled us to increase the number of rays used in the simulation. We used a grid of 27 vertical (i.e. inclination) and 2 horizontal (i.e. azimuthal) angles per octant in the FS models, i.e. 216 outgoing rays from each grid point. This small number of horizontal angles was adequate in the (2D) FS geometry. Figure 2 illustrates the difference in outgoing intensity between our chosen angular resolution and one obtained with the often used A6 standard set of the RH code with only 24 outgoing rays (cf. Bruls & von der Lühe 2001). Close to the flux sheet boundary, the difference is significant (≈7%). This difference has its origin in the very hot walls of the FS in this model and the optical thinness of the sheet’s interior. The part where the walls are really hot is very localized, being restricted to a few grid points. Because the strongly evacuated sheet interior is optically thin, the hot walls can be seen from far away, as long as the ray travels only upwards through the sheet interior. The discretization of the rays gives rise to different lighthouse-like rays through the sheet. To see the hot walls at every point within the sheet it was necessary to increase the number of vertical rays. Some locations in the interior would not otherwise have been hit by a ray originating in the hot wall, which, in turn, may affect the local atomic level populations and the corresponding line depths. We note, however, that the artificially suppressed NLTE effects at these points are only relevant when the optical thickness of the interior is not negligible.

For the 1D NLTE calculation, we used a standard set of five Gaussian angles for all calculations.

3. Results

3.1. General features

We like to compare our results for the LTE, 1D NLTE, and 3D NLTE RT computations. We first concentrate on the Fe I 525.0 nm line formed in our standard atmosphere (see Sect. 2.3.1), but later also consider the other spectral lines and models. Figure 3 shows the intensity along the x-axis of the standard flux sheet (i.e. cutting across the sheet) for three selected wavelength positions, namely in the line core, in the line flank at approximately half depth (~0.004 nm from the line center), and near the continuum (~0.008 nm from the line center). The nominal flux sheet boundary at τ_c = 1 is indicated by vertical lines. The center of the flux sheet is located at x = 0. Only half of the cut through the sheet is plotted since the other half is simply a mirror image owing to the symmetry of the model.

In the line core and flank, both 1D and 3D NLTE computations show similar, strong variations along x (Fig. 3), as reported by Bruls & von der Lühe (2001), with the rest intensity reaching a peak value of two to four times that found in the QS far away from the FS. We note that quantitative differences between the results of Bruls & von der Lühe (2001) and ourselves may be caused by the quite different atomic and atmospheric parametrizations, the most important probably being the Wilson depression (magnetic field strength). In the line wing, the variation is slightly higher than in the line core. The LTE intensity, however, shows only a minor variation along the x-axis in both the line wing and core. Far away from the FS, all three methods of computation give almost identical results, only the LTE rest intensity remaining slightly below that of the NLTE methods. In the continuum, all three methods are identical and show a sharp rise in the rest intensity close to the boundary of the FS (cf. Caccin & Severino 1979).

Finally, a large difference between the 1D and 3D NLTE calculations is found inside the FS, with the 3D profile being weakened (i.e. having a higher rest intensity than the profiles computed with the other methods).

The upper panels of Fig. 4 show the corresponding spectra for five selected locations indicated above the panels. In all plotted cases, the LTE line profile saturates, maintaining nearly the same rest intensity (as seen in Fig. 3a, while the NLTE profiles become much weaker at some locations. We note that the relative differences in the line depth (see Fig. 4) are considerably smaller than those in rest intensity (Fig. 3) owing to the generally low line core intensity.

We now attempt to understand the intensity variations plotted as a function of x in Fig. 3. Far away from the flux sheet, the line is formed in the QS part of the atmosphere just below the optically thin canopy, which lies at a height of 400 km. The difference between the LTE and NLTE profiles at that location is basically the same as one gets when using the HSRA on its own, without any FS. As we move closer to the sheet, the height of the canopy’s base decreases continually. Unit optical depth in the line is reached at lower heights (below the canopy and, for a sufficiently low canopy base, at its boundary) owing to the lower density and hence opacity of the canopy, which remains optically thin to the line right up to positions slightly inside the nominal FS boundary. The line core is therefore formed in the hotter layers of the external atmosphere and consequently has a higher line core intensity, which culminates just inside the sheet’s nominal boundary.

The 3D NLTE contribution functions in Fig. 5a illustrate the decreasing height of the line-core formation region as we move towards the sheet. We note the peculiar shape of the contribution functions far away from the sheet. In contrast to the typical case in which the contribution function first sharply increases when moving up in the atmosphere, then slowly decreases above the height of maximum contribution, these contribution functions show a steep edge at the higher end, which corresponds to the lower boundary of the canopy. Examples are found between the two highlighted contribution functions (solid lines) at x = 465 and x = 150 km.

The position of this edge moves down as we move closer to the sheet in complete accordance with the geometry of our atmosphere (see also Fig. 1). The intensity maximum of I_core for the NLTE-cases lies slightly inside the flux sheet, where we – looking down vertically – directly see the hot walls below the τ_c = 1 level of the QS part. All three contribution functions (at the line core, at half-depth, and in the continuum, i.e. Figs. 5a−c, respectively). At that position, the contribution functions are peaked almost exactly at the same height, confirming the hypothesis that we see light formed below the height at which the continuum is formed far away from the FS. We note that the spectral line does not completely disappear owing to a small difference between the line and continuum formation levels and the very strong temperature gradient at the hot walls.

Why do we not see the same effect in the line formed under LTE? As the population numbers of Fe i are generally higher in LTE (owing to missing ionization by UV radiation), the line is slightly stronger than in the NLTE case. This leads to the situation where the canopy starts to become optically thick in the line core even far from the FS. Closer to the FS, the optical thickness of the canopy increases as the height of its base decreases. Because the temperature in either the canopy or the FS at heights of around 300−400 km is only slightly cooler than in the external atmosphere at the same geometrical height, the line core intensity decreases only slowly when we move towards the sheet (i.e. the line becomes slightly stronger).

Moving away from the line core, the opacity decreases, so that in the LTE continuum the intensity shoots up, as in the NLTE line cores (Fig. 3c). The three identical continuum curves in Fig. 3c indicate that the continuum is as expected formed in LTE. In the line flanks, the LTE opacity takes on an intermediate value so that the intensity still increases when horizontally approaching the FS boundary until reaching a gentle maximum at the FS boundary (see Fig. 3b).

The main effect of 3D RT is seen within the sheet. There the LTE and 1D NLTE profiles attain nearly the same depth as they have far from the flux sheet, as expected for a model with the same atmospheres inside and outside the flux sheet. Inside the flux sheet the LTE and 1D NLTE profiles are independent of x because the local atmosphere does not change (when neglecting the wall smoothing). The 3D NLTE line profile, however, remains significantly weaker (cf. Fig. 4e) owing to the ionization of neutral iron by the UV radiation flowing in sideways from the hot walls of the flux sheet, as discussed by Stenholm & Stenflo (1977). The line-core contribution function becomes very broad, extending from the continuum forming level at −200 km up to  + 300 km in the sheet.

In a small region at the FS boundary, the 1D NLTE intensity exceeds the 3D case. In the hot walls, both the horizontal RT towards the cool sheet interior and radiation from the cool interior of the flux sheet (which are both not present in the 1D case) decrease the amount of UV radiation reaching an atom there, thus giving rise to the opposite effect with respect to that seen inside the flux sheet.

We conclude this section by restating Stenholm & Stenflo (1977) since they provided a succinct summary of the results presented so far: “even when a line in a plane-parallel atmosphere is formed in LTE (i.e. LTE is a reasonable approximation), there is no guarantee that LTE is still a good choice in more complex geometrical configurations”.

3.2. The role of UV opacity

It is important to judge the influence of enhancing the UV opacity following Bruls et al. (1992) to compare with the work of Stenholm & Stenflo (1977), who used a simple two-level atom without any enhanced UV opacity. The role of UV opacity is to partially suppress the influence of the UV radiation from the deeper, hotter layers and therefore reduce the NLTE effect of UV overionization.

The spectra in the lower panels of Fig. 4 are calculated without any enhanced UV opacity. Even far away from the sheet, the NLTE-profiles differ rather significantly from the LTE profiles when no opacity enhancement is introduced. They start to differ dramatically when approaching the sheet (at 150 km). The difference between the 1D and the 3D calculations in the sheet is also strongly enhanced due to the strong influence of the horizontal transport of UV radiation in the 3D case. The depth of the line in the 3D case is smaller by more than a factor of two. These results are more comparable to those obtained by Stenholm & Stenflo (1977). Obviously the employment of a more adequate UV opacity model reduces not just the strength of NLTE effects in general, but in particular the effect of the horizontal RT, which remains important nontheless inside the magnetic element.

We note that the influence of the UV opacity enhancement on the Fe i 630.1 nm line, which is principally more sensitive to NLTE effects, in particular horizontal transport, is slightly smaller than that on the 525.0 nm line. The 630.1 nm line is already in a partially saturated regime and therefore removing the UV-opacity enhances the NLTE effects to a lesser extent.

3.3. Behaviors of different lines

The discussion in Sect. 3.2 already indicates that the effects of NLTE and 3D RT depend strongly on the formation height, so that different lines may exhibit different behaviors. In this section, we compare the profiles of three widely employed lines (Fe i 525.0 nm, Fe i 630.1 nm, and Fe i 630.2 nm). The results for the 524.7 nm line turn out to be exactly the same as for the 525.0 nm line, which is unsurprising, given that the two lines are nearly identical, save their Landé factors (which do not play a role here, since we set B = 0 G for the line-formation computations). The 630.1 nm line is expected to be somewhat stronger, and the 630.2 nm line to be somewhat weaker than the 525.0 nm line. In addition, the two 630 nm lines have a considerably higher excitation potential, which may also influence our results.

Figure 6 shows the line core intensities along the x-direction as in Fig. 3a, but now for the three lines 525.0 nm, 630.1 nm, and 630.2 nm. We repeat the results for the 525.0 nm line from Fig. 3a for ease of comparison. The global features are the same for all three lines. Far away from the sheet, the lines are formed in a very similar way (although the 630 nm doublet is somewhat stronger in NLTE relative to its strength in LTE). Close to the flux sheet boundary, there are, however, notable differences.

A major difference of the 630.2 nm line in comparison to the other lines is the generally higher intensity level owing to its weaker opacity. The LTE curve confirms this view, as can be deduced from the increasing I_core level as we move towards the sheet: this is because the LTE line is formed (partially) in the hotter gas below the canopy and not mainly either in the canopy itself or its boundary where the cores of the other lines are formed (see also the discussion in Sect. 3.2). The significantly stronger 3D effect seen in this line is quite striking inside the FS compared to that for the 525.0 nm line.

The 3D-NLTE case of the 630.1 nm line differs quite strongly from those of the other lines. Although the intensity increases towards the sheet, there is no peak near the flux sheet boundary, which remains at almost the same level over the whole sheet. The line-core contribution functions drawn for the three lines in Fig. 7 indicate that the weaker lines have strong contributions from the continuum near the boundary (at x = 94 km), in contrast to the stronger 630.1 nm line. The contribution function in the line core is clearly shifted towards higher levels, i.e. we do not see directly the hottest parts of the walls. We note that outside the sheet (e.g. at x = 150 km), the 630.1 nm line is formed only marginally above the other lines. Inside the sheet, the effect of horizontal RT weakens this line, leading to higher intensities in the 3D-NLTE case.

Figure 6b indicates that at certain x-positions the LTE intensity of the 630.1 nm line is slightly higher than in NLTE. This may be explained by mainly two different effects (Jo Bruls, priv. comm.). The photospheric source functions for these lines usually follow the Planck curve for both LTE and NLTE. In the NLTE case, owing to the lower population of Fe I, the τ = 1 level is shifted downwards, into hotter regions, weakening the line. On the other hand – especially for the stronger 630.1 nm line formed higher up – the NLTE source function may clearly fall below the Planck function, leading to the opposite effect. Which effect wins in the end is determined by the complex interplay of many transitions and levels in the atom and is difficult to forecast.

3.4. Influence of the wall thickness

The true thickness of the walls of magnetic elements is not well-known. The minimum thickness based on the expected solar magnetic resistivity, diffusivity etc. was estimated by Schüssler (1986) to be on the order of 3−10 km. In addition, radiation is expected to smooth the transition from the surroundings to the FT/FS and convection is expected to make the wall irregular, although the profile shape of the g = 3 line at 1.56   μm suggests that the magnetic field drops rather suddenly at the edge of the magnetic element (Zayer et al. 1989). In our calculations, we used smoothings with three different Fermi functions, each defined by the distance over which one atmospheric component changes from the 90% to the 10% level. In total, thicknesses of 40 km, 20 km (standard model), and 10 km were tested. We expected the wall thickness to have a strong influence on the lines where the continuum radiation of the hot wall could be seen directly, as in the case of the 525.0 nm or the 630.2 nm lines. By changing the distance over which the atmosphere changes from the cold interior to the hot QS outside the sheet, we directly changed the temperature at the locations where the continuum radiation forms. We also expected that the effect of horizontal RT would decrease when we expanded the thickness of the flux sheet boundary, since a smoother transition brings the FS ever closer to the plane-parallel, i.e. the 1D case. The resulting intensity is determined by a complex interplay of many factors (see also Bruls & von der Lühe 2001).

Figure 8 shows the LTE, 3D NLTE, and the 1D NLTE line-core intensities of the 525.0 nm line along the x-axis for the three different wall thicknesses. As expected, the strongest influence is seen at the hot wall itself. While the LTE curves almost coincide, the 1D and 3D NLTE curves show a systematic increase in the intensity maximum with decreasing wall thickness. Interestingly, the change of the 3D intensities in the inner part of the sheet is, however, much smaller. The largest effect is displayed by the 524.7/525.0 nm lines. The other lines (not shown), especially the 630.1 nm line, are less sensitive to the wall thickness.

3.5. Variation in the radius

We next considered FS models with different widths. To change the width, the flux sheets were stretched in the x direction. Therefore, the thickness of the walls scaled with the width of the FS. All other quantities remained the same.

Core intensities at three locations along the x-axis are plotted in Fig. 9 versus the FS half-width of three spectral lines. The selected locations are the sheet center (plus signs), the position of the intensity maximum near the wall of the FS (diamonds), and the QS far away from the sheet (asterisks). The latter does not show much variation with radius as a consequence of its position far away from the sheet. For r ≥ 400 km, the intensity at the FS center almost corresponds to the 1D NLTE case. As expected, these broad flux sheets are optically thick in the horizontal directions, so that the hot walls are no longer visible at the center of the FS, hence do not influence the radiation field there.

The intensity maximum near the boundary (diamonds) displays a similar behavior. This can be understood by bearing in mind that the wall thickness increases with FS width, so that for wider tubes the walls are not so hot (see Sect. 3.4). However, horizontal transport plays a role as well, since an observer located close to one wall also sees radiation from the opposite wall. With increasing FS half-width, the strength of this irradiation decreases. For the intensity maximum at the flux sheet boundary, the effect of horizontal RT seems to saturate at R ≈ 200 km. This is consistent with the behavior seen at the center of the flux sheet where the effect saturates at double the distance, i.e. at R ≈ 400 km. Therefore, the opposite wall (at 2R) is invisible at the boundary and can no longer influence the line for broader FS.

A few test computations with FS of different half-width, but with the same wall thickness also show a decrease in the core intensity with FT size. The decrease is only slightly less steep than in Fig. 9, confirming that the main effect is due to the FS size and not due to the wall thickness, as already suggested by a comparison of Figs. 9 with 8.

3.6. Variation in the temperature stratification

We then modified the model by replacing the atmosphere in the FS interior by, in turns, the network model of Solanki (1986) (from now on referred to as NET) and a plage model according to Solanki & Brigljevic (1992; hereafter PLA). The nominal field strength of the FS (B = 1600 G) and the external QS atmosphere (HSRA) were unchanged. Although these semi-empirical models are relatively old, they agree well with more recent 3D radiation MHD simulations (Vitas et al. 2009). The NET model is the hottest among the three selected atmospheres, being hotter than both the HSRA and the PLA model at all τ. The PLA model is similar to the HSRA at heights around τ_c = 1, but has a flatter temperature gradient and is therefore hotter than the HSRA at the formation height of the 525.0 nm line core, which was taken into consideration here (see Fig. 10).

Far away from the FS, the rest intensities are similar for all models and computational methods, which is unsurprising as the canopy lies so high that the FS atmosphere is almost irrelevant. Close to or inside the FS, however, the differences are significant for both the different models and methods.

We first considered the LTE case, which is easy to understand because the different temperatures of the inner atmospheres are directly visible in the resulting rest intensities. The closer to the FS we are, the more we see the influence of the FS atmospheres. In accordance with the higher temperature in the PLA and NET models, the rest intensities increase towards the FS (right panel in Fig. 10). The highest intensity is reached near the FS boundary and for the hottest model (NET). The weakening of Fe i lines in hotter atmospheres is due to the enhanced ionization of iron.

In the NLTE case, the rest intensities vary strongly as well, but only in the interior of the FS. For the NET model, the NLTE intensity lies even below the LTE value. This can be explained by the large temperature difference between the NET model and the HSRA at the geometric height of line core formation. The downward shifted NET model in the FS is still approximately 200 K hotter than the external HSRA model. We note that this no longer applies at greater depth. In the sheet, the LTE line-core intensity is therefore even higher than the corresponding 3D NLTE intensity because the flux sheet interior is “cooled” by horizontal RT towards the now “cold” walls. The usual line weakening by the horizontal transport of UV radiation, as seen, e.g., in the HSRA-HSRA reference case, is inverted and the net horizontal transport of radiation has the opposite sign.

The surprisingly good agreement between the PLA-LTE and the PLA-3D-NLTE rest intensity leads to the conclusion that the PLA model – when calculated in LTE – compensates the most effectively for 3D and NLTE effects. The other lines considered in this paper also agree quite well (not plotted). This suggests that the PLA model of Solanki & Brigljevic (1992) also remains valid as a description of the atmosphere of magnetic elements in full 3D NLTE, at least for the spectral lines considered here.

3.7. Influence of geometry: flux tubes

We present results for a flux tube (FT) geometry, i.e. for magnetic elements with a round cross-section. Despite the lower resolution of the atmospheric model (cf. Sect. 2.3.1), the tube may be represented in a fairly smooth way because the filling factor of 0.2 leads to a considerably larger FT fraction when we look along only one dimension (e.g. x-axis). The hot walls are therefore resolved almost as well as in the FS case. All other model ingredients (atomic model, atmospheric components, radius, magnetic field, filling factor, and smoothing) remained unchanged from our standard FS model.

Figure 11 shows maps of the core intensities of three lines computed in 3D NLTE (top) and LTE (bottom). The grayscales for the 3D NLTE and the corresponding LTE images are the same for a given spectral line to facilitate the comparison. Owing to the rather different maximum intensities reached by the different lines, the greyscales were normalized individually. The maximum intensity per line – typically reached near the FT boundary in the 3D NLTE case – was used as the reference value (100%). Figure 12 displays the same intensities along a cut passing through the axis of the FT. We observed the same effects as for the standard flux sheet, although the line core intensities in the tube itself are higher than for the FS, which may be explained by the more intense irradiation from the now completely surrounding hot walls. Differences near the FT boundary originated from the different geometrical expansion of the FT in three dimensions.

This confirms that the results obtained for the (2D) FS geometry apply qualitatively to the 3D FT case as well. The horizontal transfer effects, i.e. those of 3D RT, are actually larger in FTs than for a FS (see e.g. the line weakening of the 525.0 nm line at the centers of the FT and FS).

3.8. Center-to-limb variation

Figure 13 depicts the synthetic images of a flux tube in the intensity (3D-NLTE calculation) of the 525.0 nm line core for three different viewing angles. In the following, γ = arccosμ denotes the angle between the vertical (z-axis) and the line of sight (LOS) and χ the azimuthal angle in the xy plane measured from the x-axis. In the leftmost panel, showing the vertical view (i.e. μ = 1.0), we see the strong and very localized (≈10−20 km) intensity maximum near the FT boundary and, less strikingly, but still clearly visible, the enhanced intensity in the interior of the FT. For μ = 0.93 (i.e. γ = 21.5°) and χ = 45° (middle panel), we still see half of the ring, which is now much broader (≈50 km) and brighter. It actually resembles an arched structure such as a small loop. The rest of the hot wall is hidden from view by the optically thick cooler gas of the external atmosphere in the foreground. For μ = 0.7 (γ = 45°) and χ = 0° (right panel), the hottest parts of the wall cannot be seen directly anymore, leading to a more diffuse intensity distribution, partially produced by the warm walls at greater height, partially due to the horizontal RT. We note that for all LTE calculations, no significant structures could be recognized at all. The same is true for NLTE computations at μ < 0.7 (γ > 45°).

Figure 14 shows the intensity distributions for the same (and more) viewing angles, but now for the 2D flux sheet. In the FS geometry, there may be different intensity distributions for a single μ and different azimuthal angles, as illustrated for μ = 0.68. When the ray is mainly in the x-direction, i.e. χ is small (dotted curve), we see the same smooth distribution as for the 3D FT (rightmost panel of Fig. 13). If we look at the same μ = 0.68, but now mainly along the y-direction, i.e. χ is large (dashed curve), we have a more direct view of one of the hot walls because the ray is oriented along the central plane of the flux sheet. Figure 14 shows an additional viewing angle with a very small μ of 0.26 (γ = 75°, thin solid line), where one barely recognizes any variation, since the hot wall is almost completely hidden by the front part of the FS.

3.9. Solution for the full Stokes vector

Here we discuss solutions for the full Stokes vector including Zeeman splitting (following the method described in Sect. 2.3 and the Appendix in more detail).

Figure 15 displays the run of the line-core intensities of our four selected lines for different calculation methods. Each panel pair compares the results obtained by including Zeeman splitting (panels on the right) with those without it (panels on the left). We note that the non-magnetic case (Stokes I only) was already discussed in Sect. 3.3. We now show both the 524.7 nm and 525.0 nm lines as their profiles are no longer identical in the Zeeman split case due to their different Landé factors.

The inclusion of Zeeman splitting may lead to either an intensification or attenuation of the NLTE and horizontal RT effects, but the concrete behavior of a line can be understood in terms of its strength, Landé factor, and to a certain extent its temperature sensitivity.

We first consider the 524.7 nm and 525.0 nm lines. In the absence of Zeeman splitting these two lines are practically identical, as expected (panels a and c). In the panels b and d, however, the 525.0 nm line is more strongly weakened owing to its higher Landé factor. This is true for all three curves. The Zeeman splitting essentially reduces the opacity in the line core. This may allow the hot walls to become visible to the line core, as happens in LTE for the 525.0 nm and 630.2 nm lines, where we observe strong peaks near the boundary in the Zeeman-split case.

At the FS center, the effects of 3D RT (i.e. the relative difference between the 1D and 3D intensities) may also be either enhanced or damped. The 630.1 nm line, which is the most Zeeman insensitive of the four considered spectral lines, displays little dependancy on Zeeman splitting in LTE and 1D NLTE, but a strengthening of the 3D effects. In contrast, precisely the opposite behavior is exhibited by the 630.2 nm line, where the relative difference between the three methods is very small when Zeeman splitting is included. The 524.7 nm and 525.0 nm lines are particularly affected by 3D effects. We note that these curves refer only to the line core. At other places in the line profile, even the 630.1 nm line appears to be considerably influenced by Zeeman splitting. The exact behavior depends on details of the splitting compared to the line width, as well as other factors.

Figure 16 demonstrates the variation in the Fe i 525.0 nm Stokes I profiles with distance from the x-axis of the FS, in analogy to Fig. 4, but for the solution of the full Stokes vector. As expected, far away from the FS, the profiles computed in 3D are very similar to LTE profiles (panel e). The commonly made assumption that Fe i 525.0 nm is “formed close to LTE” seems to be valid. The inclusion of the Zeeman effect has little influence because the line is formed mainly in the non-magnetic atmosphere below the canopy. The situation changes when we approach the boundary of the FS where the field strength is higher and the line is partially formed in the canopy. In addition to the effects already discussed in Sect. 3.3, the inclusion of the full Stokes vector in the formal solution produces much broader profiles and enhanced line-core intensities (especially in the 3D NLTE and 1D NLTE calculation, less so in LTE) at all locations close to or inside the FS (compare panels a) to c) of Fig. 16 with the corresponding panels of Fig. 4).

Figure 17 displays the Stokes V/I_c and Fig. 18 the Q/I_c profiles corresponding to the first four I profiles in Fig. 16. Positions further away from the FS do not show a significant polarization signals, as the line is formed in the field-free part of the atmosphere.

The 3D NLTE V/I_c profiles are weaker than the LTE ones at all locations and, inside the FS, also weaker than the 1D NLTE ones (except for the special case very close to the boundary, where the influence of horizontal RT is reversed). In the 525.0 nm line, the larger Zeeman splitting in the 3D profiles can be clearly seen (both in Stokes I and V). In general, however, owing to 3D RT, Stokes V displays smaller effects than Stokes I because these effects are largest in the line core, where Stokes V is zero. Interestingly, the form of the 525.0 nm V profiles in panel e) to differ from those calculated in either 1D or LTE or those of the 524.7 nm line. The center of gravity of the 525.0 nm 3D profile shows a systematic shift away from the line center, an effect that may be of considerable relevance to the diagnostics of the magnetic field strength, since such a profile mimics a stronger field strength and a seemingly lower filling factor than the profile in LTE. The Stokes V 3D NLTE profile shape in the FS center is consistent with an enhanced weakening of the line core inside the FS.

The Stokes Q/I_c profiles plotted in Fig. 18, clearly illustrate that the effect of 3D RT is largest near the line core. At the center of the FS, these profiles vanish for reasons of symmetry (the plot does not exactly corresponds to the FS center, so that the signals do not completely disappear). Near the flux sheet’s boundary, however, as shown in panel b, the strong weakening of the 3D NLTE Stokes Q in the line core is clearly apparent. In the line wings, unsurprisingly, the 1D and 3D profiles almost coincide. The difference in the strength of the π-components of the 524.6/525.0 nm lines is due to the Zeeman desaturation of the line core, which is stronger for the higher g line.

3.10. Variation in the field strength

The field strength has a two-fold effect. Firstly, it affects the construction of the atmospheric models by determining the Wilson depression, which in turn influences the magnitude of the horizontal transfer effects. The Wilson depression corresponding to the different field strengths is listed in Table 2. Obviously, the dependence is strongly non-linear. Secondly, the field strength directly affects the line profile through the Zeeman effect. Here we considered the effect of changing the field strength, while other atmospheric parameters were unchanged.

Figure 19a illustrates the influence of the magnetic field strength on the 525.0 nm line core intensity without Zeeman splitting. As expected, the influence of the field strength is strong even when neglecting the Zeeman splitting. The main effects – the strong peak near the flux sheet’s boundary and the enhanced intensity in the center of the FS – both increase with increasing Wilson depression. Below approximately 1500 G at z = 0 km (which corresponds to τ_c = 1 in the HSRA), the peak due to the hot walls disappears, while the intensity produced by 3D RT at the FS center remains somewhat higher than in 1D NLTE even for 1000 G, demonstrating that horizontal radiative transport is still a relevant factor.

Figure 19b depicts the intensities obtained when the line splitting due to the Zeeman effect is taken into account. The full Stokes solution provides significantly higher rest intensities close to and inside the FS compared to the solution neglecting the Zeeman effect, whereas the intensities far away are unaffected as they are mainly formed in the – non-magnetic – QS atmosphere. We note that only part of the increase in rest intensity in Fig. 19b is directly due to Zeeman splitting. Its indirect influence in terms of 3D RT is equally important, as can be seen by comparing Figs. 19a with b.

4. Discussion

After presenting, in Sect. 3, the influence of various model parameters on both the line profiles and the rest intensity, we now consider questions such as, how strongly do typical observational quantities differ if one employs true 3D RT instead of 1D or even LTE? As Stenholm & Stenflo (1977) already pointed out, not only the inclusion of NLTE effects, but also horizontal RT may play major roles when considering a flux tube geometry.

Here we have gone beyond the earlier studies of Stenholm & Stenflo (1977), Caccin & Severino (1979), and Bruls & von der Lühe (2001) by carrying out a systematic study of the changing parameters of a thin flux tube (or thin flux sheet) on the profiles of four widely used spectral lines at different levels of simplification in the RT (full 3D NLTE, 1D NLTE, and LTE).

To compare the influence of different calculation methods and models on the determination of important atmospheric parameters, such as temperature, abundances, magnetic flux, and field strength, we selected frequently used observables to quantify the effect of the different parameters. We show the results for the Fe i 525.0 nm line, which is a line widely used in solar diagnostics (e.g., Barthol et al. 2011). The effects are similar for the other lines, though they may differ in their details. Zeeman splitting and UV opacity enhancement were included in all computations presented in this section.

4.1. Core line depth and determination of temperature

Table 3 gives an overview over the line depth d_l = (I_c − r)/I_c, where r is the rest intensity at line center and I_c the continuum intensity, calculated for some of our model FSs. We consider here the depth of the spatially averaged profile, i.e. averaged over the area over which the flux sheet expands until it merges with its neighboring sheets. For a filling factor of 0.2, which is assumed for all computations presented in this section, this corresponds to integrating out to five times the distance from the FS boundary (at a height corresponding to τ = 1 in the QS).

The results obtained for d_l of Fe i 525.0 nm are summarized below and in Table 3. The standard model corresponds to B = 1600 G, R = 100 km, a filling factor of 0.2, a wall thickness of 20 km, the HSRA in both the QS and the FS interior, and an enhanced UV opacity. The line depth obtained in 3D NLTE in the standard model is printed in boldface in Table 3. The values for the QS entry were taken from a flat 1D HSRA model without any magnetic structure. They correspond well to the values reached far away from the FS. We note that because of Zeeman splitting, the line core depths do not necessarily correspond to the maximum line depth. However, in the averaged profiles, the core intensity is always very close to the minimum intensity. The main conclusions reached from these computations are:

• 1.

  The average 1D profile is always deeper than the corresponding3D profile except for the hotter PLA andNET atmospheres.

• 2.

  The values of d_l obtained for 1D and 3D RT differ by approximately 5−10%. The LTE and 3D line depths may differ by up to 30%.

• 3.

  The atmosphere describing the interior of the FS is – as expected – strongly relevant. The LTE and especially the 1D NLTE line in the hotter NET atmosphere is weaker than that of the true 3D NLTE calculation. We note that the difference in d_l between the FS atmospheres is much smaller when true 3D NLTE is used.

• 4.

  For models with weak magnetic fields (<1500 G), the effects of 3D RT are relatively small owing to the small Wilson depression and the small amount of radiation entering the FT from the walls. For strong fields, the effects increase disproportionately.

• 5.

  The influence of the variation in the geometric quantities (such as the flux sheet radius or the wall thickness, which are both not given in Table 3) on the depth of the average line profile is generally small.

In the past, the thermal structure of the magnetic elements has generally been derived in LTE (e.g. Stenflo 1975; Chapman 1977, 1979; Walton 1987; Solanki 1986; Solanki & Brigljevic 1992). To quantify the error in the temperature of the FS when the proper 3D NLTE formalism was not applied, we tried to find a 1D atmosphere that reproduces the 3D NLTE line depth obtained for the standard FS model. As the 3D NLTE line depth is weaker, we searched for a hotter model than the HSRA. By replacing the HSRA in the FS interior with an atmosphere whose whole temperature curve is shifted up by a constant amount relative to the HSRA, we found that the difference between the depth of the LTE and the 3D NLTE line profiles may lead to an overestimate of the temperature of approximately 50 K for a filling factor of 0.2.

Modern telescopes reach higher spatial resolution and can begin to resolve ever smaller magnetic elements (Lagg et al. 2010). We therefore also considered the error in the temperature obtained from observations that fully resolve the magnetic elements. If we considered the profile emerging along the center of the FS, then this temperature error could reach values as high as 300−400 K (for the standard model). When compared with the 1D NLTE profiles, the difference was only slightly smaller, confirming that 3D RT is relevant as predicted by Stenholm & Stenflo (1978). Similar temperature differences of a few 100 K were also confirmed by inversions (kindly provided to us by Michiel van Noort) of the calculated spectra. As in the vast majority of the cases employed in the past, the inversions were based on LTE forward calculations.

Although this may be an alarming result for those who determine temperatures by comparing observations with LTE calculations, one has to consider that our model is rather simple and may not correspond to a realistic scenario. Thus, we consider it important to investigate this question further using realistic 3D radiation MHD simulations, which we plan to do in an upcoming publication.

Interestingly, the line depth in 3D NLTE for flux sheets whose interiors are described by the PLA and NET atmospheres is, however, rather similar to the line depths obtained in LTE for the same models. Since these models were deduced by employing the LTE approximation, this is a heartening result, since it suggests that these models may also have at least some validity for a 3D diagnostic RT. This might help explain the similarity of the PLA model to the temperature stratification in flux sheets formed in 3D radiation-MHD simulations (Vitas et al. 2009), albeit the RT in the simulations is restricted to LTE.

4.2. Equivalent width and determination of abundances

The equivalent width is a widely used observable for determining the abundances or empirical oscillator strengths. Table 4 summarizes the equivalent widths of the spatially averaged profiles for the same set of calculations and with the same averaging method as in Sect. 4.1. All the values were normalized relative to the 3D NLTE results for the standard model, which is consequently given a value of 100% (given in boldface).

The typical equivalent width of the LTE lines is much larger than that of 3D NLTE calculations (32% in our standard model). The 1D NLTE value is about 5 to 10% larger. The LTE equivalent width is not strongly influenced by any variation in the model parameters because the LTE lines are typically saturated and do not vary strongly across our FS models. Exceptions are the calculations with hotter FS/FT atmospheres (PLA, NET), where the differences between the 3D NLTE and LTE calculation are much smaller. This is not unexpected given the results of Sect. 4.1.

Assuming that the true equivalent width is closer to that obtained in 3D NLTE, then the larger equivalent width obtained using LTE or 1D NLTE can be misinterpreted as a lower iron abundance in the magnetic features or a higher temperature than really present, possibly combined with reduced microturbulence. Sheminova & Solanki (1999) determined elemental abundances in magnetic elements under the assumption of LTE and found similar values as in the quiet Sun. However, they employed empirical model atmospheres for the magnetic elements that were derived in LTE, which likely had an incorrect temperature structure to compensate for the wrong line depths obtained in LTE. The use of these models also gives line equivalent widths close to the observed values for the correct abundances. This implies that a self-consistent LTE computations should more often lead to errors in the determined temperature than the abundances.

By using MHD simulations with magnetic flux typical of moderate active-region plage, Fabbian et al. (2010) found evidence of smaller equivalent widths for the iron lines in the visible, in good agreement with many previous results (e.g. Solanki 1993). This reduced equivalent width is of relevance to the determination of abundances. However, they employed only LTE, so that they basically see the effect of the different temperature stratifications in plage and quiet regions. We expect that the use of 3D NLTE could give rise to a sizeable additional correction.

4.3. Contrast and feature size

One of the interesting results of our calculations is the strong intensity peak near the FS/FT boundary seen, e.g., in Fig. 3. We now investigate the maximum contrast, defined as (I_max − I_QS)/I_QS, where I_QS is the intensity in the quiet Sun and I_max is the maximum intensity reached along the x-axis (generally near the boundary of the FS). We note that we assume an infinite spatial resolution of the observations here (in contrast to Sects. 4.1 and 4.2).

Although we did not find such strong contrasts as Caccin & Severino (1979), our results still show considerable intensity variation within the model. Table 5 gives the line-core contrast for the same model calculations as given in Sects. 4.1 and 4.2. The contrast may be either large (strong fields, especially in NLTE) or small (weak fields, especially for LTE), depending on whether the hot walls are in direct view of the observer. Not only do the LTE and NLTE calculations differ strongly: the 3D and 1D NLTE calculations may differ considerably as well as can be seen from Table 5.

Another interesting finding is the very narrow geometrical width (≈20 km) of these structures as pointed out by Bruls & von der Lühe (2001). If such strong transitions from the FS exterior to its interior exist, a resolution of at least 5−10 km is needed to see the maximum contrast. However, we stress that this test is no proof that such small structures will actually be observable, since horizontal energy transfer via radiation tends to smooth out the horizontal temperature gradients. Thus, when replacing the HSRA model in the FS by the empirically derived NET or PLA models, the large and narrow peaks in the line core intensity at the edge of the FS are either strongly damped or even disappear. This question warrants additional investigation with the help of 3D radiation MHD simulations.

4.4. Stokes profiles and determination of magnetic field strengths

Figures 17 and 18 in Sect. 3.9 reveal that all Stokes parameters are strongly influenced by the calculation method and the input parameters. In Fig. 17, both the strengths of the V/I_c profiles and the profile shapes differ significantly between the LTE and the 3D NLTE case. Both quantities directly influence the determination of the magnetic field strengths. In the present section, we test this using the so-called center of gravity method (Rees & Semel 1979). The results are given in Table 6 for the same models as in Table 5.

For the models with the same atmosphere inside and outside the FS (HSRA/HSRA), the field strength obtained from the 3D NLTE computations is approximately 15% higher than that deduced from the Stokes V profiles computed in LTE, except for the B = 1700 G case, where the relative difference is larger. About half of this difference is due to horizontal RT effects, while the other half comes from NLTE (see the difference in B_cog between LTE and 1D NLTE in Table 6). We note that the Stokes V profiles (Fig. 17) in LTE and 1D NLTE are quite similar but that the Stokes I profiles (Fig. 16) differ strongly and lead to this sizeable difference in B_cog.

Part of the difference in B_cog between LTE and 3D NLTE may be attributed to the different height of formation of the spectral line computed under these different assumptions.

The use of different FS atmospheres also has a strong influence: the field strength deduced from the NET model in LTE reaches almost the value obtained from the 3D NLTE calculation. For the slightly cooler PLA atmosphere, however, the difference in field strength between 3D NLTE and LTE is distinctly higher, although it remains significantly smaller than that found in models with the HSRA/HSRA combination. This result further supports the conclusion that for these two empirical FT models, calculations in LTE provide a reasonable approximation for more realistic 3D NLTE computations.

5. Conclusions

With an improved version of the 3D RH code (Uitenbroek 2000), we have calculated widely used diagnostic iron lines in various flux-sheet (FS) and flux-tube (FT) models. Our adaption of the formal solution to monotonic, parabolic Bézier interpolation of the source function (Auer 2003) is crucial to preventing considerable errors in the calculated intensities and to achieving better convergence. We have developed an enhanced Bézier interpolation scheme for the formal solution of the full Stokes vector problem (for details see Appendices).

We started from a standard expanding flux-sheet model, which consisted of a HSRA atmosphere in both the non-magnetic QS and the magnetic FS, a radius R = 100 km, field strength B = 1600 G at external τ_c = 1, wall thickness w = 20 km, and a filling factor α = 0.2. We studied the behavior of four Fe i spectral lines (524.7 nm, 525.0 nm, 630.1 nm, 630.2 nm) under three different levels of approximation for computing the radiative transfer (full 3D NLTE RT, 1.5D NLTE and 1.5D LTE), each without and with the inclusion of Zeeman splitting (following the field-free method of Rees 1969). We considered the dependence of these lines on FS/FT parameters, such as radius, field strength or Wilson depression, wall thickness, temperature stratification, and UV opacity.

We confirm that 3D RT has a significant effect on spectral line formation, in qualitative agreement with the findings of Stenholm & Stenflo (1977), who used a simple two-level atom in a cylindrical flux-tube model, and Bruls & von der Lühe (2001), who employed 2D calculations in a more realistic expanding flux-sheet geometry. In the following, we list our main conclusions.

The use of a multi-level atom seems to slightly increase the effects found by Stenholm & Stenflo (1977). The inclusion of missing UV opacity in the way proposed by Bruls et al. (1992), however, strongly reduces these effects. Taken together, we found that the influence of 3D RT is much smaller than indicated by the results of Stenholm & Stenflo (1977).

The expanding flux sheet/tube geometry can lead to large contrasts in the line core, particularly at locations where the hot walls of the flux sheet/tube become visible (Caccin & Severino 1979). For almost all the considered models, the NLTE computations produce larger contrasts than LTE, The contrast (and the strength of the 3D effects) increases strongly with increasing field strength (owing to the increasing depth of the Wilson depression).

A striking feature of the emerging line radiation is the strong peak seen at the FS/FT boundary for many of the models and spectral lines, where the flux tube walls are in the direct view of the observer. These peaks (bright rings in FTs) can be a lot narrower than the generally quoted 50 km horizontal mean-free-path for plane-parallel atmospheres even when computed in 3D NLTE. If such sharp boundaries exist in the solar photosphere, than they should remain visible even in the presence of horizontal RT. This result agrees with the findings of Bruls & von der Lühe (2001) and promises a new discovery space for telescopes reaching a spatial resolution as fine as 10 km in the solar photosphere. We note, however, that one reason why the FT or FS walls were so prominently visible is that they represent a sharp horizontal transition from optically thin to optically thick gas. We suspect that the sharp structures that do not display this transition become more smeared out by the radiation.

The resulting intensity of a spherically symmetric FT model, calculated in 3D, shows qualitatively the same effects as those of the computationally cheaper FS model. However, the influence of 3D RT is greater in the FT than the FS owing to the larger wall area radiating into a particular point in the interior of the FT.

Three-dimensional RT generally leads to a weakening of the line in the FS/FT interior. However, the inverse effect also occurs if the interior is hotter than the exterior (at equal geometric height), a situation occurring at certain heights if empirical models, such as the NET or PLA flux-tube atmospheres (Solanki 1986; Solanki & Brigljevic 1992), are used for the FS/FT interior.

The behavior of different lines may vary strongly. The most important factors are the depth of the Wilson depression, which influences the importance of horizontal RT, and the formation height of the line, which determines whether the radiation from the hot walls reaches the absorbing atoms directly.

The average difference between an LTE and a full 3D NLTE RT is on the order of 10% in the line depth and up to 30% in equivalent width. Hence, if LTE is assumed when determining the temperature in the FS/FT interior, significant errors can result in the FT temperature when 3D effects are not taken into account. The higher the spatial resolution, the larger the errors. For spatially resolved magnetic elements, the errors can reach 300−400 K. Abundances may also be underestimated in LTE or 1D NLTE, if correct temperatures are used for the analysis.

Non-LTE and 3D RT effects are strong for all Stokes parameters. The difference between Stokes profiles computed under LTE, 1D NLTE, or 3D NLTE is of similar relative magnitude to that for Stokes I. The application of either LTE or 1D NLTE is thus expected to lead to systematic errors in the determination of the magnetic field strength on the order of from 10% to 20% if the center of gravity method (Rees & Semel 1979) is applied. Since much of this error is introduced by the line being formed at different heights in LTE and NLTE, this difference in field strength is equivalent to an error in the height to which a given field strength is assigned.

By good fortune the empirical PLA and NET models turn out to give rather similar results in LTE as in 3D NLTE, so that the validity of results obtained with these models is broader than previously thought.

In summary, the present investigation provides evidence that the interpretation of observations of magnetic flux concentrations may change when 3D NLTE RT effects are taken into account. Our quantitative estimates can be considered a rough guide only, because of the simplicity of the model employed here. Hence, the natural next step of this line of investigation is to employ 3D MHD simulations and investigate commonly used diagnostics in both LTE and 3D NLTE. Such an investigation is planned.

Online material

Appendix A: Formal solution: influence of the source function interpolation method

In the short characteristics (SC) method (Olson & Kunasz 1987; Kunasz & Auer 1988), the formal solution is integrated only along the “short” path from the nearest up- and downwind points. The incoming intensities – as well as all other quantities used for the formal integration at the up- and downwind points – have usually to be interpolated from neighboring grid points. The problem of polynomial interpolation at the up- and downwind points, which may lead to spurious extrema, was widely discussed in the literature until Auer & Paletou (1994) found an elegant and efficient solution using monotonic interpolation.

In contrast, one finds very few references to the problems that arise when the source function is interpolated with a parabolic Lagrange polynomial along the short characteristic as proposed and established as standard by Auer et al. (1994). With S(τ) = S₀ + c₁τ + c₂τ², the integral over S along the short characteristics is then given by the analytical expression ${\mathit{I}}_{\mathrm{local}}\mathrm{=}{\mathit{c}}_{\mathrm{0}}{\mathit{S}}_{\mathrm{0}}\mathrm{+}{\mathit{c}}_{\mathrm{1}}{\mathit{w}}_{\mathrm{1}}\mathrm{+}{\mathit{c}}_{\mathrm{2}}{\mathit{w}}_{\mathrm{2}}\mathit{,}$(A.1)where I_local is the locally produced intensity without irradiation from outside. The w_i may be – in the case of Lagrange interpolation – expressed as simple functions of known quantities (see Auer & Paletou 1994). The first indication of these interpolation problems is given by Auer (2003), who proposes monotonic Bézier interpolation instead of simple Lagrange interpolation to prevent spurious extrema of S along the integration path.

To demonstrate that the proper interpolation of S is as important as the correct interpolation of up- and downwind quantities, we present the situation of an overshooting source function, which we call “a critical case”, taken from our fluxsheet atmosphere simulation (the standard model, see Sect. 2.3 for a description). Figure A.1 shows the values of S at the up- and downwind points (S₁,   S₂) and the integration point (S₀). The unique Lagrange parabola through the three base points is given by the solid line. One immediately sees that S clearly takes values below the base points S₀ and S₁, while it should decrease monotonically from S₀ to S₁. From ${I}\hspace{0.17em}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\hspace{0.17em}{I}\mathrm{1}{\mathrm{e}}^{\mathrm{-}\mathrm{\Delta }{\mathit{\tau }}_{\mathrm{1}}}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\underset{\mathrm{0}}{\overset{{\mathit{\tau }}_{\mathrm{1}}}{\mathrm{\int }}}\mathit{S}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }}\mathrm{d}\mathit{\tau ,}$(A.2)it directly follows that the resulting intensity must be wrong (i.e. too small), because the area under S is far too small. We note that the factor e^−τ does not change the qualitative result as it is on the order of unity over the relevant integration interval.

Bézier parabolas are defined by two base point values and a single control point, which is always located at the intersection of the two tangents of S at the base points. This control point therefore “controls” the run of S between the two base points. The unique Lagrange parabola through S₀, S₁, and S₂ is recovered in the interval from S₀ to S₁ when the control point (c₀ in Fig. A.1) lies on the corresponding tangents of the Lagrange parabola (fine dotted). One imposes monotonicity in the interval from S₀ to S₁ by constraining the control point to values between the two base point values of S. By moving the control point up to the lower of the two values (c′₀), one gets a new parabola (that no longer passes through S₂!) with a monotonic run in the integration interval (dashed curve). Therefore, the area under S is no longer artificially reduced. The difference in the locally produced intensity for our sample, taken from an existing calculation, is considerable.

We therefore used monotonic parabolic Bézier interpolation as proposed by Auer (2003) for all our calculations. We performed tests with a snapshot taken from a realistic 3D radiation hydrodynamic simulation run (i.e. for B = 0) with the MURAM code (Vögler et al. 2005). Additionally, we used a 3D version of the FAL-C model of Fontenla et al. (1993) obtained by arranging many identical 1D atmospheres into a 3D cube. In both cases, we found that about 3% of all formal solutions suffer from S overshooting in the integration interval. While the resulting intensity for the FAL-C model is almost the same as for the Lagrange interpolation (most of the critical cases in FAL-C display only quite weak overshooting because of the smoothness of the atmosphere), the resulting intensities in the snapshot may differ considerably owing to the strongly varying atmospheric parameters along the integration path. Therefore, the source function may exhibit strong overshooting.

An increasing number of publications have recognized this problem in recent years. For instance, Hauschildt & Baron (2006) used linear interpolation (instead of monotonic Bézier) for critical cases and Koesterke et al. (2008), and Hayek (2008) used cubic Bézier polynomials.

In more extreme atmospheres than our fluxsheet model (as found e.g. in realistic MHD simulations, which will be discussed in future publications), the resulting intensity may even be wrong by orders of magnitude with disastrous consequences not only for the final result, but also for the convergence of the RT computations themselves. The convergence is affected because the (usually diagonal) approximate Lambda operator used for the improvement of the iteration is generally interpolated in the same manner as the source function.

In our experience, we have found that the monotonic, parabolic Bézier interpolation, which is identical to Lagrange interpolation when S does not overshoot in the integration interval, has much better convergence properties than the pure Lagrange interpolation, which can lead to a complete divergence for non-smooth atmospheres, such as those along rays passing through the walls of a flux tube.

However, imposing either monotonicity or a special treatment in critical cases may lead to convergence problems as well because a change in the interpolation method at a certain point in the cube may influence neighboring points in such a way, that in the next iteration the interpolation method is changed back again to the original method. This may lead to an infinite flip-flop situation, so that convergence is never achieved. We encountered these flip-flop situations in a few of our unsmoothed flux-sheet models. They occurred at a convergence level of a few percent in the population numbers and were much stronger when linear interpolation was used in the critical cases (a change from linear to parabolic may have more influence on the area under S than a change from monotonic Bézier to normal Lagrange parabolas).

The problem of strongly overshooting source functions may also arise in the 1D schemes as they interpolate S in the integration interval with a Lagrange polynomial as well. In addition, we therefore implemented and used monotonic Bézier interpolation for our 1D calculations.

Appendix B: Formal solution for the full Stokes vector problem

The solution to the formal integration of the RT equation for the full Stokes vector within the short characteristics method of Olson & Kunasz (1987) and Kunasz & Auer (1988) was given by Rees et al. (1989). Their second method, which is based on a diagonal-element lambda operator (DELO) uses linear interpolation of the source function components in the integration interval. Socas-Navarro et al. (2000) improved the accuracy of the source function integration by introducing parabolic interpolation. The RH code is based on this improved method.

As mentioned in the previous section, the interpolation with a Lagrange parabola may cause great difficulty in the formal solution. However, these problems are not limited only to Stokes I. The overshooting of the parabola may affect the other Stokes components in the same way. We therefore used the monotonic Bézier interpolation for all components.

Unfortunately, the full Stokes vector problem is even more intricate: the integration with a monotonic, Bézier-interpolated source function may still lead to unphysical results in special (but not uncommon) cases. Figure B.1 illustrates a pathological situation where the integration of the source functions leads to the unphysical result of V > I. Owing to the different (in this case even opposite) curvatures of S_I and S_V (the latter being a representative of all Stokes components), S_V is larger than S_I for a large part of the integration interval (τ₀ to τ₁). The likelihood of this occurring is highest when S_V is almost as large as S_I at one point and much smaller (or negative) at neighboring point(s).

We note that the problem arises because one always has to consider a certain neighborhood of the integration interval when determining the interpolating function (except for linear interpolation where the two end points are sufficient). This very general result of interpolation theory is also responsible for the overshooting of S_I in the case of parabolic Lagrange interpolation. If the neighborhood does not behave “well”, it may produce unrealistic behavior (i.e. oscillations or overshooting) in the integration interval.

The Bézier parabola helps to reduce the influence of these difficult neighborhoods as it is defined in a more local way, i.e. by the two end points and a local control point, thus limiting the influence of points outside the integration interval. We note that the control point is influenced by the neighborhood (i.e. S₂), but in a less strict way. By imposing monotonicity (i.e. the control point has to be between the two values of the endpoints), the influence of the neighborhood is reduced exactly for these cases where overshooting would otherwise take place.

As all Stokes components contribute to Stokes I, it is reasonable to take the shape of S_I(τ) as a reference. If we replace ${\mathit{c}}_{\mathrm{0}}^{\mathrm{\left(}\mathit{V}\mathrm{\right)}}$ by ${\mathit{c}}_{\mathrm{0}\mathit{,}\mathrm{new}}^{\mathrm{\left(}\mathit{V}\mathrm{\right)}}$ in such a way that the shape of S_I(τ) predefines the shape of S_V(τ) along the integration interval, the unphysical condition of S_V > S_I may be prevented. For Bézier parabolas, this can easily be achieved with $\frac{{\mathit{S}}_{\mathrm{0}}^{\mathrm{\left(}\mathrm{QUV}\mathrm{\right)}}\mathrm{-}{\mathit{c}}_{\mathrm{0}\mathit{,}\mathrm{new}}^{\mathrm{\left(}\mathrm{QUV}\mathrm{\right)}}}{{\mathit{c}}_{\mathrm{0}\mathit{,}\mathrm{new}}^{\mathrm{\left(}\mathrm{QUV}\mathrm{\right)}}\mathrm{-}{\mathit{S}}_{\mathrm{1}}^{\mathrm{\left(}\mathrm{QUV}\mathrm{\right)}}}\mathrm{=}\frac{{\mathit{S}}_{\mathrm{0}}^{\mathrm{\left(}\mathit{I}\mathrm{\right)}}\mathrm{-}{\mathit{c}}_{\mathrm{0}}^{\mathrm{\left(}\mathit{I}\mathrm{\right)}}}{{\mathit{c}}_{\mathrm{0}}^{\mathrm{\left(}\mathit{I}\mathrm{\right)}}\mathrm{-}{\mathit{S}}_{\mathrm{1}}^{\mathrm{\left(}\mathit{I}\mathrm{\right)}}}\mathrm{·}$(B.1)In this way, we use only local information, i.e. the shape of S_I(τ) and the endpoint values of S_QUV, to define the run of S_QUV(τ).

We note that in the DELO method of Rees et al. (1989), the integration of the source functions is only part of the solution. The determination of the full Stokes vector intensities requires the solution of a system of linear equations. Even with our integration method, in very rare cases, the solution of the equation system produces unphysical results with one or more Stokes components becoming larger than Stokes I. Consequently, we changed the RH code in such a way, that the polarized Stokes component fulfills the condition $\sqrt{{\mathit{Q}}^{\mathrm{2}}\mathrm{+}{\mathit{U}}^{\mathrm{2}}\mathrm{+}{\mathit{V}}^{\mathrm{2}}}\mathrm{\le }\mathit{I}$(B.2)

after each formal solution. In all of these cases in which the relation was not fulfilled, the polarized components were reduced by a single factor f (i.e. Q_new = fQ, U_new = fU, V_new = fV) such that the equality in the equation above held, i.e. $\sqrt{{{\mathit{Q}}_{\mathrm{new}}}^{\mathrm{2}}\mathrm{+}{{\mathit{U}}_{\mathrm{new}}}^{\mathrm{2}}\mathrm{+}{{\mathit{V}}_{\mathrm{new}}}^{\mathrm{2}}}\mathrm{=}\mathit{I}\mathit{.}$(B.3)

Acknowledgments

We are grateful to Han Uitenbroek for providing his excellent full Stokes NLTE radiative transfer code and helpful discussions. We thank Jo Bruls for providing the iron model atom and many insights into NLTE effects. We thank Michiel van Noort for inverting some of our calculated profiles. R. H. appreciates the flexibility of Prof. Dr. Norbert Dillier concerning the working hours at the University Hospital of Zurich. This work has been partially supported by WCU grant No. R31-10016 funded by the Korean Ministry of Education, Science and Technology.

References

All Tables

All Figures