Combining magneto-hydrostatic constraints with Stokes profiles inversions. III. Uncertainty in the inference of electric currents

Electric currents play an important role in the energy balance of the plasma in the solar atmosphere. They are also indicative of non-potential magnetic fields and magnetic reconnection. Unfortunately, the direct measuring of electric currents has traditionally been riddled with inaccuracies. We study how accurately we can infer electric currents under different scenarios. We carry out increasingly complex inversions of the radiative transfer equation for polarized light applied to Stokes profiles synthesized from radiative three-dimensional magnetohydrodynamic (MHD) simulations. The inversion yields the magnetic field vector, ${\bf B}$, from which the electric current density, ${\bf j}$, is derived by applying Ampere's law. We find that the retrieval of the electric current density is only slightly affected by photon noise or spectral resolution. However, the retrieval steadily improves as the Stokes inversion becomes increasingly elaborated. In the least complex case (a Milne-Eddington-like inversion applied to a single spectral region), it is possible to determine the individual components of the electric current density ($j_{\rm x}$, $j_{\rm y}$, $j_{\rm z}$) with an accuracy of $\sigma=0.90-1.00$ dex, whereas the modulus ($\|{\bf j}\|$) can only be determined with $\sigma=0.75$ dex. In the most complicated case (with multiple spectral regions, a large number of nodes, Tikhonov vertical regularization, and magnetohydrostatic equilibrium), these numbers improve to $\sigma=0.70-0.75$ dex for the individual components and $\sigma=0.5$ dex for the modulus. Moreover, in regions where the magnetic field is above 300 gauss, $\|{\bf j}\|$ can be inferred with an accuracy of $\sigma=0.3$ dex. In general, the $x$ and $y$ components of the electric current density are retrieved slightly better than the $z$ component.


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 ). 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 , 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.

MHD simulations and synthetic Stokes profiles
In this paper we employ the same physical parameters as in Pastor . 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: j = µ −1 0 ∇ × B, where j is given in A m −2 when B is given in teslas, and µ 0 = 4π10 −7 kg m s −2 A −2 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, j mhd = j 2 x + j 2 y + j 2 z (rightmost column), in SI units of amperes per m 2 .
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 ) 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 1 , and atomic parameters are given in Table 1. Unless otherwise specified, normally distributed random noise -with a σ pn = 10 −3 (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).

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 . 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.

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), HeLIx + (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 From left to right: j x,mhd (x, y, τ c ), j y,mhd (x, y, τ c ), j z,mhd (x, y, τ c ), and j mhd (x, y, τ c ) . The electric currents are shown at two different optical depths: τ c = 0.1 (upper panels) and τ c = 1 (lower panels). Notes. The first column provides the label used for each spectral region (which might include more than one spectral line). Regarding the symbols: ∆λ is the pixel size in mÅ; n λ is the number of wavelengths used for that spectral line; λ 0 is the central wavelength for the electronic transition associated with the spectral line; e low and e upp are the electronic configurations of the lower and upper energy level, respectively; and E low is the excitation potential (in eV) of the lower energy level. These atomic data are adopted from Nave et al. (1994). The α and σ/a 2 0 are the velocity exponent and collision cross-section parameters (in units of Bohr's radius, a 0 ), respectively, as defined in the Anstee, Barklem, and O'Mara collision theory for the broadening of metallic lines by neutral hydrogen collisions (Anstee & O'Mara 1995;Barklem & O'Mara 1997;Barklem et al. 1998). 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.

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).

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.

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.

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 χ 2 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).

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.

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. 2). 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 −2 . This feature is seen not only in me, but in all tested Stokes inversions.
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 −2 (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 −2 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 2 .
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 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 −2 (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 2 . 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 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 j x component of the electric current density, logarithm of the absolute value of the j y component of the electric current density, logarithm of the absolute value of the j z 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. Table 3. Values of σ calculated from Fig. 4 (see text for details) at two different optical depths (τ c = 1, 0.1) for the different Stokes inversions described in Sect. 3. For each Stokes inversion type, we provide two values (e.g., 0.85/0.70). The first one (e.g., 0.85) corresponds to the standard deviation, σ, obtained from the histogram that includes all pixels in the simulations (Fig. 4), while the second one (e.g., 0.70) is the standard deviation, σ, obtained from the histogram that includes only pixels where B(τ c = 1) > 300 (Fig. 5).  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. 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).

An analysis of
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.

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: 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. 2 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.
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. 2. 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. 2). 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.

Height dependence
In Sect. 4.1 we pointed out that the inference of electric currents quickly worsens for optical depths τ c < 10 −2 . We can offer two different explanations as to why this happens. The first is that for τ c < 10 −2 (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 upperphotosphere (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 −2 and quickly drops for smaller optical depths (Kuckein et al. 2012;Felipe et al. 2016;Shchukina & Trujillo Bueno 2019).
The second reason is that, even if the spatial derivatives of the magnetic field could be determined for τ c = 10 −2 (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 −2 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 −3 , but higher up, at τ c = 10 −4 , 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 −4 , 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 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: ∂B z /∂y (or alternatively ∂B z /∂x). Middle panel: ∂B x /∂y (or alternatively ∂B y /∂x). Bottom panel: ∂B y /∂z (or alternatively ∂B x /∂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 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.

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 −4 , and σ p = 10 −3 (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 −3 , 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 −2 , 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. 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 −1 (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Å −1 , 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 −2 to be determined, and he8 allows values as low as 10 −1.42 ≈ 0.038 A m −2 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.

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 Table 4. Powers of ten, f , that yield the floor values as 10 − f , at two different optical depths (τ c = 1,τ c = 0.1) for the different inversions carried out in Sect. 5. 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. Table 5. Standard deviation in the inference of the electric current density, σ, at two different optical depths (τ c = 1, 0.1) for the different Stokes inversions described in Sect. 5.