© The Authors 2023

1. Introduction

Electric currents, j, in the solar atmosphere play a very important role. They are crucial for the energy balance through Joule and ambipolar heating (Priest 1999; Cheung & Cameron 2012; Khomenko et al. 2021). They can also be employed as proxies of regions where the magnetic field is highly non-potential and therefore likely to cause transient and explosive events in the solar chromosphere and corona (Priest & Forbes 2002; Green et al. 2018). Consequently, measuring electric currents has been long considered one of the most important goals in solar physics (Jurčák et al. 2019). The inference of electric currents is possible via magnetic field extrapolations (Louis et al. 2021). However, the extrapolated magnetic field differs both in magnitude and fine structure from the actual measurements (Vissers et al. 2022), thereby casting doubts on the electric currents inferred in this fashion.

It is therefore more compelling to use the measured magnetic field instead. Most previous determinations of electric currents based on the measured B were carried out through Stokes inversion techniques that employ the Milne-Eddington (ME) approximation (Solanki et al. 2003; Wang et al. 2017). However, ME inversion techniques have been used to calculate only the vertical component of the electric current density, j_z, instead of the full vector, j. This occurs because ME inversion yields constant values along the z direction (Auer et al. 1977; Landolfi & Landi Degl’Innocenti 1982) and therefore lacks the information needed to fully evaluate j_x and j_y. Non-ME Stokes inversion codes are indeed capable of inferring j, but they do so by calculating a z scale that is based on hydrostatic equilibrium (HE; Socas-Navarro 2005). To our knowledge, the only work in which the three components of the electric current density, j, were determined from the measured magnetic field B, not using HE, is Puschmann et al. (2010b). This work is instead based on the assumption of magnetohydrostatic (MHS) equilibrium (Puschmann et al. 2010a).

Despite being of great importance, the study by Puschmann et al. (2010a,b) suffers from a number of shortcomings (Borrero et al. 2019). In Borrero et al. (2021) we present a novel Stokes inversion code that also employs MHS equilibrium and solves most of those shortcomings. The resulting method was subsequently employed to determine electric currents in the solar photosphere and compared to results from magnetohydrodynamic (MHD) simulations (Pastor Yabar et al. 2021), thereby allowing us to study, for the first time, the accuracy in the inference of j. However, a missing piece of the puzzle is whether this accuracy actually improves results that employ MHS instead of HE and, if so, by how much. In addition, we study the effects of, among other factors, spectral resolution, photon noise, and the number of spectral lines inverted. To this end, in this work we again employ MHD simulations to calculate simulated Stokes profiles in a number of spectral lines (Sect. 2). These Stokes profiles are then subjected to a number of Stokes inversion techniques with increasing degrees of complexity (including HE, MHS equilibrium, etc.; Sect. 3). Results are discussed and compared in Sect. 4. The effects of photon noise and limited spectral resolution are studied in Sect. 5. Finally, Sect. 6 summarizes our findings.

2. MHD simulations and synthetic Stokes profiles

In this paper we employ the same physical parameters as in Pastor Yabar et al. (2021). These parameters come from one single snapshot that resulted from the three-dimensional radiative MHD simulations by Rempel (2012). The original simulation box has a size of 4096×4096×768 cells, with the first two values corresponding to the solar surface, XY, and the third one, Z, perpendicular to the surface (i.e., antiparallel to the Sun’s gravity direction). The grid spacing is Δx × Δy × Δz = 12 × 12 × 8 km in each spatial direction. The simulation setup is that of a large sunspot surrounded by granulation and different magnetic knots. From here we selected a smaller region of 512×512 cells on the XY plane. This region includes a couple of pores embedded in granulation (see Appendix B in Pastor Yabar et al. 2021). The magnetic field from the simulations, B_mhd, was employed to determine the electric current density by applying Ampere’s law: $\mathit{j}={\mu }_0^{-1}\mathrm{\nabla }\times \mathit{B}$, where j is given in A m⁻² when B is given in teslas, and μ₀ = 4π10⁻⁷ kg m s⁻² A⁻² is the vacuum magnetic permeability. Figure 1 illustrates the logarithm of the absolute value of the three components of electric current density, j_mhd: j_x, mhd (first column), j_y, mhd (second column), and j_z, mhd (third column), as well as its modulus, $\Vert {\mathit{j}}_{\mathrm{mhd}}\Vert =\sqrt{j_x^2+j_{\mathit{y}}^2+j_z^2}$ (rightmost column), in SI units of amperes per m².

Fig. 1. Logarithm of the absolute value of the electric current density in the MHD simulations. From left to right: jx, mhd(x, y, τc), jy, mhd(x, y, τc), jz, mhd(x, y, τc), and ∥jmhd(x, y, τc)∥. The electric currents are shown at two different optical depths: τc = 0.1 (upper panels) and τc = 1 (lower panels).

The physical parameters from the MHD simulations were employed to solve the radiative transfer equation for polarized light, using the forward module of the FIRTEZ code (Pastor Yabar et al. 2019) so as to obtain the Stokes vector as a function of wavelength, 𝕀^syn(x, y, λ), in a number of spectral lines. We refer to these as “synthetic” observations. Although FIRTEZ’s solver for the radiative transfer equation considers the physical parameters to be constant within each grid cell, in reality they are expected to vary internally. The effect of sub-grid variations (i.e., at scales smaller than Δz <  8 km) is not investigated in this paper.

The selected spectral lines, number of wavelengths, spectral sampling¹, and atomic parameters are given in Table 1. Unless otherwise specified, normally distributed random noise – with a σ_pn = 10⁻³ (in units of the continuum intensity), where pn refers to photon noise – was added to 𝕀^syn(x, y, λ) to mimic the effects of photon noise. The synthetic observations are not spatially convolved with a point spread function, so the equivalent horizontal spatial resolution of the observations is that of the original MHD simulation: Δx = Δy = 12 km. This is similar to the resolution that new-generation solar telescopes will achieve (Collados et al. 2013; Elmore et al. 2014; Jurčák et al. 2019; Rimmele et al. 2020).

3. Inference of electric currents via Stokes inversion

In this section we perform a series of inversions of the synthetic observations, 𝕀^syn(λ), described in Sect. 2 to recover the physical parameters of the solar atmosphere: T_inv(x, y, z), ρ_inv(x, y, z), P_g, inv(x, y, z), B_inv(x, y, z), and v_z, inv(x, y, z). The inversions were all performed using the FIRTEZ-dz code (Pastor Yabar et al. 2019). We note that some other well-known Stokes inversion codes are able to carry out some of the inversions described in this section. We will mention which codes can be used and when. After the inversion, the horizontal components of the magnetic field (B_x, B_y) were corrected to remove the well-known 180° ambiguity by applying the non-potential magnetic field calculation (Georgoulis 2005). After this, Ampere’s law was used to determine the electric current density, j_inv, which was then compared to the currents from the numerical simulations, j_mhd (see Fig. 1). In the following sections (Sects. 3.1–3.6), Stokes inversions are carried out with increasing degrees of complexity to study if and how the inferred electric current density improves.

3.1. Milne-Eddington-like inversion (me)

As a first step we performed an inversion of only the first spectral range in Table 1: the Fe I line pair at 630.2 nm that was typically observed by the Advanced Stokes Polarimeter (ASP; Elmore et al. 1992) and later became the lines of choice of the spectropolarimeter on board Hinode (Lites et al. 2001; Ichimoto et al. 2007).

The inversion here was carried out along the z direction independently for each (x, y) column with the following numbers of nodes: 32 for T(z), 16 for v_z(z), 1 for B_x(z), 1 for B_y(z), and 1 for B_z(z). Because only one node is given to the each of the three components of the magnetic field, FIRTEZ-dz retrieves values that are constant in z, and therefore these results are comparable to a traditional ME inversion, such as those carried out with the ASP/HAO code (Lites & Skumanich 1990; Lites et al. 1993), HELI X⁺ Lagg et al. (2004, 2009), or VFISV (Borrero et al. 2011). A summary and comparison of these three codes can be found in Borrero et al. (2014). However, unlike with ME Stokes inversion codes, we still allow for a z-dependent v_z and T_z. The temperature in particular is important as it means that we solve the Saha and Boltzmann equations (i.e., we assume local thermodynamic equilibrium) instead of parameterizing the source function with two free parameters, as done by ME inversion codes. In order to do so, the gas pressure and density must also be known. This was accomplished by assuming HE.

During the inversion, different weights were given to each of the four Stokes parameters: w_i = 1, w_q = w_u = 4, and w_v = 2. We refer to this inversion as me.

3.2. Inversion under hydrostatic equilibrium (he1)

The next step was to drop the assumption made by ME codes and allow for the three components of the magnetic field to vary with height. In this case, the number of nodes was increased to: four for B_x(z), four for B_y(z), and eight for B_z(z). The T(z) and v_z were again given 32 and 16 nodes, respectively. All other inversion parameters (weights, HE, etc.) were the same as in Sect. 3.1. We refer to the inversion described in this section as he1. Examples of Stokes inversion codes that are capable of performing this kind of analysis are: SIR (Ruiz Cobo & del Toro Iniesta 1992), SPINOR (Frutiger et al. 1999; van Noort 2012), NICOLE (Socas-Navarro et al. 2015), and SNAPI (Milić & van Noort 2018).

3.3. Increased number of spectral windows (he2)

In the he2 case, the Stokes inversion was carried out in the same fashion as for he1 (HE; Sect. 3.2) but with a greater number of spectral regions, and thus a greater number of spectral lines, such that all lines in Table 1 were included in the inversion. The idea behind this is twofold: (a) the addition of redundant information about the magnetic field in the region where the Fe I line pair at 630.2 nm is formed should improve the determination of the derivatives of the magnetic field components with respect to the vertical axis, z, in this region, and (b) the addition of spectral lines that convey information about the magnetic field outside the region where the Fe I line pair at 630.2 nm is formed should allow us to determine the electric current density over a wider range of heights (or optical depths). All inversion codes mentioned in Sect. 3.2 (he1) can also perform the inversions described here.

3.4. Increased number of nodes (he3)

In the he3 case, the Stokes inversion was carried out in the same fashion as in he2 (HE plus multiple spectral ranges or lines; Sect. 3.3) but now with the number of nodes increased to 16 for B_x(z), 16 for B_y(z), and 32 for B_z(z). The idea behind this inversion is that in regions where the kinetic energy is much higher than the magnetic energy, the MHD simulations feature very large magnetic field variations along the z direction. This occurs because, in this case, the convective motions in the solar atmosphere braid the magnetic field lines into non-potential configurations where the magnetic field can present very large variations over small spatial scales. These variations lead, in general, to very asymmetric Stokes parameters (Landolfi & Landi Degl’Innocenti 1996; Borrero et al. 2007) that can only be reproduced by inversion if enough free parameters are given to the three components of the magnetic field. Examples of such large z variations, and how Stokes inversion techniques perform with a low number of nodes, can be found in the literature (see, e.g., Fig. 2 in Bellot Rubio 2006 and Fig. 5 in Borrero et al. 2014). Of course, it is always a wise policy to increase the number of free parameters only if enough information is contained in the synthesized spectral lines, and therefore this new degree of complexity was only added in the he2 case and not he1 because the former considers more spectral lines. Again, all inversion codes mentioned in Sect. 3.2 (he1) can also perform the inversions described here.

3.5. Effect of vertical Tikhonov regularization (he4)

In the he4 case, the Stokes inversion was carried out in the same way as in he3 (HE, plus multiple spectral ranges or lines, plus an increased number of nodes; Sect. 3.4), but with a Tikhonov-like regularization included along the z direction. This regularization penalizes solutions where the components of the magnetic field present very large variations along the z direction, unless those variations are truly needed in order to properly fit the synthetic Stokes profiles, 𝕀^syn(x, y, λ). This was done via a modification of the χ² merit function that is being minimized by the Stokes inversion code (de la Cruz Rodríguez et al. 2019). One might wonder whether limiting the possible vertical variations in the three components of the magnetic field works against the previous Stokes inversion (he3), where we allowed for larger variations along the z-axis, and whether, because of this, he4 is comparable to he2. This is in principle a valid concern; however, we need to take into consideration the fact that inversion he2 limited the vertical variations everywhere, whereas he3 allowed for larger vertical variations everywhere. As we pointed out in Sect. 3.4, the large numbers of free parameters might be needed only in regions where the kinetic energy is higher than the magnetic energy. Therefore, in regions where the opposite occurs (i.e., umbra, pores, etc.), this might not be necessary and could even be counterproductive. To this end, the regularization included in he4 allows the inversion code FIRTEZ-dz to choose between either of the aforementioned situations based on whether the fit to the synthetic Stokes profiles truly improves when large z variations are included; otherwise, it favors smoother solutions along the z direction. Other codes capable of carrying out Stokes inversion with vertical regularization are STiC (de la Cruz Rodríguez et al. 2019) and TIC (Li et al. 2022).

3.6. Inversion under magnetohydrostatic equilibrium (mhs)

In this new inversion, which we call mhs, we included all previous improvements (multiple spectral lines or ranges), increased number of nodes, and vertical z Tikhonov regularization), and, in addition, we dropped the assumption of HE and instead used MHS equilibrium. This was done following the iterative approach described in Borrero et al. (2021, Fig. 2) but including all terms of the divergence of the Lorentz force (see Eq. (B.7) in that paper). This is the most complex inversion tested in this paper. Only the FIRTEZ-dz Stokes inversion code is capable of including MHS constraints.

4. Discussion

4.1. Summary of results

Let us first briefly consider Fig. 2. This figure shows the same plot as Fig. 1 but as inferred from the me inversion (Sect. 3.1). Visual comparison between these two figures reveals that the Stokes inversion is in general capable of grasping the overall spatial distribution of the electric current density, although the magnitude and horizontal thickness of the current sheets are clearly underestimated by the me inversion. We note that Fig. 2 also includes j_x and j_y despite the ME-like inversion retrieving ∂B_y/∂z = ∂B_x/∂z = 0 (see Eq. (1)). A more detailed comparison can be made from Fig. 3. This new figure presents, from left to right, density plots of the logarithm of the absolute value of j_x, j_y, j_z, and ∥j∥ in the MHD simulations (horizontal axis) versus those inferred from the me inversion (vertical axis). For better visualization of low (blue/purple) and high (yellow) density regions, we employ a logarithmic color scale. The more the yellow regions are concentrated along the diagonal black line, the better the inference of the electric currents is. These plots are in qualitative agreement with those presented in Pastor Yabar et al. (2021) in that the larger the value of the electric current density, the better the inference, in particular for values ⪆0.1 A m⁻². This feature is seen not only in me, but in all tested Stokes inversions.

Fig. 2. Same as Fig. 1 but showing the electric currents as obtained from the ME-like inversions described in Sect. 3.1 (he1).

Fig. 3. Logarithmic scatter-density plots of the electric current density from the MHD simulations (horizontal axis) and the values inferred by the me Stokes inversion (vertical axis) at τc = 1 (lower panels) and τc = 0.1 (upper panels). From left to right: logarithm of the absolute value of the ∥jx∥ component of the electric current density, logarithm of the absolute value of the ∥jy∥ component of the electric current density, logarithm of the absolute value of the ∥jz∥ component of the electric current density, and logarithm of the absolute value of the modulus of the electric current density, ∥j∥. The horizontal dashed black line represents the so-called floor value, below which the Stokes inversion is not capable of correctly inferring the electric currents (see text for details). The solid thick diagonal corresponds to the one-to-one retrieval along which the Stokes inversion perfectly retrieves the electric currents. All pixels in the selected region described in Sect. 2 are included in this figure.

In addition, it is worth noting that, on the low end of the electric current density, the inversion yields an almost constant value of ∥j_inv∥≈10^−1.3 ≈ 0.05 A m⁻² (see the horizontal thick dashed line in the rightmost panel of Fig. 3), even if the electric current in the original MHD simulation is much lower than this value. In Pastor Yabar et al. (2021) we referred to this value as the “floor” value, and it corresponds to the smallest value of electric current density that can be inferred by the Stokes inversion. Although this floor is visually clearest in the modulus of the electric current density vector, ∥j∥ (rightmost panel in Fig. 3), it can been seen in all three individual components of the electric current vector as well. The floor values for all Stokes inversions carried out in this work are presented in Table 2. In this table we provide f, which we define as the power of ten that yields the floor value as 10^−f. With this definition, the larger the value of f, the better a given inversion type is at inferring smaller values of the electric current density. The numbers in this table were obtained by calculating the first moment of the distribution of log∥j_inv∥ in the regions where ∥j_mhd∥< 0.01 A m⁻² as well as for the three individual components of the electric current density vector. According to Table 2, f is smaller for ∥j_z∥ (leftmost panel) than for ∥j_x∥ and ∥j_y∥, and therefore the latter components are less affected by the floor value than the former. This also makes the vertical component of the electric current density, j_z, the main contributor to the floor value of the modulus ∥j∥. The floor values are also slightly smaller at τ_c = 1 than at τ_c = 0.1, meaning that it is possible to detect slightly lower electric currents at the latter optical depth (τ_c = 0.1) than at the former (τ_c = 1). Floor values are, however, of limited relevance for our work because such small electric current density values play a negligible role in the energy balance of the solar photosphere, as Ohmic heating varies with ∥j∥².

A more quantitative analysis can be made via Fig. 4, where we plot histograms of the logarithm of the absolute value of the ratio between the currents present in the MHD simulations and the currents inferred from the Stokes inversion: j_x (top left; panel a), j_y (top right; panel b), j_z (bottom left; panel c), and ∥j∥ (bottom right; panel d). Figure 4 corresponds to τ_c = 1. Very similar results are obtained for τ_c = 0.1 (not shown). Results for τ_c <  10⁻² deteriorate significantly (see Sect. 4.3).

Fig. 4. Histograms of the logarithm of the absolute value of the quotient of the components of the electric current density from the MHD simulations and as obtained from the Stokes inversion: (a) ∥jx∥ (upper left); (b) ∥jy∥ (upper right); (c) ∥jz∥ (lower left); and (d) ∥j∥ (lower right). The histograms correspond to the Stokes inversions described in Sect. 3: red for me, yellow for he1, green for he2, cyan for he3, blue for he4, and purple for mhs.

Pixels on the positive side of the horizontal axis in Fig. 4 correspond to regions where the inversion underestimates the true values from the MHD simulations, whereas on the negative side of the horizontal axis the inversion overestimates the values from the MHD simulations. All histograms in this figure, in particular those corresponding to j_z and ‖j‖ (panels c and d), feature a negative skewness, with a long tail on the negative side of the horizontal axis. This tail is produced by the floor value described above, where the inversion overestimates the values of the electric current density. In addition, it can be noted that the peak of the different histograms is usually shifted toward positive values, meaning that the inversion tends to underestimate the values of the current density in the MHD simulations. This underestimation happens for values of log(∥j_inv∥) >  0 A m⁻² (see the ∥j∥ retrieval in the rightmost panels in Fig. 3). This bias becomes less important as the complexity of the inversion setup increases, as demonstrated by the fact that the peak of the histograms in he4 and mhs are mostly centered at zero (see Fig. 4).

In order to evaluate the ability of the different inversions to retrieve the correct values of electric currents in the MHD simulations, we considered the standard deviation, σ, of the histograms². The smaller the value of σ, the better the retrieval of the corresponding component of the electric current density. Values of σ at two different optical depths, τ_c = 0.1, 1, in all considered inversions are presented in Table 3. For each inversion we provide two values of σ as XXX/YYY. The first value, XXX, corresponds to the standard deviation of the histogram in all pixels of the simulation (Fig. 4), whereas the second value, YYY, is the standard deviation obtained from the histograms using only those pixels where the magnetic field is stronger than 300 gauss (Fig. 5). This second estimation of σ is provided because, as pointed out in Pastor Yabar et al. (2021), the electric current density is determined much more reliably in regions that harbor stronger magnetic fields. We emphasize here that the values given for the standard deviation of the histograms in Table 3 are given for log[∥j_mhd/j_inv∥]. This means that a standard deviation of, for example, 0.7 dex implies that the values of the electric current inferred through the inversion are, about 68% of the time, between j_inv ∈ [10^−0.7, 10^0.7]≈[1/5, 5]j_mhd. In other words, it implies that the electric current can be inferred within a factor of 5 of the original values (higher or lower) in the MHD simulations.

Fig. 5. Same as Fig. 4 but employing only regions where ∥B(τc = 1)∥ >  300 gauss.

An analysis of Table 3 reveals that the values of σ typically decrease from me to mhs. This implies that the inference of all three components of the electric current density improves as the complexity of the performed Stokes inversion increases. The overall improvement in σ from the simplest to the most complex inversion is about 25–50%. This is not due to a particular inversion setup; they each produce a similar and steady enhancement in the accuracy of the retrieval of the components of the electric current density. The only notable exception seems to be he4, where, as anticipated in Sect. 3.5, inferences of j_z at τ_c = 0.1 worsen with respect to he3. Table 3 also confirms that the retrieval of j_z is always less accurate than for j_x and j_y and that currents can be determined more reliably at τ_c = 1 than at τ_c = 0.1, although at τ_c = 0.1 it is possible to detect lower values of the components of the electric current density (see Table 2).

If we include all pixels in the simulation (using only the first number in Fig. 4), we can summarize our results by saying that, under the simplest Stokes inversion described in Sect. 3 (me), the three components of the electric current density can be determined with a standard deviation of σ = 0.90 − 1.00 dex. In the case of the most advanced Stokes inversion (mhs), the three components of the electric currents can be determined with an accuracy of σ = 0.70 − 0.85 dex. The total electric current density can be retrieved with σ = 0.75 dex for me and σ = 0.50 dex for mhs.

Table 3 also shows that the inference of the electric current density clearly improves when only regions of strong magnetic fields are considered (∥B(τ_c = 1)∥ >  300 gauss). This can be readily seen by comparing the widths and heights of the histograms in Fig. 5 with those in Fig. 4. Results in this case can be summarized as follows. Under me, the three components of the electric current density can be determined with a standard deviation of σ = 0.85 − 0.90 dex, whereas for mhs σ = 0.5 − 0.6 dex. When referring to the total current density, the accuracy ranges between σ = 0.65 dex for he1 and σ = 0.30 dex for mhs.

A visual impression of the overall improvement between me and mhs can be drawn by comparing Fig. 3 with Fig. 6, where we present the logarithmic scatter-density plots of the electric current density from the MHD simulations (horizontal axis) and the values inferred by the mhs Stokes inversion (vertical axis) at τ_c = 1 (lower panels) and τ_c = 0.1 (upper panels). These two figures illustrate the improvement between the simplest and the most complex inversion that we carried out. However, we emphasize that the improvement is not solely due to the use of MHS equilibrium (as opposed to HE) but also the cumulative effect of all the inversions described in Sect. 3.

Fig. 6. Same as Fig. 3 but showing the electric currents as obtained from the inversions under MHS equilibrium described in Sect. 3.6 (mhs).

4.2. Spatial derivatives

In Sect. 4.1 we have seen that j_x and j_y are always better retrieved than j_z. In order to understand the reason, we need to investigate our ability to infer the individual spatial derivatives that enter the calculation of the three components of the electric current density:

$$\begin{array}{cc}& j_x\propto (\partial B_y/\partial z)-(\partial B_z/\partial y)\hfill \\ \hfill & j_y\propto (\partial B_x/\partial z)-(\partial B_z/\partial x)\hfill \\ \hfill & j_z\propto (\partial B_x/\partial y)-(\partial B_y/\partial x).\hfill \end{array}$$(1)

Figure 7 shows histograms of the logarithm of the ratio between the derivatives in the MHD simulations and the derivatives as inferred from the different Stokes inversions described in Sect. 3. The histograms are similar to those in Fig. 4 and so are their interpretation: the positive region on the horizontal axis corresponds to pixels where the spatial derivative in the MHD simulation is larger than that inferred through the inversion, whereas the negative region on the horizontal axis corresponds to pixels where the spatial derivative inferred through the inversion is larger than the derivative in the MHD simulations. In this figure we show only three of the six derivatives in Eq. (1) because, statistically speaking, ∂B_z/∂y (top panel) behaves identically to ∂B_z/∂x, ∂B_x/∂y (middle panel) behaves identically to ∂B_y/∂x, and finally ∂B_y/∂z (bottom panel) behaves identically to ∂B_x/∂z. Therefore, these three panels suffice for the argumentation that follows.

Fig. 7. Histograms of the logarithm of the absolute value of the quotient of the spatial derivatives of the components of the electric current density from the MHD simulations and as obtained from the Stokes inversion. Top panel: ∂Bz/∂y (or alternatively ∂Bz/∂x). Middle panel: ∂Bx/∂y (or alternatively ∂By/∂x). Bottom panel: ∂By/∂z (or alternatively ∂Bx/∂z). The histograms correspond to the Stokes inversions described in Sect. 3: me (red), he1 (yellow), he2 (green), he3 (cyan), he4 (blue), and mhs (purple).

The first obvious conclusion we can draw from Fig. 7 is that, just as we saw with the electric current density, the inference of the spatial derivatives of the magnetic field steadily improves as the complexity of the inversion increases. Histograms are narrower and more centered at the zero value as we go from me (Sect. 3.1) to mhs (Sect. 3.6). The second conclusion is that, from best to worst, the order in which spatial derivatives are best inferred is: ∂B_z/∂x (or ∂B_z/∂y; top panel in Fig. 7), ∂B_x/∂y (or ∂B_y/∂x; middle panel), and finally ∂B_y/∂z (or ∂B_x/∂z; bottom panel). The reason as to why this is the case can be understood in terms of the sensitivity of the three components of the magnetic field to the different Stokes profiles. The B_z is generally better retrieved than either B_x or B_y because signals in Stokes V are typically larger than in Q and U even if the vertical and horizontal components of the magnetic field have the same value (see, e.g., Fig. 5 in Borrero & Kobel 2013). This explains why ∂B_z/∂x is better determined than both ∂B_x/∂y and ∂B_x/∂z. Now, ∂B_x/∂y is more reliably inferred than ∂B_x/∂z because for the former we employ data (i.e., the Stokes vector) from two different adjacent pixels, whereas in the latter we employ the data from just one. The increased amount of information available allows the horizontal derivative of the horizontal components of the magnetic field to be determined more accurately than the vertical derivative of the horizontal components of the magnetic field.

We switch our attention now to Eq. (1). Here we see that j_x and j_y are obtained from a combination of two spatial derivatives, one of which can be determined very accurately (∂B_z/∂y or ∂B_z/∂x) while the other is rather poorly determined (∂B_y/∂z or ∂B_x/∂z) and is actually not retrieved at all in the ME-like case (Sect. 3.1). Meanwhile, j_z is obtained from the combination of two spatial derivatives (∂B_x/∂y and ∂B_y/∂x), both of which can be determined with an accuracy that is somewhere in between the derivatives needed to determine j_x and j_y. This raises the question as to what makes the inference of j_x and j_y more reliable than j_z. We believe there are two main reasons. To begin with, we find that in the original MHD simulations (Sect. 2), at τ_c = 1 and τ_c = 0.1, ∥∂B_z/∂y∥> ∥∂B_y/∂z∥ in 55% of the analyzed horizontal area. Therefore, for the correct retrieval of j_x, it is slightly more important to accurately determine the first of these two spatial derivatives, which we actually do. An identical argument can be made for the case of j_y and its two spatial derivatives (Eq. (1)). Interestingly, this discrepancy does not appear between ∥∂B_x/∂y∥ and ∥∂B_y/∂x∥: neither one dominates over the other, and therefore j_z is not affected.

In addition, it can be seen that the inversion tends to overestimate ∥∂B_z/∂y∥, as indicated by the long tail toward negative values in the top panel of Fig. 7. At the same time, the inversion tends to underestimate ∥∂B_y/∂z∥, as demonstrated by the fact that the peak of the histogram in the bottom panel of Fig. 7 always lies on the positive side. Because of this, the overestimation of one derivative involved in the calculation of j_x (and of j_y) partially compensates for the underestimation of the other derivative, such that a final, accurate value is yielded.

4.3. Height dependence

In Sect. 4.1 we pointed out that the inference of electric currents quickly worsens for optical depths τ_c <  10⁻². We can offer two different explanations as to why this happens. The first is that for τ_c <  10⁻² (above the mid-photosphere) the determination of the spatial derivatives is worse than for τ_c ∈ [0.1, 1]. This is illustrated in Fig. 8, which is to be compared with Fig. 7. Unreliable results for the mid-photosphere and above might be caused by the spectral ranges selected for our analysis (Table 1). Our original intention was to include lines with sensitivity to the magnetic field in a wide variety of heights and optical depths (see Fig. 10 in Borrero et al. 2016). In particular, this was the rationale for including Si I 1082.7 nm. This was guided by the fact that, for the temperature T(z) and line-of-sight velocity v_z(z), this spectral line features a strong sensitivity to the mid- and upper-photosphere (Bard & Carlsson 2008; González Manrique et al. 2020). Unfortunately, this does not seem to be the case for the magnetic field, as the sensitivity peaks at τ_c = 10⁻² and quickly drops for smaller optical depths (Kuckein et al. 2012; Felipe et al. 2016; Shchukina & Trujillo Bueno 2019).

Fig. 8. Same as Fig. 7 but at τc = 10−2.

The second reason is that, even if the spatial derivatives of the magnetic field could be determined for τ_c = 10⁻² (Fig. 8) as accurately as for τ_c = 1 (Fig. 7), the determination of the components of the electric current density vector would still worsen because, unlike at τ_c = 1 (where ∥∂B_z/∂y∥> ∥∂B_y/∂z∥ in 55% of the analyzed area), at τ_c = 10⁻² the opposite happens: ∥∂B_z/∂y∥< ∥∂B_y/∂z∥ in about 65% of the area. And consequently, since ∥∂B_y/∂z∥ is less accurate than ∥∂B_z/∂y∥, the resulting j_x would also be less reliable.

The tendency for ∥∂B_z/∂y∥< ∥∂B_y/∂z∥ also continues at τ_c = 10⁻³, but higher up, at τ_c = 10⁻⁴, neither dominates over the other. Unfortunately, at this height the sensitivity to the magnetic field provided by the selected spectral ranges in Table 1 is already negligible. Moreover, at this height it is not even clear if the MHS approach employed here to determine the gas pressure and the z scale (Sect. 3.6) is valid. The reason is that at τ_c <  10⁻⁴, plasma velocities are comparable to the sound speed, and therefore the role of the advection term in the momentum equation cannot be neglected. In this case it would be more appropriate to employ a magneto-hydro-stationary approach instead, which unfortunately has not yet been developed.

It is worth mentioning that the fact that some derivatives dominate over others at different heights is a property of the MHD simulations and of the selected region within the entire simulation domain. We believe this region is large enough to offer a good statistical sample of what actually happens in the solar atmosphere. However, the real Sun (or different regions on the Sun) might behave somewhat differently. In this regard, our study does not offer a full answer. Instead, what our analysis does is to establish under which circumstances the electric current density can be correctly determined in the solar photosphere.

5. Effects of photon noise and spectral resolution

In order to study the effect of photon noise, we set up a new inversion in which we analyzed only the spectral lines in the first spectral range in Table 1 without vertical regularization, using HE, and employing four, four, and eight nodes for B_x(z), B_y(z), and B_z(z), respectively. Three different levels of photon noise were employed: σ_p = 0, σ_p = 3 × 10⁻⁴, and σ_p = 10⁻³ (values given in units of the quiet Sun continuum intensity). We refer to these inversions as he5, he6, and he1, respectively. We note that the third of these inversions, with σ_p = 10⁻³, was already described in Sect. 3.2. Results indicate that as the level of noise increases from he5, to he6, and to he1, the floor values, f, become smaller (see Table 4). This means that as the noise increases, the smallest value of the electric current that the inversion can detect increases: he5 can detect electric current densities, ∥j∥, as small as 10^−1.53 ≈ 0.029 A m⁻², whereas the smallest electric current density that he1 can detect is 10^−1.36 ≈ 0.043. The largest increase is seen in j_z, thus indicating that the noise mostly affects the linear polarization profiles Q and U and thereby the determination of B_x and B_y (Borrero & Kobel 2011, 2012). Interestingly, the values of σ do not significantly change as the noise increases (see Table 5).

The fact that the noise does not worsen the determination of the electric current density and the fact that the noiseless case (he5) still yields no clear improvement points at the model’s limitations are the main sources of errors in the determination of the electric currents. By limitations of the model, we refer to the number of nodes, the interpolation method employed between those nodes, and whether they are sufficient to reproduce the z variations of the physical parameters that are present in the MHD simulations. In addition, we need to consider the limitations of radiative transfer itself, whereby variations along the line of sight (i.e., in the z direction) of the physical parameters at scales much smaller than the mean-free path of the photons at a given wavelength do not significantly influence the emerging Stokes spectra (del Toro Iniesta & Ruiz Cobo 2016). This in turn explains why the vertical derivatives of the components of the magnetic field are retrieved the worst (see the bottom panel in Fig. 7). This interpretation is further reinforced by the fact that the more complexity (nodes or free parameters) we allow during the inversion, the better the retrieval of the electric current density is.

We also studied the effect of spectral resolution by performing an inversion that is identical to he1 but convolving the data across the spectral dimension with a transmission profile characterized by a Gaussian function with a full width half maximum (FWHM) of 40 mÅ and later with a FWHM of 80 mÅ. These new inversions are referred to as he7 and he8. The spectral sampling was still 20 mÅ pixel⁻¹ (because only the first spectral range in Table 1 was inverted). We note that, although he1 has an infinite spectral resolution (i.e., the convolving function can be viewed as a zero-width Dirac δ), he7 and he8 possess a spectral resolution of 0.025 and 0.0125 mÅ⁻¹, respectively (i.e., inversely proportional to the width of the convolving function; see the footnote in Sect. 2). Table 4 shows that the floor values barely change as the spectral resolution degrades. The only noticeable change occurs in the floor value for the j_z component of the electric current density, which becomes slightly larger (f = 1.36, 1.41, 1.42) as the spectral resolution degrades. In this way, the he1 inversion allows values of j_z as low as 10^−1.37 ≈ 0.042 A m⁻² to be determined, and he8 allows values as low as 10^−1.42 ≈ 0.038 A m⁻² to be inferred. A similar effect was seen in Sect. 4.1, where the least complex inversion (he1) allowed the lowest values of the electric current density to be correctly inferred. Table 5 also shows that the error in the determination of the electric current density increases very modestly as the spectral resolution degrades, with the notable exception of j_z, where σ increases by as much as 5–10% from he1 to he8.

6. Conclusions

We have studied the accuracy to which electric currents in the lower solar photosphere can be inferred via the application of Stokes inversion codes to the radiative transfer equation for polarized light. Increasingly complex inversions were applied to a set of synthetic observations in multiple spectral lines obtained from radiative three-dimensional MHS simulations. The Stokes inversions range from simple ME-like inversions to regularized z-dependent inversions that include MHS constraints.

While the simplest inversions are capable of retrieving the three components of the electric current density (j_x,j_y,j_z) within a factor of 10 of the values present in the MHD simulations (σ ≈ 1.0 dex), the most complex one improves this to within a factor of 5 (σ ≈ 0.7 dex) of the values present in the MHD simulations. Concerning the modulus ∥j∥ of the electric current density, the accuracy is σ ≈ 0.7 dex for the ME-like inversion (factor of 5) and σ ≈ 0.5 dex (factor of 3) for the MHS inversion. The improvement is not caused by one particular kind of inversion but rather by the cumulative effect of all new layers of complexity (number of spectral lines, number of nodes, using MHS instead of HE, etc.).

Our results suggest that the main sources of errors are the limitations of the inversion models themselves. This is also supported by the fact that neither photon noise nor spectral resolution play an important role in the determination of the electric current density.

The most complex inversion tested in this paper (Sect. 3.6) allows the modulus electric current density to be determined in regions of relatively strong magnetic fields (B >  300 gauss) within a factor of 2 (σ = 0.2 dex) of the currents present in the simulations. This might be accurate enough to allow us to determine regions in the solar photosphere where Joule or Ohmic heating is occurring.

Acknowledgments

This work has received funding from the Deutsche Forschungsgemeinschaft (DFG project number 321818926) and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (SUNMAG, grant agreement 759548). JMB acknowledges travel support from the Spanish Ministry of Economy and Competitiveness (MINECO) under the 2015 Severo Ochoa Program MINECO SEV-2015-0548 and from the SOLARNET project that has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement no 824135. The Institute for Solar Physics is supported by a grant for research infrastructures of national importance from the Swedish Research Council (registration number 2017-00625). Special thanks to K.D. Leka for suggestions that led to the development of the work presented in this paper. This research has made use of NASA’s Astrophysics Data System.

References

All Tables

All Figures

	Fig. 1. Logarithm of the absolute value of the electric current density in the MHD simulations. From left to right: jx, mhd(x, y, τc), jy, mhd(x, y, τc), jz, mhd(x, y, τc), and ∥jmhd(x, y, τc)∥. The electric currents are shown at two different optical depths: τc = 0.1 (upper panels) and τc = 1 (lower panels).
In the text

	Fig. 2. Same as Fig. 1 but showing the electric currents as obtained from the ME-like inversions described in Sect. 3.1 (he1).
In the text

Fig. 3. Logarithmic scatter-density plots of the electric current density from the MHD simulations (horizontal axis) and the values inferred by the me Stokes inversion (vertical axis) at τc = 1 (lower panels) and τc = 0.1 (upper panels). From left to right: logarithm of the absolute value of the ∥jx∥ component of the electric current density, logarithm of the absolute value of the ∥jy∥ component of the electric current density, logarithm of the absolute value of the ∥jz∥ component of the electric current density, and logarithm of the absolute value of the modulus of the electric current density, ∥j∥. The horizontal dashed black line represents the so-called floor value, below which the Stokes inversion is not capable of correctly inferring the electric currents (see text for details). The solid thick diagonal corresponds to the one-to-one retrieval along which the Stokes inversion perfectly retrieves the electric currents. All pixels in the selected region described in Sect. 2 are included in this figure.

In the text

Fig. 4. Histograms of the logarithm of the absolute value of the quotient of the components of the electric current density from the MHD simulations and as obtained from the Stokes inversion: (a) ∥jx∥ (upper left); (b) ∥jy∥ (upper right); (c) ∥jz∥ (lower left); and (d) ∥j∥ (lower right). The histograms correspond to the Stokes inversions described in Sect. 3: red for me, yellow for he1, green for he2, cyan for he3, blue for he4, and purple for mhs.

In the text

	Fig. 5. Same as Fig. 4 but employing only regions where ∥B(τc = 1)∥ >  300 gauss.
In the text

	Fig. 6. Same as Fig. 3 but showing the electric currents as obtained from the inversions under MHS equilibrium described in Sect. 3.6 (mhs).
In the text

Fig. 7. Histograms of the logarithm of the absolute value of the quotient of the spatial derivatives of the components of the electric current density from the MHD simulations and as obtained from the Stokes inversion. Top panel: ∂Bz/∂y (or alternatively ∂Bz/∂x). Middle panel: ∂Bx/∂y (or alternatively ∂By/∂x). Bottom panel: ∂By/∂z (or alternatively ∂Bx/∂z). The histograms correspond to the Stokes inversions described in Sect. 3: me (red), he1 (yellow), he2 (green), he3 (cyan), he4 (blue), and mhs (purple).

In the text

	Fig. 8. Same as Fig. 7 but at τc = 10−2.
In the text