© ESO 2021

1. Introduction

As instrumentation improves, inversions of the increasingly large numbers of polarized Stokes profiles measured in observations of small-scale magnetic fields in the solar photosphere continues to be fraught with risks, limitations, and challenges (see Bellot Rubio & Orozco Suárez 2019 for a full review). At spatial resolutions currently achievable, asymmetric and multi-lobed polarized Stokes vectors are common (Khomenko et al. 2003; Kiess et al. 2018) in the internetwork (IN), indicating that we are not resolving the magnetic structures. To accurately reproduce these profiles, multiple components or gradients in optical depth are necessary, with both solutions providing very similar fits to the observations (Martínez González et al. 2016). To retrieve accurate inclinations, one must measure both linear polarization and circular polarization above the noise level. As the amplitudes of these signals typically are very close to this noise level, the signal-to-noise ratio (S/N) has a huge impact on the distribution of inclinations retrieved (Asensio Ramos 2009). Even more challenging is determining the azimuthal angle of the magnetic vector, which requires both Stokes Q and U above the noise level and utilization of disambiguation techniques (Metcalf et al. 2006). That is, to fully constrain the magnetism of the IN, full-vector spectropolarimetry is required.

Recent observations by Campbell et al. (2021) with the GREGOR infrared spectrograph integral field unit (GRIS-IFU) instrument (Collados et al. 2012) demonstrate that even with a very highly magnetically sensitive spectral line in the near infrared, observed with the largest solar telescope in Europe (Schmidt et al. 2012; Kleint et al. 2020), it is possible to observe linear polarization in one frame, only for it to vanish in the next with a relatively short cadence. Therefore, with IFUs, a balance must be achieved between total integration time, field of view (FOV) size, and cadence. These instrumental limitations are further compounded by Zeeman mixing induced by insufficient spatial and spectral resolution or, indeed, the evolution of the target during exposures. Observing the dynamics of small-scale magnetic fields has previously proved very challenging (Wiegelmann et al. 2013; Campbell et al. 2021; Gošić et al. 2021). It remains a key science goal to reveal and study the dynamics of magnetism in increasingly larger areas of the IN, as the emergence of flux in these regions can account for the transport of a large amount of energy to the solar surface (Trujillo Bueno et al. 2004). Previous studies with very large FOVs have shown there are enormous regions of the IN apparently devoid of magnetic flux that can be revealed as harbouring a significant amount of transverse magnetic flux with sufficient S/N and spatial resolution (Lites et al. 2008) or, indeed, with the Hanle effect (Stenflo 1982).

One method for ensuring results derived from observations are well constrained is to compare them to semi-empirical models (Borrero & Kobel 2011) or realistic magnetohydrodynamic (MHD) simulations (Vögler et al. 2005; Khomenko et al. 2005; Khomenko & Collados 2007; Lagg et al. 2016; Lites et al. 2017; Milić et al. 2019; Smitha et al. 2020). With the advent of next generation high resolution telescopes, our understanding of how the magnetic field is organized in the IN photosphere is likely to experience a significant leap as observations reveal these regions of the solar atmosphere at unprecedented spatial resolution. In this study, we use a snapshot of a realistic MURaM simulation to evaluate the extent to which we can retrieve accurate information about the magnetic vector in the IN photosphere using inversion techniques. We synthesize the full Stokes vector by solving radiative transfer under the assumption of local thermodynamic equilibrium (i.e. the forward problem, see del Toro Iniesta & Ruiz Cobo 2016 for a full review), only to then invert the synthetic spectra to establish the extent to which the original atmosphere differs according to the inversions. We then investigate how this model atmosphere may appear had it been observed at the spatial resolutions achievable by the GREGOR/GRIS-IFU. Similar processes have been employed to determine the effects of instrumentation on the analysis of wave phenomena (Keys et al. 2021). We critically examine the results reported by Campbell et al. (2021) at this resolution, and also seek to help guide the approach taken by future observers considering this problem at higher resolution (e.g. with the Daniel K. Inouye Solar Telescope (DKIST); Rimmele et al. 2020).

2. MHD cube and synthesis

We use a snapshot of a realistic 3D MHD model simulation of solar magnetoconvection (Vögler et al. 2005) for synthesizing the Stokes profiles. The MHD cube, described in detail by Nelson et al. (2013), has a spatial resolution of 25 km in the plane of the solar surface and 14 km in depth. The cube extends 12 Mm × 12 Mm in the x and y directions and 1.4 Mm in the z direction, or 480 × 480 × 100 grid points, respectively. The simulation started with a plane parallel solar atmosphere and then purely hydrodynamical convection was allowed to develop until it reached a steady state. A ‘checkerboard’, mixed polarity magnetic field (of magnetic field strength, B, of ±50 G) was introduced. The snapshot was taken 30 min after this point, allowing the magnetic field to disperse and the cancellation of opposite polarities resulted in a reduction of magnetic energy that resembles the conditions of the photospheric IN. The result is a cube with an abundance of mixed polarities, mostly concentrated in the intergranular lanes (IGLs).

2.1. Transformation of the atmosphere to an optical depth scale and line synthesis

We use the SIR code (Stokes inversion based on Response Functions; Ruiz Cobo & del Toro Iniesta 1992) for synthesizing and inverting the full Stokes vector. We assume disk centre observations for all the studies we perform in this work (i.e. μ = 1, where μ = cos(θ) and θ is the heliocentric angle). We use the abundance values of the different atomic species given in Asplund et al. (2009). We use the same atomic data as used by Campbell et al. (2021) and invert the same spectral lines. We do not degrade the Stokes profiles during the synthesis with any macroturbulence correction, and the microturbulence is null as in the original MURaM snapshot.

We transform the 3D cube, originally in gas pressure and a geometrical height scale, to an atmosphere in optical depth using the subroutines found inside the SIR source code distribution. The transformation is a requirement of the SIR code for doing synthesis and inversions. To evaluate the optical depth at 5000 Å from geometrical height and gas pressure we integrate the equation

$$\begin{array}{c}\hfill \mathrm{d}{\tau }_{5000\mathrm{\AA }}=-\rho k_{5000\mathrm{\AA }}\mathrm{d}z,\end{array}$$(1)

where ρ is the density, k_5000 Å is the continuum opacity at 5000 Å, and z is the geometrical height. The k_5000 Å is evaluated following the recipe of Wittmann (1974). The main difficulty of this computation is that we need to evaluate all the opacity contributions to k_5000 Å and for that we need to know the density of each ion that has an important contribution to the continuum opacity (e.g., H⁻, H, He⁻, He, ${\mathrm{H}}_2^{-}$, ${\mathrm{H}}_2^{+}$, and metals such as Cl⁻, C, Na, Mg). Those densities are computed by solving the Saha equation under the local thermodynamic equilibrium (LTE) approximation of SIR. However, the Saha equation needs, as input, the electron density to compute the atomic populations of the different ions, but the electron density depends on the ionization rate (i.e. on the results of the Saha equation for every ion). Thus, the problem needs to be solved iteratively to achieve a solution that is self-consistent. SIR does this using an algorithm proposed by Mihalas (1967), that is also described in Wittmann (1974).

2.2. Degradation and Comparison with GREGOR observations

In order to compare real observations to our synthetic profiles produced from the MHD cube, we first have to account for instrumental effects by degrading the spatial and spectral resolution of the synthetic profiles. We followed a similar procedure as Milić et al. (2019), but with the added benefit of having real observations with which to compare our synthetic spectra. First, we convolved the synthetic unpolarized and polarized spectra with the spatial point spread function (PSF) of the GREGOR telescope. To estimate this PSF, we first considered the root-mean-square continuum intensity contrast, defined as

$$\begin{array}{c}\hfill \delta I_{\mathrm{rms}}\approx \frac{\sqrt{\frac{1}{N}\sum _{x,y}(I(x,y)-⟨I_{x,y}⟩{)}^2}}{⟨I_{x,y}⟩},\end{array}$$(2)

where I(x, y) are continuum intensity values of individual pixels and ⟨I_x, y⟩ is the mean over all pixels, N. The peak δI_rms recorded by Campbell et al. (2021) with the GREGOR/GRIS-IFU was 3.1% on 6 May 2019 with exceptional atmospheric seeing conditions. We used a least-squares procedure to iteratively convolve the synthetic cube with a two-dimensional Lorentzian kernel of variable scale parameter defining the width of the kernel, Γ, in Fourier space until the δI_rms of the convolved cube reached 3.1%. The Γ value was determined as between 4−5 pixels in the MHD cube, making the full width at half maximum (FWHM) about 0.53″ or 380 km on the solar surface. We then re-binned the degraded cube to a spatial sampling of 100 km/pixel, which is approximately equivalent to the 0.135″/pixel sampling of the GRIS-IFU in the x-direction. We refer to these degraded profiles as version 1 (V1). The theoretical diffraction limited resolution of GREGOR is 0.27″ in this near infrared spectral window, but observations have not achieved this spatial resolution – 0.4″ was reported by Lagg et al. (2016), but with a lower δI_rms. We therefore followed the method described by Milić et al. (2019) to consider the impact of stray light, which further degrades the δI_rms. Here we are defining stray light as an additional unpolarized background added to the synthetic profiles, which the adaptive optics cannot correct. We added a very broad Gaussian to our spatial PSF, so that it becomes a Voigt profile. We model our spatial PSF as:

$$\begin{array}{c}\hfill I^{\mathrm{stray}}(x,y)=(1-k)I^{\mathrm{spat}}(x,y)+k[I^{\mathrm{spat}}(x,y)\ast G({\sigma }_{\mathrm{stray}})],\end{array}$$(3)

where I^stray(x, y) is the spatially degraded and re-sampled cube of synthetic profiles, I^spat(x, y), further degraded by stray light, k is the stray light fraction and G(σ_stray) is a Gaussian with standard deviation, σ_stray, of 2″. Based on measurements by Borrero et al. (2016) we used a k value of 0.3. We then convolved the Voigt spatial PSF as before, keeping the G(σ_stray) fixed and adjusting the Γ value of the Lorentzian until we achieved a δI_rms of 3.1% again. With this approach, the FWHM is measured as 0.3″, much closer to the diffraction limited resolution of GREGOR. We refer to the synthetic profiles degraded (and re-sampled) at every wavelength with only the Lorentzian with FWHM of 0.3″ as version 2 (V2), and the synthetic profiles degraded by the Voigt function, which additionally includes G(σ_stray), as version 3 (V3). We also refer to the undegraded synthetic spectra as version 0 (V0).

Although we produced the synthetic cube with the same dispersion as GREGOR/GRIS, we still need to consider the spectral resolution. The spectral PSF of GREGOR/GRIS has not been empirically measured, but estimations have been made with close agreement (Borrero et al. 2016; Campbell et al. 2021). We convolved the synthetic cube in the wavelength domain with a Gaussian of FWHM 172 mÅ. This also significantly reduced the amplitude of the polarized Stokes profiles. We applied this spectral degradation to V1, V2, and V3, and henceforth when we refer to any version (except V0) the inclusion of spectral degradation is implied. Finally, we added photon noise to each Stokes vector in wavelength space. At each spectral point, we added a value randomly derived from a normal distribution with a standard deviation of σ_n = 8 × 10⁻⁴I_c. This is the same as the noise level in Stokes Q achieved by the GRIS-IFU in recent observations (Campbell et al. 2021) and this will allow us to compare their results with the MHD simulations. We also separately considered a degraded cube with σ_n = 3 × 10⁻⁴I_c, which is similar to the noise level achieved by Hinode/SP in ‘sit-and-stare’ mode (Lites et al. 2008) and slightly below the noise level achieved by GRIS (Lagg et al. 2016). The latter σ_n value is meaningful as it allowed Lites et al. (2008) to reveal hidden magnetism in the IN, and we have calculated this as the lowest noise level the GRIS-IFU could achieve without the cadence between frames becoming too large to be sensible. This value should also be achievable by the IFU at DKIST. In summary, V0 is the undegraded synthetic profiles. V1 is the synthetic profiles degraded to a spatial resolution of 0.53″, re-sampled to 0.135″/pixel in x and y and spectrally degraded to 172 mÅ. V2 is the synthetic profiles degraded to a spatial resolution of 0.3″, spatially re-sampled and spectrally degraded in the same manner. Finally, V3 is the same as V2 but with additional stray light degradation.

Figure 1 demonstrates the impact of the spatial PSF and stray light on the synthetic intensities. Much of the spatial detail differentiating IGLs from granules is lost. The distribution of continuum normalized intensities at GREGOR resolution is consistent with the distribution recorded by the GRIS-IFU for both V1 and V3. At the same time, V2 shows how if one wants to accurately model the spatial PSF based on measurements of δI_rms in ground-based observations, the spatial resolution will be dramatically higher if the impact of stray light is accounted for. This is important as we assume the stray light is unpolarized, meaning that if we neglect it we will under-estimate the spatial resolution and degrade the polarization signals by too large a degree, significantly reducing their amplitude while smearing the polarization over greater areas. Before instrumental degradation, the synthesized cube has 77.5% of profiles with Stokes V signals and 19.5% with either Stokes Q or U displaying a maximum value in the 15648.52 Å line greater than 3.6 × 10⁻³I_c (i.e. the 4.5σ_n level of the GRIS-IFU observations). Additionally, 18.7% of profiles have a maximum in either Stokes Q or U and Stokes V simultaneously above this amplitude. These values are reported in Table 1, along with the values for the various degraded versions at amplitudes of 3.6 × 10⁻³I_c and 2.4 × 10⁻³I_c (the latter being the 3σ_n level of the GRIS-IFU observations). The statistics for V2 and V3 are reasonably similar to that reported by Campbell et al. (2021) in terms of linear polarization, however with the synthetic profiles we have much more circular polarization than is measured in observations.

Fig. 1. Synthetic data produced from synthesis of MHD simulation outputs using SIR, shown at a continuum wavelength of 1565.47 nm. The top left panel shows the original Stokes I synthetic map (V0), while the top right and middle left panels show the spatially degraded map at an effective resolution of 380 km (V1) and 215 km (V2). The middle right panel shows the same as V2 but further degraded by stray light (V3). The bottom panel shows a histogram of the synthetic continuum intensities in each case, in addition to the distribution from real GREGOR/GRIS-IFU observations from 6 May 2019.

3. Inversions

3.1. Inversions before instrumental degradation

In order to test how well we can infer the thermodynamic, kinematic, and magnetic properties of the solar atmosphere we inverted the synthetic spectra with the least-squares SIR code (Ruiz Cobo & del Toro Iniesta 1992). We adopt a similar Monte Carlo approach to Campbell et al. (2021) in that we repeat the inversion procedure 30 times per pixel with randomized model parameters in the inclination of the magnetic field with respect to the observer’s line-of-sight (LOS), γ, the azimuth of the magnetic field in the plane perpendicular to the observer’s LOS, ϕ, the modulus of the magnetic field, B, and the line-of-sight velocity, v_LOS. We followed this strategy to reduce as much as possible the probability of converging at a solution that is located in a local, as opposed to global, minimum in the χ² hyper-surface, and this also reduces the influence of the initial guess model on the final solution. For the remaining atmospheric parameters, namely temperature, T, electron pressure, gas density, and gas pressure, we initialised the inversion code with an initial FAL-C quiet Sun model (Fontenla et al. 2006), which describes each of the parameters characterizing the stratification of the atmosphere in the log of the optical depth at 5000 Å, log τ_5000 Å. We assume the magnetic filling factor, α, is 1, because in MURaM each pixel has no unresolved magnetic components. In this way, we are performing the experiment under the assumption that the resolution of the MHD cube represents the absolute scale of the Sun’s atmosphere. We use the parallelized wrapper to SIR developed and made available by Gafeira et al. (2021). By inverting multiple lines, with varying strengths, opacities, and effective Landé g-factors, we can probe the atmosphere more reliably than if we inverted only one line. Table 2 summarises the nodes used in the first inversion scheme (scheme C), which is a simple Milne-Eddington-like inversion, allowing us to retrieve depth-averaged parameters from the atmosphere. We chose 4 nodes in T. In principle, the exact number of nodes in T must be high enough to match the realistic stratification of the quiet Sun but low enough to avoid oscillatory solutions being introduced, which we do not observe in any cases in our results. We do not allow the microturbulent and macroturbulent velocities, v_mic and v_mac, which account for the presence of small-scale and large-scale turbulent motions, respectively, to vary as free parameters. In the second inversion scheme (scheme G), we increase the nodes in B, γ, and v_LOS to 2, to evaluate both the extent to which we can retrieve accurate atmospheric parameters at more than one optical depth and the extent to which the inclusion of a simple linear gradient can improve the χ² values achieved. Using response functions, presented by Campbell et al. (2021) and Quintero Noda et al. (2021), it is expected that the greatest response to B, γ, and v_LOS is to be found at optical depths between log τ_5000 Å = −1.0 and −0.5. Meanwhile, the maximum response to T is expected to be found at k_5000 Å = 0.5.

Figure 2 shows, for the synthesized and inverted profiles under scheme C, the wavelength-integrated, net linear (LP) and circular (CP) polarization, which we define as follows:

$$\begin{array}{c}\hfill \text{LP}=\frac{{\int }_{{\lambda }_b}^{{\lambda }_r}{[Q^2(\lambda )+U^2(\lambda )]}^{\frac{1}{2}}\mathrm{d}\lambda }{I_{\mathrm{c}}{\int }_{{\lambda }_b}^{{\lambda }_r}\mathrm{d}\lambda },\end{array}$$(4)

$$\begin{array}{c}\hfill \text{CP}=\text{sgn}(V_{\text{b}})\frac{{\int }_{{\lambda }_b}^{{\lambda }_r}|V(\lambda )|\mathrm{d}\lambda }{I_{\mathrm{c}}{\int }_{{\lambda }_b}^{{\lambda }_r}\mathrm{d}\lambda },\end{array}$$(5)

Fig. 2. Comparison of the synthesized (left) and inverted (middle) maps of the LP (upper row) and CP (lower row) according to scheme G. For both LP and CP the synthesized and inverted parameters are plotted against each other with the colour scale on a log10 scale (right). LP and CP, defined in Eqs. (4) and (5), respectively, are computed by integrating across the 15 648.52 Å line. The solid, green line indicates the exact linear relationship. The dashed, green ellipse indicates the 99% confidence level.

where sgn(V_b) represents the sign of the blue Stokes V lobe and λ_r and λ_b are the red and blue wavelength limits of integration, respectively. It is clear that there is strong correlation for the each of these parameters between the synthetic and inverted profiles. The LP and CP are highly sensitive to the shape and amplitude of the lobes in Stokes Q and U in the former case and Stokes V in the latter case. As scheme C can only produce symmetrical polarized Stokes profiles, a small difference is therefore expected.

Figure 3 shows various atmospheric parameters for the synthesized and inverted profiles under scheme C and G. Values of T and γ compare very closely for both schemes, whereas values of B and v_LOS appear slightly under-estimated by the inversion in scheme C. As is clear from the maps for all parameters, the topology of the parameters is well retrieved by the inversion. Figure 4 further demonstrates this, showing the distributions of B, γ, and v_LOS at log τ_5000 Å = −0.5 for the input MHD atmosphere and those values returned by scheme C and G inversions. T is also shown but at log τ_5000 Å = 0.5. Both C and G inversions reproduce the B distribution well, but scheme C produces a γ distribution that more closely resembles the MHD cube. In order to quantify how well constrained the parameters are, we calculate the percentage difference for each parameter. For temperature, we define the percentage difference, ΔT, as

$$\begin{array}{c}\hfill \mathrm{\Delta }T=\frac{T^{\mathrm{MHD}}-T^{\mathrm{SIR}}}{T^{\mathrm{MHD}}},\end{array}$$(6)

Fig. 3. Comparison of the true MHD (first column) and inverted (second column) maps of the T (top row), B (second row), γ (third row), and vLOS (bottom row). The maps of inverted parameters are shown as determined by inversion scheme G (except for γ, which is shown as determined by scheme C). The colour table for B is extended higher than the maximum value (1915 G) to accommodate those with colour vision deficiency. The MHD and inverted parameters, for both scheme C (third column) and G (fourth column) are plotted against each other with the colour scale on a log10 scale. The solid, green line indicates the exact linear relationship. The dashed, green ellipse indicates the 99% confidence level. T is shown at log τ5000 Å = 0.5, while the other parameters are shown at log τ5000 Å = −0.5.

where T^MHD is the value from the MHD cube transformed to the optical depth scale and T^SIR is the value returned by SIR under scheme C or G. The values of T^MHD and T^SIR are retrieved at the same optical depth. We similarly define ΔB, Δγ, and Δv_LOS. Positive values indicate that SIR has under-estimated the parameter, while negative values indicate an over-estimate. Figure 4 also shows the percentage differences for both scheme C and G. It is important to emphasize that because the colour tables in Figs. 3 and 4 are on a log₁₀ scale, the vast majority of pixels are contained in the brightest bins in the histograms. Additionally, in Table 3, we calculate the mean, standard deviation, and median of ΔT, ΔB, Δγ, and Δv_LOS across the populations in each case. When calculating ΔB and Δγ we only consider pixels with a polarized Stokes vector with a maximum amplitude greater than 2.4 × 10⁻³I_c to restrict analysis to magnetic fields that would realistically be detectable. For Δv_LOS, we consider only pixels with values > 100 ms⁻¹ to avoid very small values entering into the denominator, but 97% of the FOV meets this criteria. In Fig. 4, for each parameter, the median is overplotted as incrementally calculated across the populations in the same bins as the distributions are shown in the histograms. There is a small difference in percentage differences returned by the inversion schemes for T, but as this parameter had the same number of nodes in both schemes the improvement in scheme G is likely a by-product of the other parameters being better constrained. For γ, the values appear slightly better constrained by scheme C. It is further evident from the median values of ΔB and Δv_LOS in Table 3 that scheme C under-estimates these parameters on average, but the reduction in the medians and the shifting of the distributions in Fig. 4 shows that scheme G resolves this issue, although there is still a spread of values. For B in particular, stronger fields are better constrained than those in the weak field regime. As Fig. 5 illustrates, scheme G produces consistently lower χ² values than scheme C particularly in IGLs and in IGL-granule boundaries, where v_LOS and B gradients are likely to be present.

Fig. 4. Distributions of values (left) and histograms between the ‘true’ MHD values and the differences between ‘true’ MHD and inverted values (middle and right, for scheme C and G, respectively) in B, γ, vLOS, and T (upper to lower row, respectively). Values are retrieved at log τ5000 Å = −0.5, except for T which is retrieved at log τ5000 Å = 0.5. In the distributions, values from the MHD cube (orange lines) and returned from the inversions under scheme C (red lines) and G (blue lines) are shown. The histograms of B, γ, ΔB, and Δγ are shown for pixels with maximum amplitude in at least one polarized Stokes vector > 2.4 × 10−3Ic. The median error as a function of each parameter is shown by the green, dashed lines. The percentages in all histograms are weighted with respect to the total number of pixels, N, i.e. including those with no polarization. The colour scale is shown on a log10 scale.

Fig. 5. Map of the normalized difference between χ2 values produces by scheme C and scheme G inversions.

We have established that the inclusion of simple gradients in the inversion schemes can improve the χ² statistically, but in order to examine the extent to which we can retrieve accurate atmospheric parameters at more than one optical depth we calculated the Pearson correlations, P_c, between the MHD values and the values returned by scheme G inversions. First, the values returned for scheme C were 1.00, 0.99, 0.98, and 0.90 for T, v_LOS, B, and γ, respectively. Figure 6 shows the P_c for scheme G at a range of optical depths. Of course, the results resemble very clearly the response functions, but we are explicitly testing the ability of the SIR code to retrieve the atmospheric parameters from realistic 3D MHD model outputs. The results confirm the potential for these lines to constrain the thermodynamic and kinematic parameters deep in the atmosphere, as also reported by Milić et al. (2019), with T responsive to a broader range of optical depths than v_LOS. In the deepest layers, the P_c is highest in T as expected, since local thermodynamic equilibrium is firmly established in these layers (del Toro Iniesta 2003). T also has two peaks in P_c, one very deep in the atmosphere, where Stokes I is very responsive to T in continuum wavelengths, and one higher, where wavelengths with spectral lines carry information about T. On the other hand, the magnetic parameters do not appear to be as responsive to the deepest layers. For B, the P_c is highest between log τ_5000 Å = −1.5 and −0.5, but the correlation is also significant as high in the atmosphere as log τ_5000 Å = −2.0 and as deep as log τ_5000 Å = 0.0. The B values are less strongly correlated at log τ_5000 Å = −2.0 and at higher layers. For γ, the P_c is highest between log τ_5000 Å = −1.0 and −0.5, but the correlation is weaker than all other atmospheric parameters, highlighting the difficulty of constraining γ in quiet Sun scenarios. It must be pointed out that these correlations could potentially be improved if nodes could be placed at specific optical depths based on an analysis of the response functions (as has been investigated by Socas-Navarro 2004) and associated errors. This functionality is not yet available with SIR but will be in the future (private communication with the developer).

Fig. 6. Pearson correlation coefficients for values of B, γ, vLOS, and T between the MHD cube and scheme G estimates at a range of optical depths. The colour scale saturates at 0.7.

A sample profile from a magnetic feature with significant linear polarization is shown in Fig. 7. In the MHD cube, this pixel had a very complex stratification in B, in addition to gradients in γ and v_LOS. The asymmetry in Stokes V – which can be an observational indication of gradients in v_LOS – is present but very small. There are also small asymmetries in the σ-lobes of the linear polarization. Therefore, although scheme G does produce a better fit, and the approximation produced by SIR to the stratification of both v_LOS and γ is reasonable, it is also clear that scheme C produced a very good fit without the inclusion of additional nodes. A much larger number of nodes would be required to model the stratification of B. The T is very well constrained by both schemes. One can observe that scheme C and G naturally agree in their values of B, γ, and v_LOS at optical depths where the response functions peak, at optical depths between log τ_5000 Å = −1.0 and −0.5.

Fig. 7. Sample synthetic profile with the inverted profiles and associated atmospheres. Stokes I, Q, U, and V (grey, solid lines) are shown for the 15 648.52 Å line (left panels), normalised by the average continuum intensity, along with the inverted profiles from scheme C (blue, dotted) and G (red, dashed) lines. The dashed, grey line shows the rest wavelength. The retrieved atmospheric parameters (namely T, B, γ, and vLOS) are shown as a function of log τ5000 Å (right panels). The stratification of the atmospheric parameters in the MHD model is also shown (grey, solid line).

3.2. Inversions at GREGOR resolutions

At the resolution of the GREGOR telescope, we can assume that we are not able to resolve the small-scale magnetic structures we observe. For that reason, when inverting data observed at this resolution, we have to assume the resolution element (in the x − y plane) contains more than one model. We then quantify the contribution of each model by the magnetic filling factor, α. Additionally, we include v_mic and v_mac as free parameters to allow SIR to account for the spectral PSF. Table 4 describes the inversion that we employ for inverting the spatially and spectrally smeared synthetic cube. As in Milić et al. (2019), we observe in the degraded cube that the inferred velocities are lower and that downflows in IGLs are observed over a smaller area. In this section of the study, we focus our attention on the magnetic parameters – namely B and γ. We are particularly concerned with the influence of noise on the retrieval of accurate magnetic parameters, especially γ. We cannot compare the original MHD values to the inverted parameters in specific pixels after degradation because the spatial and spectral PSF describe only their impact on observable quantities (i.e. on the Stokes vector) and it is highly likely that the impact on atmospheric parameters is non-linear. Therefore, instead, our approach is to first invert the degraded profiles and use the inverted atmosphere as a baseline. We compare the results of inversions of both V1 and V3, representing lower and upper estimates of the spatial resolution of the GRIS-IFU, and then attempt to establish the influence stray light has on the on retrieved parameters by comparing V2 with V3. We then add noise to the polarized Stokes vectors and invert again. We then apply the same principle component analysis (PCA) procedure as Khomenko et al. (2003), Martínez González et al. (2016), Campbell et al. (2021) and invert once more. We can then compare the atmospheric models of the noiseless, noisy and de-noised inversions and examine the impact of S/N on our results. In this way, we are also critically examining the effectiveness and correctness of the PCA procedure.

Figure 8 shows the distributions of B, γ, and αB returned by the inversion under S1 for the noiseless and noisy V3 datacubes. On immediate inspection of the B distribution, compared to Fig. 4, we can observe the general impact of the smeared polarization profiles on inferred values; the strongest magnetic elements have been smeared over larger areas, resulting in a higher proportion of pixels with B values in excess of 100 G. The distribution of B and αB is very similar to that reported by Martínez González et al. (2016) and Campbell et al. (2021). Regarding the γ distribution, we know the vast majority of magnetic fields in the MHD cube have inclined (transverse) components, but as the amplitude of the relatively weak linear polarization signals are reduced by smearing, the degraded inverted atmosphere is much more vertical. In particular, the population of the highly inclined fields are impacted to a much larger extent than the intermediately inclined fields. Nevertheless, there remains a significant proportion of fields with transverse components.

Fig. 8. Distributions of B (left), γ (middle), and αB (right) returned from the inversions under S1 at GREGOR resolutions (V3). The inversions were run without the addition of noise (red), and then repeated with noise levels of 8 × 10−4Ic (blue) and 3 × 10−4Ic added (grey). The histograms of the inversions with noise are shown for pixels with maximum amplitude in at least one Stokes vector > 3σn. The lower row shows the difference between the values returned from inversions of the noiseless case and the values returned when noise is added. The percentages in all histograms are weighted with respect to the total number of pixels, i.e. including those with no polarization.

Figure 8 also shows the percentage differences between the γ, B, and αB values for inversions of the noiseless and noisy inversions at the two noise levels. Using the noiseless case as the baseline, we define the percentage differences between the B values, δB as,

$$\begin{array}{c}\hfill \delta B=\frac{B^{\mathrm{noiseless}}-B^{\mathrm{noise}}}{B^{\mathrm{noiseless}}}\end{array}$$(7)

where B^noiseless is the value returned by SIR under S1 with no noise degradation and B^noise is the value when a given level of noise is added, quantified by σ_n. We similarly define δγ and δαB. It is clear that noise introduces an error into each of these parameters, and that with a lower noise level (higher S/N) we access polarization in a larger fraction of the FOV. With a higher S/N, a larger proportion of fields with strengths < 200 G and magnetic flux densities < 25 G are accessible. The distribution of δB in particular is asymmetric, suggesting that the error induced by noise results in an over-estimation more often than an under-estimation. However, the distributions of δγ and δαB are more symmetric. Table 5 shows the median of the absolute values in δB, δγ, and δαB with both noise levels for V3 when only pixels with maximum amplitude in at least one Stokes vector > 3σ_n in the 15 648.52 Å line are included, while Table 6 shows the same for 4.5σ_n. It is clear that the inclusion of noise results in deviations from the noiseless case. It is also clear that the median in |δB|, |δγ|, and |δαB| is slightly reduced with a larger noise threshold, but this improvement is small. In terms of B, while the median is about 20% or less, the deviation can be as large as 50% or higher, which is not a negligible deviation. However, the deviation in αB is smaller. Similarly, a deviation in γ as large as 25−50% is not uncommon, although the median is less than 10%. This validates the decision to characterize the field inclinations in broad terms (i.e. not in terms of precise values, but instead in terms of whether they appear highly vertical, highly inclined, or intermediate, as in Campbell et al. 2021). At the same time, we found no evidence that PCA alters the values in any statistically significant way – the distributions for the inverted atmosphere after the application of PCA looked identical – showing that when a sensible number of eigenvectors is selected it can be safely applied to the Stokes vector without modification of the inferred atmospheres.

As the next step in our analysis, we inverted the degraded profiles again at the two noise levels but set any individual Stokes vector that did not have a maximum amplitude > 3σ_n before application of PCA to zero. Figure 9 shows the map of γ inverted in this way alongside the noiseless case for V1 and V3. The maps at σ_n = 8 × 10⁻⁴I_c look very similar to the regions of interest (ROIs) reported by Campbell et al. (2021) in terms of inclined fields flanked by opposite polarity vertical fields. The importance of adequate S/N is revealed by the map inverted with a noise level of σ_n = 3 × 10⁻⁴I_c, as our results suggest a much larger number of the inclined fields in particular should be measurable at this S/N compared to the map inverted with a noise level of σ_n = 8 × 10⁻⁴I_c. At the same time, by comparing maps at a given S/N between the figures, one can observe the importance of spatial resolution also; an increase in the spatial resolution of 0.23″ enables a much larger number of inclined fields to be revealed, and the magnetic structures are revealed in increasing complexity.

Fig. 9. Maps of inclinations returned from S1 inversions for the spatially and spectrally degraded synthetic profiles, with no noise (left) and two different noise levels (8 × 10−4Ic, middle, and 3 × 10−4Ic, right). Only Stokes vectors with a maximum amplitude > 3σn are retained after PCA. Black areas are representative of pixels which do not reach this threshold. The lower row shows the maps for V1 (FWHM = 0.53″) while the upper row shows the same for V3 (FWHM = 0.3″, with stray light included in the estimation and application of spatial PSF). The arrow in the upper left panel indicates the spatial location of the sample profile shown in Fig. 11.

While throughout this section of analysis we have been comparing V1 to V3 (as opposed to V1 to V2), and interpreting differences in terms of the difference in spatial resolution, it is also true that V3 additionally includes stray light as an additional variable. Figure 10 shows the distributions of inverted atmospheric parameters for V2 and V3, where the only difference is the absence or inclusion of stray light, respectively. Evidently, the only atmospheric parameter impacted by the unpolarized stray light component is the T. The T values for V2 have a wider spread. This is a predictable difference, as the inclusion of stray light reduces the continuum intensity contrast, which T is very sensitive to. As Milić et al. (2019) also found, v_LOS is invariant whether stray light is included or not. We additionally report that B and γ are also invariant.

Fig. 10. Impact of unpolarized stray light on the retrieval of the atmospheric parameters. Distributions of B, T, vLOS, and γ returned from the inversions under S1, for V2 (red lines) and V3 (blue lines) degraded synthetic profiles (i.e. without and with stray light, respectively). T is shown at log τ5000 Å = 0.5.

Figure 11 shows a sample profile from the degraded cube (V3), whose location is marked by the arrow in Fig. 9. This profile is symmetric enough that scheme S1 is able to provide a very good fit to the full Stokes vector. The profile and the inferred depth-averaged atmospheric parameters resemble closely the vectors and values recorded and inferred in linear polarization features by Campbell et al. (2021), in terms of being a highly inclined (γ = 97°) weak field (B = 337 G) with a αB value of 143 G. In principle, as both linear Stokes vectors are above the noise threshold, the azimuth should also be constrained (without resolution of the 180° ambiguity). The same profile in V1 shows a more highly inclined (γ = 100°) weaker field (B = 300 G) with a lower αB value of 106 G. We cannot directly compare these values to MHD values. However, it is worth noting that there is a mix of opposite polarities in this region in the undegraded cube, showing that the spatial and spectral PSFs fundamentally alter how we observe and interpret the magnetic vector when we are not able to resolve the small-scale magnetic structures. Further, the amplitudes of each polarized Stokes vector is 2−3 times smaller, highlighting the importance of adequate S/N in allowing us to constrain the magnetic vector in the IN at GREGOR resolutions, as well as the need for higher resolution telescopes.

Fig. 11. Sample V3 profile, showing the full Stokes vector for the 15648.52/15652.87 Å line pair in the left four panels, along with the S1 synthetic profiles derived from the SIR inversions. The horizontal (dot-dashed) lines show the 3σn noise thresholds for the polarized Stokes vectors, while the vertical (dashed) lines denote the rest wavelengths of each spectral line. On the right, the retrieved atmospheric parameters from SIR is shown in the table, while the temperature as a function of optical depth is shown in the lower right plot with the original FALC input model. The arrow in Fig. 9 indicates the spatial location of the sample profile.

4. Discussion

We have used a snapshot of a realistic MURaM simulation to show the extent to which we can accurately constrain information about the magnetic vector in the IN photosphere using inversion techniques. For inversions with constant parameters (scheme C, i.e. a Milne-Eddington atmosphere), we find the depth-averaged parameters are well constrained in congruence with previous studies (see, for instance, Orozco Suárez et al. 2010). The results with simple gradients (scheme G) in B, γ, and v_LOS show mixed results. The inclusion of a gradient in the B and v_LOS seems to improve the χ² and reduce the statistical error, suggesting that including gradients in v_LOS and B may be warranted. However, the inclusion of gradients in the γ in fact makes the statistical error slightly worse. The extent to which additional nodes can be added without the solutions becoming unconstrained or, even, oscillatory in optical depth is a question worthy of future investigation, as is the impact of the placement of nodes at specific optical depths (Socas-Navarro 2004; Socas-Navarro et al. 2008). We find high correlations between the inverted and MHD parameters, and the depths at which B, γ, v_LOS, and T are highly correlated resembles the response functions very closely. It is clear that both v_LOS and T are responsive to very deep layers of the atmosphere but we do not find evidence that B is responsive to the same layers. The weakest correlation is observed in γ and this parameter appears to be retrievable in the most narrow range of optical depths. Future studies that use SIR may be able to improve on these results if nodes can be placed at specifically selected optical depths, but ultimately more spectral lines will be required to probe the atmosphere at a greater range of optical depths. The recent publication of an expanded list of empirically determined atomic line parameters of spectral lines contained in the spectral region observed by GRIS (Trelles Arjona et al. 2021) is worthy of further investigation to examine whether the inclusion of any of the spectral lines could improve the inversions. Otherwise, future multi-instrument telescopes (e.g. DKIST) will provide an opportunity to sample a more varied selection of spectral lines with sensitivities to a larger range of optical depths.

We then investigated how this model atmosphere would appear had it been observed at the spatial and spectral resolutions achievable by the GREGOR/GRIS-IFU and compared the results to observations reported by Campbell et al. (2021). We found that the synthetic spectra closely resemble the real observations, with the caveat that a higher fraction of the FOV contains circular polarization signals when degraded than we observe in real observations. We studied a linear polarization feature (see Fig. 11) which resembles those reported by Campbell et al. (2021) in terms of being flanked by opposite polarity vertical fields, resembling ‘loop-like’ structures, but also in terms of its magnetic flux density (αB = 143 G), which is very close to the values estimated in similar features observed by the GREGOR/GRIS-IFU. This shows that MURaM simulations are capable of showing close agreement with real observations, and the symbiosis between realistic MHD simulations and observations continues to prove essential. The highly significant impact of spatial and spectral degradation on the amplitudes of polarized Stokes profiles, and consequently our ability to measure them with modern polarimeters, is clear. Instrumental effects can result an atmosphere with nearly one fifth of pixels containing linear polarization signals high enough to be observed, to consequently have a reduced amplitude such that the signals are measurable in less than a few percent of the FOV, depending on the S/N and tolerance for noise. At the same time, we found that B, γ, and v_LOS are invariant, and that the contrast in T is reduced, when an unpolarized stray light contribution is included.

As the linear polarization signals are typically so close to the noise level at GREGOR resolutions, the inclination and azimuth values returned from inversions are highly susceptible to noise contamination in the IN photosphere. As Campbell et al. (2021) found, this makes it very difficult to observe the temporal evolution of weak and inclined fields; their linear polarization signals can vanish between frames. The authors elected not only to restrict analysis to pixels which had just one polarization signal greater than a noise threshold, but also to set any individual Stokes vector with a maximum amplitude lower than this noise threshold to zero (as in Fig. 9). This was intended to restrict analysis to only those Stokes vectors that were actually measured, and is necessary as inversion codes like SIR will insert spurious Stokes vectors into pure photon noise. Consider the case where we have measured only Stokes V above a noise threshold. If we have not confidently measured at least Stokes Q or U, the only information we would have about the linear polarization signals is that their amplitude is not larger than the noise threshold. However, the amplitude of the Stokes Q or U signals could be any value between the noise threshold and zero. Depending on the relative amplitude of the corresponding Stokes V profile, this can result in large errors in γ (as in Fig. 8). This difficulty is even further encountered when trying to measure azimuths accurately, as one would need to measure both linear Stokes vectors confidently. Faced with this problem, we are left with a number of unpalatable options. We could:

1. Attempt to bin pixels spatially, spectrally, or temporally to reduce the noise level, as in Lagg et al. (2016), Campbell et al. (2021), at the risk of fundamentally altering the atmospheres due to Zeeman mixing of Stokes signals,

2. Set the Stokes vectors that are not measured above an appropriate confidence level to zero, and accept that we can only determine whether fields are predominantly vertical or horizontal in most pixels, or

3. Invert the profiles with noise, knowing that the inversion will produce an inflated number of transverse fields beyond what has actually been measured in the observations, that the retrieved magnetic parameters may have large errors induced by noise, and thus any temporal variation may be spurious.

If inversions are conducted with both the latter approaches, we can apply lower and upper limits, respectively, on the number of transverse fields. The errors introduced by noise, if similar in magnitude to those in Fig. 8 and Table 5, might be acceptable for retrieving statistical results. However, when trying to constrain the magnetic vector in small-scale features by following the dynamics in given pixels or groups of pixels between frames the errors could be much larger than the medians we report here. Therefore, it naturally follows that perhaps future observers should target a higher S/N. Future observers attempting to observe these magnetic fields with larger telescopes will also have to consider the implications of these results. Hopefully, with higher spatial resolution the measured polarization signals will have a higher amplitude, as there will be less mixing of opposite polarities. For this reason, next generation solar telescopes like DKIST and the European Solar Telescope (Matthews et al. 2016) will revolutionize our understanding of IN magnetic fields.

Acknowledgments

We are grateful to have benefited from the insight and expertise of the anonymous referee who helped improve this manuscript. R.J.C. additionally thanks all speakers, organizers and participants of the SOLARNET-funded Summer School held in September 2019 in the Universitä della Svizzera italiana, Lugano, Switzerland. In particular, Dr. Juan Manuel Borrero and Dr. Adur Pastor Yabar are thanked for useful discussions that helped guide the approach taken in this study. R.J.C. thanks Dr. Robert Ryans for IT support and assistance in utilizing QUB’s high performance computing (HPC) facilities. Dr. Ricardo Gafeira is thanked for helping us make use of his parallelized SIR wrapper, which is highly recommended. R.J.C. acknowledges support from the Northern Ireland Department for the Economy (DfE) for the award of a PhD studentship. This research has received financial support from the European Union’s Horizon 2020 research and innovation program under grant agreement No. 824135 (SOLARNET). M. M. acknowledges support from the Science and Technology Facilities Council (STFC) under grant No. ST/P000304/s & ST/T00021X/s. The 1.5-meter GREGOR solar telescope was built by a German consortium under the leadership of the Leibniz-Institute for Solar Physics (KIS) in Freiburg with the Leibniz Institute for Astrophysics Potsdam, the Institute for Astrophysics Göttingen, and the Max Planck Institute for Solar System Research in Göttingen as partners, and with contributions by the Instituto de Astrofísica de Canarias and the Astronomical Institute of the Academy of Sciences of the Czech Republic.

References

All Tables

All Figures

Fig. 1. Synthetic data produced from synthesis of MHD simulation outputs using SIR, shown at a continuum wavelength of 1565.47 nm. The top left panel shows the original Stokes I synthetic map (V0), while the top right and middle left panels show the spatially degraded map at an effective resolution of 380 km (V1) and 215 km (V2). The middle right panel shows the same as V2 but further degraded by stray light (V3). The bottom panel shows a histogram of the synthetic continuum intensities in each case, in addition to the distribution from real GREGOR/GRIS-IFU observations from 6 May 2019.

In the text

Fig. 2. Comparison of the synthesized (left) and inverted (middle) maps of the LP (upper row) and CP (lower row) according to scheme G. For both LP and CP the synthesized and inverted parameters are plotted against each other with the colour scale on a log10 scale (right). LP and CP, defined in Eqs. (4) and (5), respectively, are computed by integrating across the 15 648.52 Å line. The solid, green line indicates the exact linear relationship. The dashed, green ellipse indicates the 99% confidence level.

In the text

Fig. 3. Comparison of the true MHD (first column) and inverted (second column) maps of the T (top row), B (second row), γ (third row), and vLOS (bottom row). The maps of inverted parameters are shown as determined by inversion scheme G (except for γ, which is shown as determined by scheme C). The colour table for B is extended higher than the maximum value (1915 G) to accommodate those with colour vision deficiency. The MHD and inverted parameters, for both scheme C (third column) and G (fourth column) are plotted against each other with the colour scale on a log10 scale. The solid, green line indicates the exact linear relationship. The dashed, green ellipse indicates the 99% confidence level. T is shown at log τ5000 Å = 0.5, while the other parameters are shown at log τ5000 Å = −0.5.

In the text

Fig. 4. Distributions of values (left) and histograms between the ‘true’ MHD values and the differences between ‘true’ MHD and inverted values (middle and right, for scheme C and G, respectively) in B, γ, vLOS, and T (upper to lower row, respectively). Values are retrieved at log τ5000 Å = −0.5, except for T which is retrieved at log τ5000 Å = 0.5. In the distributions, values from the MHD cube (orange lines) and returned from the inversions under scheme C (red lines) and G (blue lines) are shown. The histograms of B, γ, ΔB, and Δγ are shown for pixels with maximum amplitude in at least one polarized Stokes vector > 2.4 × 10−3Ic. The median error as a function of each parameter is shown by the green, dashed lines. The percentages in all histograms are weighted with respect to the total number of pixels, N, i.e. including those with no polarization. The colour scale is shown on a log10 scale.

In the text

	Fig. 5. Map of the normalized difference between χ2 values produces by scheme C and scheme G inversions.
In the text

	Fig. 6. Pearson correlation coefficients for values of B, γ, vLOS, and T between the MHD cube and scheme G estimates at a range of optical depths. The colour scale saturates at 0.7.
In the text

Fig. 7. Sample synthetic profile with the inverted profiles and associated atmospheres. Stokes I, Q, U, and V (grey, solid lines) are shown for the 15 648.52 Å line (left panels), normalised by the average continuum intensity, along with the inverted profiles from scheme C (blue, dotted) and G (red, dashed) lines. The dashed, grey line shows the rest wavelength. The retrieved atmospheric parameters (namely T, B, γ, and vLOS) are shown as a function of log τ5000 Å (right panels). The stratification of the atmospheric parameters in the MHD model is also shown (grey, solid line).

In the text

Fig. 8. Distributions of B (left), γ (middle), and αB (right) returned from the inversions under S1 at GREGOR resolutions (V3). The inversions were run without the addition of noise (red), and then repeated with noise levels of 8 × 10−4Ic (blue) and 3 × 10−4Ic added (grey). The histograms of the inversions with noise are shown for pixels with maximum amplitude in at least one Stokes vector > 3σn. The lower row shows the difference between the values returned from inversions of the noiseless case and the values returned when noise is added. The percentages in all histograms are weighted with respect to the total number of pixels, i.e. including those with no polarization.

In the text

Fig. 9. Maps of inclinations returned from S1 inversions for the spatially and spectrally degraded synthetic profiles, with no noise (left) and two different noise levels (8 × 10−4Ic, middle, and 3 × 10−4Ic, right). Only Stokes vectors with a maximum amplitude > 3σn are retained after PCA. Black areas are representative of pixels which do not reach this threshold. The lower row shows the maps for V1 (FWHM = 0.53″) while the upper row shows the same for V3 (FWHM = 0.3″, with stray light included in the estimation and application of spatial PSF). The arrow in the upper left panel indicates the spatial location of the sample profile shown in Fig. 11.

In the text

	Fig. 10. Impact of unpolarized stray light on the retrieval of the atmospheric parameters. Distributions of B, T, vLOS, and γ returned from the inversions under S1, for V2 (red lines) and V3 (blue lines) degraded synthetic profiles (i.e. without and with stray light, respectively). T is shown at log τ5000 Å = 0.5.
In the text

Fig. 11. Sample V3 profile, showing the full Stokes vector for the 15648.52/15652.87 Å line pair in the left four panels, along with the S1 synthetic profiles derived from the SIR inversions. The horizontal (dot-dashed) lines show the 3σn noise thresholds for the polarized Stokes vectors, while the vertical (dashed) lines denote the rest wavelengths of each spectral line. On the right, the retrieved atmospheric parameters from SIR is shown in the table, while the temperature as a function of optical depth is shown in the lower right plot with the original FALC input model. The arrow in Fig. 9 indicates the spatial location of the sample profile.

In the text