© The Authors 2024

1. Introduction

Stokes inversion codes are one of the most useful tools to study the magnetic field in the lower solar atmosphere (photosphere and chromosphere). These codes typically operate by solving the radiative transfer equation for polarized light using an atmospheric model that includes the physical parameters of the plasma (temperature, density, magnetic field, etc.) in order to produce a theoretical or synthetic Stokes vector, 𝕀^syn, in absorption and emission spectral lines. The synthetic Stokes profiles are then compared with the observed one 𝕀^obs and, in case of a mismatch, the physical parameters of the atmospheric model are modified in order to produce a better fit to the aforementioned observations. For a review of the different methods available to achieve this, we refer the reader to del Toro Iniesta & Ruiz Cobo (2016).

Although originally proposed more than two decades ago by Socas-Navarro (2001), in recent years there has been a renewed interest in combining the results from Stokes inversion codes with physical constraints arising from magnetohydrodynamic (MHD) and Maxwell’s equations (Puschmann et al. 2010; Riethmüller et al. 2017; Borrero et al. 2019), as we are now able to retrieve the magnetic field in the three-dimensional Cartesian domain B(x, y, z) in a way that is consistent with the magneto-hydrostatic equilibrium (Pastor Yabar et al. 2021; Borrero et al. 2021). In spite of this success, the inferred magnetic field remains inconsistent with Maxwell’s second equation that states that the magnetic field must be solenoidal (also known as Gauss’ law for the magnetic field): ∇ ⋅ B = 0. Because of this, the magnetic field inferred from the application of Stokes inversion codes to spectropolarimetric observations is sometimes put into question. In particular, Balthasar (2018) has pointed that Stokes inversion applied to spectropolarimetric data from sunspots yields a vertical variation of B_z that is a factor of five to ten larger than the horizontal variation of B_x or B_y: dB_z/dz ≈ −3 G km⁻¹, whereas dB_x/dx ≈ 0.5 G km⁻¹. On the other hand, some authors argue that what spectropolarimetry allows to be measured is H instead of B and that ∇ ⋅ H is not necessarily zero (Bommier 2020). Alternatively, minimizing the divergence of the magnetic field vector can be used to solve the 180-degree ambiguity in the components of the magnetic field perpendicular to the observer’s line of sight (Metcalf 1994; Metcalf et al. 2006) or to establish a common geometrical height scale (Löptien et al. 2018).

In this work, we put forward the idea that the reason the null divergence condition of the magnetic field, as inferred from spectropolariemtric observations, is not satisfied is because of the intrinsic uncertainties in the determination of the different components of the magnetic field. If correct, this interpretation would imply that, within the error margins in the inference of B_x, B_y, and B_z, it should be possible to propose a fully solenoidal magnetic field (∇ ⋅ B = 0) that can fit the observed polarization signals. To test this hypothesis, we applied the FIRTEZ Stokes inversion code (Pastor Yabar et al. 2019) to spectropolarimetric observations of a sunspot located at the disk’s center observed with the Hinode/SP instrument (Sect. 2), and we show that within the errors of the inversion code, we can obtain a solenoidal magnetic field (Sect. 3) that is also capable of fitting the observed polarization signals (Sect. 4) in about 80% of the analyzed three-dimensional domain. Correlations between the fits to the Stokes parameters and the inferred magnetic field are studied in Sect. 5. The effects of imposing solenoidal magnetic fields in the inference of ∇ × B (i.e., electric currents) is investigated in Sect. 6. Finally, in Sect. 7, we present our conclusions and implications for the study of magnetic fields in the solar photosphere via Stokes inversion techniques.

2. Observations and Stokes inversion

In this work, we analyze the sunspot AR 10933 observed at disk center on January 6, 2007, between 00:00 and 01:25 UT with the Hinode spectropolarimeter (SP; Lites et al. 2001; Ichimoto et al. 2007). The Hinode/SP is attached to the Solar Optical Telescope (SOT; Suematsu et al. 2008; Tsuneta et al. 2008; Shimizu et al. 2008 on board the Japanese satellite Hinode (Kosugi et al. 2007). The observed data comprise the Stokes vector 𝕀^obs = (I, Q, U, V) in a spectral region around 630 nm that contains two Fe I lines across a total of N_λ = 112 wavelength positions with a wavelength sampling of about 21.5 mÅ. The atomic parameters for these spectral lines can be found in Borrero et al. (2014, see their Table 1). The SP is a slit-spectrograph where a given region is scanned spatially. For each slit position, the light is integrated for a total of 4.8 s, yielding a photon noise level of about σ = 10⁻³ in units of the quiet Sun continuum intensity. A map of the continuum intensity I_c, normalized to the average quiet Sun continuum intensity I_c, qs, can be seen in Fig. 1. This map includes n_x = 350 and n_y = 300 pixels on the horizontal (x, y) plane, with an spatial sampling of about 0.16 arcsec (i.e., dx = dy = 120 km at disk center).

Fig. 1. Map of the continuum intensity at 630 nm in AR 10933 observed at disk center on January 6, 2007, by SP on board the Hinode spacecraft. Physical parameters on the (x, z) slice indicated by the horizontal red line are studied in Sect. 3.

The observed Stokes vector is analyzed using the FIRTEZ Stokes inversion code (Pastor Yabar et al. 2019) with the following number of free parameters: eight for the temperature T, eight for the line-of-sight velocity (i.e., or vertical velocity at disk center) v_z, and four for each of the three components of the magnetic field B_x, B_y, and B_z, for a total of F = 28 free parameters. The free parameters are located equidistantly along the z direction, which in turn is discretized along n_z = 128 points. The grid size along this direction is dz = 12 km. During the inversion, all four Stokes parameters are given the same weight (w_i = w_q = w_u = w_v = 1) in the definition of the χ² merit function that is being minimized:

$$\begin{array}{c}\hfill {\chi }^2=\frac{1}{4N_{\lambda }-F}\sum _{k=1}^4\sum _{j=1}^{N_{\lambda }}{[{\mathbb{I}}_k^{\mathrm{obs}}({\lambda }_j)-{\mathbb{I}}_k^{\mathrm{syn}}({\lambda }_j)]}^2{\left(\frac{w_k}{\sigma }\right)}^2,\end{array}$$(1)

where the denominator 4N_λ − F refers to the number of degrees of freedom, that is, the difference between the number of data points and the number of free parameters. The index k refers to each of the four Stokes parameters I, Q, U, and V so that 𝕀₁ = I, 𝕀₂ = Q, and so forth. The term σ refers to the photon noise. We note that in this work, we are not performing any regularization, and therefore no additional terms appear in Eq. (1). As a result of the inversion, FIRTEZ provides the aforementioned physical parameters in the three-dimensional Cartesian domain (x, y, z). This is possible because FIRTEZ calculates the gas pressure P_g employing magneto-hydrostatic equilibrium instead of hydrostatic equilibrium (Borrero et al. 2019, 2021). This is done after the Stokes inversion runs at every (x, y) pixel by the MHS module that iteratively solves the following Poisson-like equation:

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2(ln{P}_{\mathrm{g}})=-\frac{\mathit{ug}}{K_{\mathrm{b}}}\frac{\partial }{\partial z}\left[\frac{\mu }{T}\right]+\frac{1}{c}\mathrm{\nabla }·\left[\frac{\mathit{j}\times \mathit{B}}{P_{\mathrm{g}}}\right]·\end{array}$$(2)

We note that FIRTEZ has an internal implementation of the non-potential field calculation method (Georgoulis 2005), so it provides the magnetic field without the 180° ambiguity in the B_x and B_y components of the magnetic field. The inversion performed under the conditions described above will be referred to as MHS (or simply as non-solenoidal) inversion, and the inferred magnetic field will be referred to as B_MHS. A graphical illustration of how the FIRTEZ code operates can be found in Fig. 2 in Borrero et al. (2021). Following this inversion, we performed a second inversion using as the initial guess the results from the MHS inversion, where B_z is no longer inverted, thus having only F = 24 free parameters. However, the Stokes V was still fitted: w_v = 1. Instead of inferring B_z from the Stokes inversion, it was obtained by considering that the magnetic field must be solenoidal:

$$\begin{array}{c}\hfill \frac{\partial B_z}{\partial z}=-(\frac{\partial B_x}{\partial x}+\frac{\partial B_y}{\partial y})·\end{array}$$(3)

We numerically solved this differential equation via finite differences, where the spatial derivatives were approximated with a second-order centered difference formula for the inner grid points along the z direction and using second-order forward and backward difference formulas at the boundaries. Given B_x and B_y, this led to a linear system of equations that could be easily solved via a matrix inversion, yielding B_z(z). The process was then repeated for every (x, y) point, thereby giving B_z(x, y, z). We note that since Eq. (3) is a first-order differential equation, only one boundary condition was needed. We selected neither the upper (z = z_max = n_z  *  dz = 1536 km) nor the lower (z = 0) boundaries because these are typically located in a region where the spectral lines carry no useful information about the magnetic field. Instead, we used as a boundary condition the value of the vertical component of the magnetic field from the previous MHS (non-solenoidal) inversion at a location where the employed spectral lines observed by Hinode/SP provide the largest sensitivity: B_z, MHS(log τ_c = −1.5). From there, we performed two integrations of Eq. (3): one downward from z(log τ_c = −1.5) until z = 0 and another upward from z(log τ_c = −1.5) until z = z_max. Because at z(log τ_c = −1.5) we employed second-order forward or backward finite differences, which are known to be slightly less accurate than centered differences, the accuracy to which ∇ ⋅ B = 0 can be achieved at this height is slightly worse than above and below this point: ∥∇⋅B∥≈10⁻⁴ vs ∥∇⋅B∥≈10⁻⁶. Next, we note that according to Eq. (3), once B_z(z) is obtained, any value can be added to B_z(z) as long as it is constant with z. We took advantage of this fact to further refine B_z(z) so that its average value between log τ_c = [0, 3] (i.e., the region where the line is formed) is the same as in the MHS inversion. Together with our selection of the boundary condition, this ensured that the Stokes V signals produced by the newly obtained B_z(z) did not deviate too much from the observed ones since B_z, MHS(z) had been previously obtained in order to fit the observed Stokes V signals.

The calculation of B_z in the fashion described above does not take place at every iteration step of the Stokes inversion itself but rather at the end, once it has been run for all (x, y) pixels. This is done because it is only at this stage that the 180-degree ambiguity in (B_x, B_y) has been corrected (see Fig. 2 in Borrero et al. 2021) and that B_x(x, y, z) and B_y(x, y, z) are known, thus allowing the right-hand side of Eq. (3) to be determined. This choice has the added benefit that the MHS module receives a fully solenoidal magnetic field, and therefore, the electric currents (j ∝ ∇ × B) and the Lorentz force (L ∝ j × B) used to establish the magnetic-hydrostatic equilibrium arise from a more realistic magnetic field (see also Sect. 6). Now, owing to the existence of several cycles between the Stokes inversion module and the MHS module, B_z(z) is adjusted several times to reflect the possible changes in B_x and B_y that might occur and that affect the right-hand side of Eq. (3). If the linear polarization signals (Q and U) are strong, convergence in B_z(z) is achieved after only one or two cycles between the Stokes inversion and the MHS module. If Q and U are weak, convergence cannot be ensured because B_x and B_y obtained from the Stokes inversion change substantially after each cycle, and therefore the right-hand-side of Eq. (3) changes somewhat randomly. In these cases, our choice of boundary condition still ensures that the average value B_z over the region where the spectral line is formed stays the same, but the gradient ∂B_z/∂z is not well constrained. The second inversion described here is hereafter referred to as MHS + ∇ ⋅ B = 0 or simply as “solenoidal inversion”.

As proof that the magnetic field obtained by the second inversion is indeed solenoidal, we present in Fig. 2 the histogram of ∇ ⋅ B for the first MHS inversion (blue) and the second MHS + ∇ ⋅ B = 0 inversion (red). As can be seen, the second inversion retrieves a magnetic field that satisfies the null divergence condition up to the numerical precision given by the order of the finite difference used to solve Eq. (3). A logarithmic scale is used because the red curve appears as an almost perfect δ-Dirac in a linear scale. We note that unlike Puschmann et al. (2010) and (Löptien et al. 2018), we do not minimize ∇ ⋅ B but rather strictly impose it. Figure 3 illustrates the spatial distribution of ∥∇⋅B∥ across the sunspot for the MHS (non-solenoidal) inversion (top) and the solenoidal inversion (bottom) at a height of 100 km above the height where the continuum is formed: z = z(τ_c = 1)+100 km. As one can see, the value of the divergence of the magnetic field vector in the second inversion is about five to six orders of magnitude smaller than in the first inversion. Similar results were obtained at different heights as well.

Fig. 2. Histograms of the divergence of the magnetic field from the first MHS inversion (blue) and second MHS + ∇ ⋅ B = 0 inversion (red).

Fig. 3. Horizontal (x, y) maps of the logarithm of the absolute value of the divergence of the magnetic field at a height of 100 km above the continuum forming layer. Top panel: results from the first MHS inversion. Bottom panel: results from the second MHS + ∇ ⋅ B = 0 inversion. Units of ∇ ⋅ B are in Gauss per kilometer.

3. Results: Inference of B_z(z)

In order to compare the topology of the magnetic field, in particular of B_z since this is the one that the solenoidal inversion determines differently from the MHS (non-solenoidal) inversion, we present in Fig. 4 the vertical component of the magnetic field in the (x, z) plane for the slice indicated by the red line in Fig. 1. The reason for choosing this region is because it is in the penumbra, where the magnetic field of a sunspot is most inhomogeneous, and therefore it is here that the two inversions performed in this work could potentially lead to more different results. Interestingly, both inversions show the well-known spine and intraspine structure of the penumbral magnetic field (Lites et al. 1993; Martinez Pillet 1997; Borrero et al. 2008), whereby weak horizontal fields (i.e., intraspines) are embedded in stronger and more vertical fields (i.e., spines). This is also often referred to as the uncombed penumbral structure (Solanki & Montavon 1993; Borrero et al. 2007; Tiwari et al. 2013). In this case, the inclusion of the constraints imposed by Gauss’ law for magnetism does not lead to significant changes in the overall structure of the magnetic field.

Fig. 4. Vertical component of the magnetic field, Bz(x, z), along the slice indicated by the red line in Fig. 1. Only regions between log τc ∈ [1, −4] (approximate region of sensitivity of the spectral lines observed by Hinode/SP) are shown. The dashed black line shows the location of the z(τc = 1) level (i.e., Wilson depression). Results from the MHS inversion are shown in the top panel, whereas results from the MHS + ∇ ⋅ B =  0 inversion are displayed in the bottom panel.

A more detailed comparison can be carried out by focusing on two positions along the x axis that roughly correspond to the locations of the center of a spine (x ≈ 18.3 Mm) and the center of an intraspine (x ≈ 19.5 Mm). In Fig. 4, these two locations are marked by the vertical dashed white lines and labeled as s (red; spine) and i (blue; intraspine). Figure 5 displays B_z(z) for these two locations. Along the intraspine, both the MHS inversion (solid blue) and solenoidal inversion (dashed steel blue) yield dB_z/dz > 0. In fact, both results agree extremely well within the 3σ error bars from the MHS inversion (shaded cyan area). We emphasize here that FIRTEZ finds the standard deviation σ in the inference of the different physical parameters by diagonalizing the modified Hessian matrix ℍ′ corresponding to the atmospheric model that best fits the observed Stokes vector. For instance, for a physical parameter denoted by index j (see Appendix B or Chapter 11.2.1 in Sánchez Almeida 1997; del Toro Iniesta 2003, respectively),

$$\begin{array}{c}\hfill {\sigma }_j^2=\frac{{\chi }^2}{F}{({\mathbb{H}}^{\prime -1})}_{\mathit{jj}},\end{array}$$(4)

Fig. 5. Vertical stratification of the vertical component of the magnetic field, Bz(z), along the two positions in Fig. 4 (see vertical white dashed lines). The red and orange colors correspond to the results along the spine from the MHS and MHS + ∇ ⋅ B = 0 inversion, respectively. The shaded pink areas indicate the 3σBz confidence level from the MHS inversion. Similarly, the blue and steel blue colors correspond to the results along the intraspine from the MHS and MHS + ∇ ⋅ B = 0 inversion, respectively. The shaded cyan areas indicate the 3σBz confidence level from the MHS inversion.

where the modified Hessian matrix ℍ′ is nearly the same as the regular Hessian matrix except that the diagonal elements are multiplied by a variable factor that controls whether the inversion is far or close to a minimum, as prescribed by the Levenberg-Marquardt algorithm (Levenberg 1944; Press et al. 1986). In the case of the spine, the MHS (non-solenoidal) inversion yields a vertical component of the magnetic field B_z that is almost constant with height (solid red line in Fig. 5), whereas the solenoidal inversion retrieves a B_z(z) that decreases by about 1000 G in 500 km (dashed orange line) and lies well beyond the 3σ error bars from the MHS inversion. We note that the vertical gradient in this case is about dB_z/dz ≈ −2 G km⁻¹. This large vertical gradient has been deemed to be at odds with the ∇ ⋅ B condition for the magnetic field (Balthasar 2018), and this is probably correct if it should be present through the entire penumbra. However, our inversion demonstrates that locally, this large gradient is not only consistent with but actually demanded by the null divergence condition of the magnetic field.

Based on the results from the previous two examples (spine and intraspine), we subsequently aimed at determining how often the vertical component of the magnetic field provided by the MHS + ∇ ⋅ B = 0 inversion is consistent with B_z as inferred from the regular MHS (non-solenoidal) inversion. In order to investigate this, we show in Fig. 6 the cumulative histogram of the grid cells at different log τ_c levels as a function of ΔB_z/σ_Bz, where σ_Bz is again the standard deviation in the determination of the vertical component of the magnetic field in the MHS inversion (Eq. (4)). These histograms take into account only those regions where the horizontal component of the magnetic field $B_{\mathrm{h}}=\sqrt{B_x^2+B_{\mathit{y}}^2}$ at log τ_c = −1.5 is larger than 300 G. This was done so as to avoid regions where B_z is not well constrained via the horizontal magnetic fields (see discussion in Sect. 2). Figure 6 shows that in 80% of the three-dimensional (x, y, z) domain, B_z from the solenoidal inversion is within B_z ± 3σ_Bz of the regular MHS (non-solenoidal) inversion. From a statistical point of view, we can interpret this result in the following way: The magnetic field inferred from the solenoidal inversion is, overall, similar to the magnetic field provided by the non-solenoidal inversion, and therefore previous results obtained regarding the topology of the magnetic field in sunspots based on the inversion of spectropolarimetric data (see Borrero & Ichimoto 2011, for a review of those results) will likely stand. However, the 3σ tolerance is large enough, and the fraction of the three-dimensional domain where this is not satisfied is also large enough (≈20%) that this statement should be confirmed in a future work and in a case-by-case basis.

Fig. 6. Cumulative histogram of the number of pixels, shown at three different optical depth levels (red: log τc = 0; blue: log τc = −1.5; green: log τc = −3), as a function of the ratio of the difference between Bz obtained from both inversions to the standard deviation σBz. Only pixels where Bh(log τc = −1.5) > 300 G are considered (see text for details).

4. Results: Fits to the polarimetric signals

In order to accept the results from the MHS + ∇ ⋅ B = 0 inversion as reliable, it is customary to show that they can fit the polarization signals to a degree that is at least comparable to the fits produced by the regular MHS inversion. Otherwise, whatever physical realism we might gain from including the constraint imposed by Maxwell’s equation will be defeated by the fact that it is unable to explain the observed Stokes vector.

To answer this question, in the context of the selected spine and intraspine, we present in Fig. 7 the observed (circles) and best-fit profiles (color lines) for both MHS and solenoidal inversions. In addition, we also provide the χ² value for each Stokes parameter separately. At first glance it is clear that both inversions fit the observed Stokes profiles very well and that the fits are almost indistinguishable for both kinds of inversions. We emphasize, however, that B_z was not being explicitly inverted in the case of the MHS + ∇ ⋅ B = 0 inversion (see Sect. 2), and therefore, it is somewhat surprising that Stokes V is fitted with such accuracy also in this inversion. Perhaps this was to be expected in the case of the intraspine (blue and steel blue lines in the right panels of Fig. 7) because B_z(z) was almost identical in the MHS and solenoidal inversions (see Fig. 5), but it comes as a surprise that, qualitatively, the fits to Stokes V are also so close in the case of the spine (red and orange lines in the left panels in Fig. 7) since both inversions feature a rather different B_z(z), that is, they differ beyond the 3σ level. This has been possible thanks to our choice of boundary conditions when integrating Eq. (3) (see discussion on Sect. 2). We note, however, that quantitatively the fits to Stokes V are slightly worse in the solenoidal inversion than in the MHS (non-solenoidal) inversion: χ_v = 3.81 vs. χ_v = 2.92, respectively.

Fig. 7. Comparison between observed and fitted Stokes profiles. Left panels: observations (black circles) and fits (color lines) at the location of the spine. Right panels: observations (black circles) and fits (color lines) at the location of the intraspine. From top to bottom: Stokes I, Q, U, and V. Solid color lines (red and blue) are the fits produced by the MHS inversion, whereas dashed color lines (orange and steel blue) correspond to the fits produced by the MHS ∇ ⋅ B = 0 inversion. The χ2 values for each case are also provided. The wavelength scale is given with respect to the first Fe I line at λ0 = 6301.5012 Å.

To investigate whether the results for the two aforementioned examples also apply to the full sunspot, we show in Fig. 8 histograms of the difference between the χ² achieved by the MHS and the solenoidal inversion. Once more, we restricted ourselves to regions where the horizontal component of the magnetic field B_h(log τ_c = −1.5) is larger than 300 G. The plot in the figure shows that, overall, the fit to Stokes I (red) improves after applying the solenoidal inversion. This comes at the expense of a worsening of the fit to Stokes V (blue). This was to be expected since B_z(z) is not inverted at all during the solenoidal inversion. Regarding the linear polarization profiles (Stokes Q and U), both inversions fit the observations equally well. Overall, the total χ² (black) also improves after application of the MHS + ∇ ⋅ B = 0 inversion, with the mean χ² over the entire map being about half of the mean in the MHS inversion. We therefore concluded that, in general, the MHS + ∇ ⋅ B = 0 inversion fits the observed polarization signals at least with the same degree of accuracy as the MHS inversion in spite of featuring a smaller number of free parameters (i.e., B_z is not inverted). Because of these reasons, along with the fact that the magnetic field is more physically grounded, one could consider the MHS + ∇ ⋅ B = 0 inversion preferable to the MHS (non-solenoidal) inversion. We emphasize, however, that the former is not a substitute for the latter because the MHS inversion is needed in order to initialize the solenoidal one. Had we not used the boundary condition at B_z(log τ_c = −1.5) from the MHS inversion and later shifted B_z(z) by a constant value to match the average from the MHS inversion between log τ_c = [0, −3] (see Sect. 2), the solenoidal inversion would not have been able, by itself, to fit the observed Stokes profiles.

Fig. 8. Histogram of the differences of the χ2 normalized to the average ${\overline{\chi }}^2$ achieved by the MHS + ∇ ⋅ B = 0 inversion and by the MHS inversion. Color lines refer to Stokes I (red), Stokes V (blue), Stokes Q (green), Stokes U (orange), and all four Stokes parameters (black). Only pixels where Bh(log τc = −1.5) > 300 G are considered (see text for details).

5. Correlations in B_z(z) and ${\chi }_{\text{v}}^2$

In Sect. 4, we showed that the solenoidal inversion yields a fit to Stokes V that is on average worse than the fit by the MHS inversion (see blue line in Fig. 8). As we have explained, this was to be expected because in the former inversion, B_z is not inverted. In addition to this, in Sect. 3, we showed that in about 20% of the three-dimensional (x, y, z) domain studied, the solenoidal inversion retrieves a vertical component of the magnetic field that is beyond the 3σ_Bz uncertainty bars from the MHS inversion. This begs the question as to whether these two effects are correlated, that is, whether in those regions where B_z(z) from the MHS + ∇ ⋅ B = 0 inversion differs significantly from the B_z(z) in the MHS inversion, the fit to the circular polarization by the solenoidal inversion worsens.

In order to address this question, we present in the left panel of Fig. 9, the spatial distribution of the difference between the χ² in Stokes V from the solenoidal and MHS (non-solenoidal) inversions normalized to the average of the two. Positive regions (red) are those where the MHS inversion fits the circular polarization better. As expected from Fig. 8, positive values should dominate over negative ones. In the right panel of Fig. 9, we show a map of the sunspot where colors indicate the percentage of pixels along the z direction, but inside the line-formation region (log τ_c ∈ [0, −3]), where the vertical magnetic field inferred from the MHS + ∇ ⋅ B = 0 inversion falls within [B_z, MHS − 3σ_Bz, B_z, MHS + 3σ_Bz] in the MHS inversion. The gray regions are those where 80 to 100% of the grid cells along the vertical direction meet this criteria, meaning that B_z(z) is very similar in both the solenoidal and MHS (non-solenoidal) inversions. The blue and orange regions are those where the opposite happens. Comparing the left and right panels in Fig. 9 answers the above question in the affirmative: Regions where B_z(z) from the solenoidal inversion differ most from that of the MHS inversion are the same regions where the fit to the observed Stokes V signals worsens the most. These regions are predominantly located in the outer penumbra, close to the boundary with the quiet Sun. It is important to emphasize that in spite of producing worse fits to Stokes V, the overall χ² improves in the MHS + ∇ ⋅ B = 0 inversion (see black lines in Fig. 8).

Fig. 9. Correlations between Bz(z) and ${\chi }_{\text{v}}^2$. Left panel: difference in the goodness of the fit to Stokes V (circular polarization) by the solenoidal inversion and the MHS inversion normalized to the average of the two. Right panel: percentage of the region (along the z direction) where the observed spectral lines are formed in a manner in which the vertical component of the magnetic field inferred from the solenoidal and MHS inversions agree with the 3σ threshold.

6. Electric currents

An important question to be addressed here is whether the electric currents as inferred from the regular MHS inversion and from the MHS + ∇ ⋅ B = 0 inversion, and given by ∇ × B, differ significantly. In order to do so, we present in Fig. 10 a comparison of the three components of the rotational of the magnetic field at z = z(τ_c = 1)+100 km (i.e., 100 km above the continuum level). These figures show that both inversions yield electric currents that are spatially very consistent with each other. Similar levels of consistency were also obtained at other heights, provided that they are located within the region where the line is formed, namely, z = z(τ_c = 1) or z = z(τ_c = 1)+200 km (not shown). In spite of this spatial consistency, currents determined from each inversion can present large differences of up to one order of magnitude when looking at individual locations. One might argue that the electric currents arising from the solenoidal inversion are to be preferred over the electric currents provided by the MHS inversion since they are obtained from a more physical magnetic field. However it remains to be seen which of the two inversions actually infers the electric currents more reliably. In order to conduct such an investigation, one needs to apply the aforementioned inversions to synthetic Stokes profiles produced by MHD simulations where the electric currents are known beforehand (Pastor Yabar et al. 2021; Borrero & Pastor Yabar 2023). This, however, falls outside the scope of the present paper.

Fig. 10. Three components of ∇ × B obtained from the MHS inversion (top panels) and MHS + ∇ ⋅ B = 0 inversion (bottom panels) at a height that corresponds to 100 km above the continuum level: z = z(τc = 1)+100 km. For visualization purposes, [∇×B]z (right panels) has been multiplied by two. The color scale saturates at ±4 G km−1. Only pixels where Bh(log τc = −1.5) > 300 G are considered (see text for details).

7. Conclusions

We have developed a new method that allows Stokes inversion codes for the radiative transfer equation to infer a magnetic field on the solar atmosphere that verifies the null divergence condition imposed by Maxwell’s equation: ∇ ⋅ B = 0. Unlike previous methods (Puschmann et al. 2010; Löptien et al. 2018), we do not minimize the divergence of the magnetic field vector, but instead we strictly impose it. This is done by determining B_z from B_x and B_y such that the magnetic field is divergence free. The accuracy of the solenoidal solution is only limited by the order of the finite differences formula employed. When applied to spectropolarimetric observations of a sunspot recorded with the SP instrument on board Hinode, we easily achieved ∥∇⋅B∥≤10⁻⁵ − 10⁻⁶ G km⁻¹. This method is now implemented and available to users of the FIRTEZ Stokes inversion code¹.

We have also demonstrated that this new inversion method is capable of fitting the observed Stokes profiles with a very similar degree of fidelity as the inversion method where the magnetic field is not solenoidal. Interestingly, the vertical component of the magnetic field in the solenoidal case agrees in about 80% of the analyzed (x, y, z) domain with the non-solenoidal one within the 3σ confidence level regarding the determination of the B_z. This means that, the overall magnetic topology in the solenoidal inversion is similar to the non-solenoidal one. We interpret these results as proof that regular (i.e., non-solenoidal) Stokes inversions fail to fulfill the ∇ ⋅ B = 0 criterion mostly as a consequence of the uncertainties in the determination of the individual components of the magnetic field. In about 20% of the (x, y, z) domain, the solenoidal magnetic field differs from the non-solenoidal one beyond the 3σ threshold. These deviations are correlated with a worsening of the ability of the solenoidal inversion to fit the observed Stokes V signals compared to the non-solenoidal inversion. It is in these regions where our newly developed method can uncover new insights about the topology of the magnetic field. We leave these investigations for a future work.

An important limitation of the method described in this paper is that it requires strong polarization signals, mainly in the linear polarization, because B_z is obtained from B_x and B_y. This is certainly the case of sunspots, but it remains to be seen how this method performs in other regions, such as network, internetwork, and plage regions. Possible improvements to the current implementation of our method might be obtaining a different component of the magnetic field out of the other two or even adjusting which component is being fixed by the ∇ ⋅ B = 0 condition on a pixel-by-pixel basis, depending on the polarization signals observed there.

Acknowledgments

Authors would like to thanks Dr. Ivan Milić for suggestions and carefully reading the manuscript. J.M.B. would like to thank K.D. Leka, H. Balthasar, and V. Bommier for bringing up to our attention, during the SPW9 workshop held in Göttingen (Germany) during 2019, some of the issues addressed in this work. BRC has been partially funded by AEI/MCIN/10.13039/501100011033/PID2021-125325OB-C55, co-funded by European Regional Development Fund (ERDF) “A way of making Europe”. A.P.Y. is funded/co-funded by the European Union (ERC, MAGHEAT, 101088184). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. This research has made use of NASA’s Astrophysics Data System. Hinode is a Japanese mission developed and launched by ISAS/JAXA, collaborating with NAOJ as a domestic partner, NASA and STFC (UK) as international partners. Scientific operation of the Hinode mission is conducted by the Hinode science team organized at ISAS/JAXA. This team mainly consists of scientists from institutes in the partner countries. Support for the post-launch operation is provided by JAXA and NAOJ (Japan), STFC (UK), NASA, ESA, and NSC (Norway).

References

All Figures

	Fig. 1. Map of the continuum intensity at 630 nm in AR 10933 observed at disk center on January 6, 2007, by SP on board the Hinode spacecraft. Physical parameters on the (x, z) slice indicated by the horizontal red line are studied in Sect. 3.
In the text

	Fig. 2. Histograms of the divergence of the magnetic field from the first MHS inversion (blue) and second MHS + ∇ ⋅ B = 0 inversion (red).
In the text

	Fig. 3. Horizontal (x, y) maps of the logarithm of the absolute value of the divergence of the magnetic field at a height of 100 km above the continuum forming layer. Top panel: results from the first MHS inversion. Bottom panel: results from the second MHS + ∇ ⋅ B = 0 inversion. Units of ∇ ⋅ B are in Gauss per kilometer.
In the text

Fig. 4. Vertical component of the magnetic field, Bz(x, z), along the slice indicated by the red line in Fig. 1. Only regions between log τc ∈ [1, −4] (approximate region of sensitivity of the spectral lines observed by Hinode/SP) are shown. The dashed black line shows the location of the z(τc = 1) level (i.e., Wilson depression). Results from the MHS inversion are shown in the top panel, whereas results from the MHS + ∇ ⋅ B =  0 inversion are displayed in the bottom panel.

In the text

Fig. 5. Vertical stratification of the vertical component of the magnetic field, Bz(z), along the two positions in Fig. 4 (see vertical white dashed lines). The red and orange colors correspond to the results along the spine from the MHS and MHS + ∇ ⋅ B = 0 inversion, respectively. The shaded pink areas indicate the 3σBz confidence level from the MHS inversion. Similarly, the blue and steel blue colors correspond to the results along the intraspine from the MHS and MHS + ∇ ⋅ B = 0 inversion, respectively. The shaded cyan areas indicate the 3σBz confidence level from the MHS inversion.

In the text

	Fig. 6. Cumulative histogram of the number of pixels, shown at three different optical depth levels (red: log τc = 0; blue: log τc = −1.5; green: log τc = −3), as a function of the ratio of the difference between Bz obtained from both inversions to the standard deviation σBz. Only pixels where Bh(log τc = −1.5) > 300 G are considered (see text for details).
In the text

Fig. 7. Comparison between observed and fitted Stokes profiles. Left panels: observations (black circles) and fits (color lines) at the location of the spine. Right panels: observations (black circles) and fits (color lines) at the location of the intraspine. From top to bottom: Stokes I, Q, U, and V. Solid color lines (red and blue) are the fits produced by the MHS inversion, whereas dashed color lines (orange and steel blue) correspond to the fits produced by the MHS ∇ ⋅ B = 0 inversion. The χ2 values for each case are also provided. The wavelength scale is given with respect to the first Fe I line at λ0 = 6301.5012 Å.

In the text

	Fig. 8. Histogram of the differences of the χ2 normalized to the average ${\overline{\chi }}^2$ achieved by the MHS + ∇ ⋅ B = 0 inversion and by the MHS inversion. Color lines refer to Stokes I (red), Stokes V (blue), Stokes Q (green), Stokes U (orange), and all four Stokes parameters (black). Only pixels where Bh(log τc = −1.5) > 300 G are considered (see text for details).
In the text

Fig. 9. Correlations between Bz(z) and ${\chi }_{\text{v}}^2$. Left panel: difference in the goodness of the fit to Stokes V (circular polarization) by the solenoidal inversion and the MHS inversion normalized to the average of the two. Right panel: percentage of the region (along the z direction) where the observed spectral lines are formed in a manner in which the vertical component of the magnetic field inferred from the solenoidal and MHS inversions agree with the 3σ threshold.

In the text

Fig. 10. Three components of ∇ × B obtained from the MHS inversion (top panels) and MHS + ∇ ⋅ B = 0 inversion (bottom panels) at a height that corresponds to 100 km above the continuum level: z = z(τc = 1)+100 km. For visualization purposes, [∇×B]z (right panels) has been multiplied by two. The color scale saturates at ±4 G km−1. Only pixels where Bh(log τc = −1.5) > 300 G are considered (see text for details).

In the text