Improvement of the Helioseismic and Magnetic Imager (HMI) Vector Magnetic Field Inversion Code

A variety of magnetic structures can be observed at the solar surface, from the most impressive ones such as sunspots and filaments, to the smallest ones such as individual elements of the network or faculae. Studying the light coming from these features, we can ascertain some of their physical properties. Arguably, spectropolarimetry is the best tool to infer the magnetic, dynamic, and thermodynamic properties of the solar atmosphere where the observed spectral lines are formed. Since the early 1970s various approaches have been used to solve the radiative transfer equation to find a model that closely matches an observed quantity by optimizing a merit function. Solving the radiative transfer equation for polarized light poses strong theoretical and technical challenges leading to different degrees of approximations depending on the studied case (see del Toro Iniesta & Ruiz Cobo, 2016, for deeper insight into this topic).

One of the simplest approximations for the solution of the radiative transfer equation that keeps a rigorous enough treatment of the magnetic field in the Zeeman regime is the Milne-Eddington (ME) approach (Landolfi & Landi Degl’Innocenti, 1982). The ME model assumes that the physical parameters in the atmosphere are independent of the optical depth except for the source function, which is assumed to vary linearly. While this is an oversimplified model, it has been widely used and leads to reasonable results in regions where field gradients are not too great.

Using a simple model generally imposes some limitations. For example, in the ME approximation the substructure of the resolution element is not treated, and that may influence the inferred atmospheric parameters. As shown by Solanki (1993), there are different scenarios in which a single-component magnetic field model may be unable to properly model observed features: 1.- network or plage where magnetic and non-magnetic areas are mixed in the same resolution element; 2.- penumbral filaments where two magnetic structures with different physical properties coexist in the same pixel or areas where the magnetic field depends on height. In general, these scenarios are modeled using more than one ME component for each resolution element weighted by a factor calling filling factor. A better solution for these scenarios would be given by: 1.- two-component model with one field-free component; or 2.- two magnetic components.

This work focuses on the first case: a two-component model with one non-magnetic atmospheric component, as we aim to deal with spatially unresolved magnetic features outside of active regions. For this model, there are different scenarios depending on the magnetic field feature (weak, intermediate and strong) and the size of the structure compared to the spatial resolution (for a more detailed discussion see Solanki, 1993). An example of this can be found in Harvey et al. (1972) who solves some discrepancies found during the interpretation of magnetic structures outside of sunspots by considering unresolved magnetic elements. They used Stokes I and V parameters observed in the visible spectral range near disk center together with a simplified two-component model (one magnetic and one non-magnetic weighted by a factor) to obtain the atmospheric parameters assuming a ME atmosphere during the inversion process. They concluded that, under some circumstances, it is unrealistic to treat the unresolved magnetic region as a single magnetic field vector occupying the whole pixel.

The filling factor provides rich information about the magnetic field, and fitting it in an inversion affects the resultant field strength, field orientation, and flux distribution, i.e. the entire magnetic configuration. It has also been suggested recently that lack of determination of filling factor leads to a mismatch in the line-of-sight and transverse magnetic field components inferred by inversion processing, and that this ultimately causes the hemisphere bias found in the east-west field component in the moderate field strength regions (Liu et al., 2021). This bias, seen in the full-disk vector magnetograms from both the Helioseismic and Magnetic Imager (HMI; Scherrer et al. 2012, Schou et al. 2012) on board the Solar Dynamics Observatory (SDO; Pesnell et al. 2012) satellite and the Vector Stokes Magnetograph (VSM) on the ground-based Synoptic Optical Long-term Investigation of the Sun (SOLIS; Keller et al., 2003), has been described and discussed in Pevtsov et al. (2021). It is suggested that this systematic bias may be caused by a combination of the noise in measurements and the magnetic filling factor that is not resolved accurately. Thus it is timely to improve the Very Fast Inversion of the Stokes Vector (VFISV, Borrero et al. (2011); Centeno et al. (2014)) code used to invert HMI data, so that it can properly determine the filling factor.

In this work we investigate a more complex model for magnetic features outside of active regions that is based on VFISV and can determine filling factor from HMI observations. The goals of this study are (1) to improve the inversion approach for HMI data to derive more accurate magnetic field data products generally, and (2) to mitigate the systematic hemispheric bias in the east-west component of magnetic field in moderate-field regions.

The paper is organized as follows. Section 2 describes the old and new models used to infer the atmospheric parameters with VFISV inversion code. In section 3, we evaluate the inferred parameters with both models using a Bayesian analysis. In section 4, we show examples of mitigation of the hemisphere bias in moderate-field regions in vector magnetic field inverted using the new version of VFISV. And finally, in section 5, we discuss the results of our work and show the conclusions.

The increase in the model complexity of ${\rm\mathbf{I}}_{\rm m}^{\rm syn}$ might imply strong degeneracies among model parameters. It is also important to recognize that the two-component model (${\rm\mathbf{I}}^{\rm syn}_{\rm 2c}$) simply fits to more parameters than the one-component model (${\rm\mathbf{I}}^{\rm syn}_{\rm 1c}$), so in order to properly ascertain which model better fits the observations, we need to account for the difference in the number of degrees of freedom. Among the various approaches that might be taken to quantify the performance, we have chosen to use a Bayesian approach. This approach allows for full characterization of the posterior distribution of the various model parameters (and uncertainties) as well as their possible degeneracies, and it also allows a rigorous comparison of models with different number of free parameters. This report does not discuss the specifics of Bayesian inference, we simply summarize the most relevant features for our situation and refer the interested reader to more complete review (for instance Gregory, 2005).

According to Bayes theorem (Bayes & Price, 1763), provided a model (${\mathbf{I}}_{j}^{\rm syn}$) and an observation (${\mathbf{I}}^{\rm obs}$) the posterior distribution of the model parameters is given by:

In the Bayesian approach the aim is to determine the posterior distributions for the model parameters, so some likelihood functional and a prior must be provided. Regarding the prior, we set an informative prior in which all but one of the model parameters $\mathcal{M}$ are equally probable in a certain range of the parameter space while impossible outside it. For $\eta_{0}$ the prior is normally distributed. The specific values for each parameter prior are gathered in Table 1. This configuration is chosen following the actual implementation of VFISV.

The likelihood function is set to be uncorrelated and of Gaussian form with zero mean and an observing uncertainty of $\sigma_{2,3}=47$ kDN/s, as derived from full-disk close-to-continuum polarimetric data. In addition, Stokes I and V are allowed worse fitting by means of weights, so that they have an effective uncertainty of $\sigma_{1}=120$ and $\sigma_{4}=69$ kDN/s, respectively. These weights are chosen in order to mimic the VFISV weightings used in the standard pipeline processing. In mathematical form, this likelihood is expressed as:

Provided with functional forms for the likelihood and prior, the estimation of the posterior distribution of the model parameters ($\mathcal{M}$) as well as the evidence (${\rm E}_{j}$) is done using PyMultiNest (Buchner et al., 2014), which internally uses the Multinest library (Feroz & Hobson, 2008; Feroz et al., 2009, 2019). The sampling process is extremely computationally demanding, so we have used it over a handful of representative cases in order to evaluate the model parameter inference ability of the two-component model, some of which are shown below.

The comparison of two models is essentially an extension of the formalism detailed in Kass & Raftery (1995); Arregui et al. (2015). Accordingly, the two models can be compared by using the Bayes factor (${\rm BF}$) which, for the two models considered here, is given by:

The two-component model presented here aims at being able to deal with cases in which a strong, spatially unresolved magnetic structure is embedded in a (nearly) field-free plasma. Clear candidates for this model are plage and network regions. Such regions host small-scale $\approx 0^{\prime\prime}$.1-100km (Lagg et al., 2010) relatively strong magnetic fields where a model with two components, one magnetic ($\mathbf{\rm I}_{\rm m}^{\rm syn}$) and one field-free ($\mathbf{\rm I}_{\rm nm}^{\rm syn}$), may perform better.

An example is shown in Figure 1 in which we present the result of the Bayesian approach for a single plage pixel ($\mathbf{I}^{\rm obs}$). In the lower left panels, the one-dimensional and two-dimensional posterior distributions of the model parameters for the one-component model ($\mathbf{I}^{\rm syn}_{1c}$ in red) and the two-component model ($\mathbf{I}^{\rm syn}_{2c}$ in blue) are shown. From the marginal posterior distribution for each model (panels (a) to (j)) it is evident that both models are able to infer a unique parameter set (at least those related to the magnetic field vector), i.e. the various model-parameter one-dimensional marginal distributions are relatively narrow for most of the model parameters. The two-dimensional marginal distributions (panels (1) to (36)) show that some degeneracies exist between some of the model parameters. For instance, ${\rm S}_{0}$ and ${\rm S}_{1}$ are inverse linearly related, i.e. if ${\rm S}_{0}$ increases its value, ${\rm S}_{1}$ becomes proportionally smaller, keeping a compatible likelihood. Another example is the degeneracy between $\Delta\lambda_{\rm D}$ and $\eta_{0}$ that follows a more complex relation. These two examples are known to occur for the HMI spectropolarimetric observations (Centeno et al., 2014). Some additional degeneracies seem to exist. For instance $\eta_{0}$ and ${\rm S}_{0}$ and ${\rm S}_{1}$ (panels (21) and (27)). In particular, as $\eta_{0}$ increases, so does ${\rm S}_{0}$, while ${\rm S}_{1}$ gets smaller. Another degeneracy is present between $\Delta\lambda_{\rm D}$ and ${\rm S}_{0}$ and ${\rm S}_{1}$ (panels (20) and (26)). Despite the presence of these degeneracies, the magnetic field parameters and line-of-sight velocity seem to be fairly decoupled ($\psi$ does show a small degeneracy with ${\rm S}_{0}$ and ${\rm S}_{1}$, but its value stays confined to $<12^{\circ}$ with 99.73% certainty). To this point both models $\mathbf{I}^{\rm syn}_{1c}$ and $\mathbf{I}^{\rm syn}_{2c}$ are quite similar in their behaviors, though the actual values for the distributions vary.

However, when considering $\mathbf{I}^{\rm syn}_{2c}$ the additional parameter ($\alpha_{\rm mag}$) (panels (29) to (36) and (i)) must be considered. It is interesting to see that there is still some degree of degeneracy among $B$, $\gamma$ and $\alpha_{\rm mag}$. This is related to the fact that when the magnetic field is not in the strong-field regime, the model is not sensitive to each parameter individually, but to $\phi=\alpha_{\rm mag}\,B\,\cos{\gamma}$ (the line-of-sight component of the mean magnetic flux density, see, for instance, Asensio Ramos et al., 2007). Thus, there exists an interplay between $\alpha_{\rm mag}$, $B$, and $\gamma$ while keeping the same likelihood. However, this degeneracy remains relatively small, for instance, the magnetic field strength is within a 300 G range with 99.73% certainty instead of extending to very low values of magnetic field strengths. This result is due to the fact that there is some amount of linear polarization that poses additional constraints in the inference even when the magnetic field strength is weak (see, for instance, Pastor Yabar et al., 2020, for a similar case in which linear polarization signals allow inferring the magnetic field orientation).

The discussion so far has focused on the model parameter inference, but what we really want to know is whether one of the parametric models is preferred. To do so we follow Kass & Raftery (1995) as shown in Table 2. For this plage case, the natural logarithm of the evidence for $\mathbf{I}^{\rm syn}_{1c}$ is $-230.31\pm 0.03$ and for $\mathbf{I}^{\rm syn}_{2c}$ it is $-197.42\pm 0.03$. Thus, the evidence ($2\log{BF^{2c,1c}}=65.78\pm 0.08$) strongly favors $\mathbf{I}^{\rm syn}_{2c}$ in the plage, i.e. the result for the two-component model is to be preferred. We also note that the orientation and strength of the magnetic field vector are different for the two models considered here.

For comparison, Figure 2 shows the same model comparison for a nearby quiet-Sun pixel. Similarly to the plage case, the non-magnetic model parameters other than ${\rm V}_{\rm los}$ (${\rm S}_{0}$, ${\rm S}_{1}$, $\eta_{0}$, $\Delta\lambda_{D}$) show degeneracies. For the quiet-Sun case, we also see clear degeneracies between $B$ and $\gamma$ for the $\mathbf{I}^{\rm syn}_{1c}$ case and between $B$, $\gamma$, and $\alpha_{\rm mag}$ for $\mathbf{I}^{\rm syn}_{2c}$ (see panels (1), (29) and (30)). This is the expected behavior in the weak field regime and in the absence of linear polarization signals, where the inferred parameter is the line-of-sight magnetic flux density $\phi={\rm B}\,\cos{\gamma}$ (or $\phi=\alpha_{\rm mag}\,{\rm B}\,\cos{\gamma}$ in the $\mathbf{I}^{\rm syn}_{2c}$ case). The azimuth parameter ($\Psi$) is hardly-defined as there is no significant linearly polarized signal. For the quiet-Sun case, we find only a slight preference for the $\mathbf{I}^{\rm syn}_{2c}$ as $2\,\log{\rm BF}^{2{\rm c},1{\rm c}}=3.85\pm 0.07$. It is important to note here that this behavior is expected for any inversion case outside the strong-field regime without linear polarization signals. In particular, for plage regions close to disk center (i.e. with low or no polarization signals at all) there is no unique determination of B and $\gamma$ (and $\alpha_{\rm mag}$ for the inversion with filling factor as a free parameter).

Figure 3 shows the line-of-sight magnetic flux density calculated as $\phi=\alpha_{\rm mag}\,{\rm B}\,\cos{\gamma}$ for the selected pixels in Figure 1 (plage, top panel) and Figure 2 (quiet-Sun, bottom panel). The magnetic flux values obtained in quiet-Sun pixels are similar, whether the Stokes parameters are inverted with one component or with two. This is because the magnetic field strength and the inclination (and the filling factor too for the two-component model) are not inferred independently, i.e. the actual inferred parameter is the mean line-of-sight magnetic flux density. In plage pixels the magnetic flux density inferred using the one-component and two-component models are different, i.e. the two-component inversion leads to a different magnetic field orientation and greater line-of-sight magnetic flux density. The value is $\approx 30\,Mx/cm^{2}$) smaller when the pixel is inverted with only one component (the median is $\approx 470\,Mx/cm^{2}$ compared to $\approx 500\,Mx/cm^{2}$).

There is a hemispheric bias in medium-strength field regions seen in vector magnetograms from both VSM/SOLIS and HMI/SDO as reported by Pevtsov et al. (2021). The bias refers to a change in the sign of the heliographic east-west component of the magnetic field when a region crosses central meridian. Figure 4 shows an example. The top panels show the three components of magnetic field on the solar disk, $B_{r}$ (radial component), $B_{t}$ (north-south component), and $B_{p}$ (east-west component), observed at 17:12 UT, November 17, 2015. The bottom panels are the vector components measured at 17:00 UT, November 20, 2015, about three days later. The data are derived by applying VFISV to Stokes [I, Q, U, V] with a constant filling factor ($\alpha_{\rm mag}=1$). The 180^∘ ambiguity of the azimuth in the inverted magnetic field is resolved using the “minimum energy” algorithm (Metcalf, 1994; Leka et al., 2009). Two fairly stable regions, A and B, are identified in the $B_{p}$ panels on the right. Looking first at the radial-component panels on the left, Region A in the north has a negative-polarity feature toward the right with a less concentrated medium-strength positive polarity pattern toward the left. Region B in the south is entirely medium-strength and has a positive radial-field component throughout. Now, consider the panels that show $B_{p}$. The east-west component of the southern medium-strength magnetic field in Region B changes sign from negative (black, upper-right panel) to positive (white, lower-right panel) after it rotates past central meridian. Similarly, the medium-strength portion of Region A in the northern hemisphere changes sign from negative to positive. The sign of the most leading negative feature changes from positive to negative, as well.

One possible cause for the presence of this hemispheric bias is the fact that a pixel might not be fully occupied by a magnetic atmosphere, i.e. typical magnetic structures outside of sunspots are smaller than the area of an HMI pixel. If this is the case, when performing an inversion assuming that the whole pixel is magnetic (as is the case for current HMI data products that depend on the ME inversion) might lead to an improper magnetic field vector inference, which might be the cause of the sign reverse in $B_{p}$. This limitation can be alleviated by allowing the inversion to also invert a filling factor that specifies the fraction of the pixel that is to be considered magnetic. It is important to note here that the inclusion of the filling factor, despite being linear, has a different relation to the longitudinal and transverse magnetic field components. For instance, let us take a simple model of the weak field approximation (Jefferies & Mickey, 1991). Measurement of the line-of-sight field is proportional to the filling factor ($\alpha_{\rm mag}$), while the transverse field is proportional to $\sqrt{\alpha_{\rm mag}}$. Thus, when $\alpha_{\rm mag}$ is not accurately determined, the various spherical components of the magnetic field vector might reflect spurious effects (such as the hemispheric bias) that are due to the fact that the magnetic field vector is not accurately inferred. It is important to note that in the original VFISV, the filling factor is set to be 1. This is close to the value found in umbrae and penumbrae, but significantly larger than the one found in plages and quiet-Sun regions. Consequently, the effect of the filling factor will mostly affect the inference of the plage - quiet-Sun regions while the penumbrae and umbrae inference will remain unchanged. While we avoid any further discussion on the specifics of the relation between the inference of the magnetic field vector and the presence of the hemispheric bias here, the reader is referred to Liu et al. (2021) for a thorough discussion on this topic.

The new version of VFISV with $\alpha_{\rm mag}$ inverted reduces this mismatch, and thus mitigates the hemisphere bias. Figure 5 shows $B_{p}$ in region A (Figure 4) observed three days apart. Left panels are from the two-component VFISV; right are from original VFISV. To avoid projection effects, $B_{p}$ is mapped onto heliographic coordinates with the X-axis in degrees and Y-axis in sin(Latitude). The region is in the east hemisphere on 17 November (top panels) and west hemisphere on 20 November (bottom panels). The sign change has been mitigated substantially.

$B_{p}$ in region B (Figure 4) is shown in Figure 6. Left panels are from the two-component VFISV and those on the right from the original. Again, the sign change in $B_{p}$ has been mitigated. More quantitatively, we calculate the percentage of pixels that have $B_{p}<0$ for this region. Only the pixels with field strength greater than 300 Mx cm^-2, $3\sigma$ of the HMI vector data (Hoeksema et al., 2014), are included. The percentage for the original data changes from $(84.8\pm 2.2)\%$ in the east hemisphere (November 17) to $(5.1\pm 1.3)\%$ in the west (November 20). In contrast, it is $(44.3\pm 3.9)\%$ in the east hemisphere and $(43.5\pm 4.0)\%$ in the west for the new data. The error refers to a confidence interval at a 99% confidence level. The hemisphere bias has been mitigated substantially when the filling factor is inverted with the two-component VFISV.

We can roughly estimate how many pixels in the full-disk vector magnetograms may have a wrong sign in $B_{p}$. We examine the 17:12 UT 2015 November 17 full disk vector magnetograms in which no active regions are on the disk. About $(22.1\pm 0.2)\%$ pixels have the opposite signs in $B_{p}$ between the original and new inverted data. The pixels included are over the entire Sun’s disk with a confidence index of disambiguation greater than 50. The confidence index denotes confidence level of disambiguation. In high-confidence pixels the minimum energy method is applied; in intermediate confidence pixels the minimum energy result is smoothed using an acute angle method (Leka et al., 2009; Hoeksema et al., 2014). Pixels with confidence greater than 50 are either high confidence or intermediate. If we assume the new version of VFISV produces true magnetic field, this implies that about 24% pixels over the disk in the original magnetograms have the wrong sign in $B_{p}$.

Vector magnetic field in the spherical coordinates, $B_{r}$, $B_{\theta}$, and $B_{\phi}$, have contributions from observed line-of-sight and transverse fields. $B_{r}$ may have this hemispheric bias problem, as well. Shown in Figure 7 are normalized histograms of $\delta B_{r}$ for regions A (left panel) and B (right) for the November 17 data (black curves) when the two regions are in the east hemisphere and November 20 data (red curves) when the regions are in the west hemisphere. Here, $\delta B_{r}=|B_{r}^{old}|-|B_{r}^{new}|$. Only the pixels with $\alpha_{\rm mag}B>300.0$ are included. For both regions, $\delta B_{r}$ becomes larger when the regions move to the west, implying this hemispheric bias problem also exists in $B_{r}$ component. Thus correction for this bias is necessary and important.

We have presented a modification of the original VFISV inversion code to invert the Stokes parameters with a two-component ME model for HMI data in order to better handle intermediate-strong small-scale (unresolved) magnetic structures embedded in non-magnetized (or weakly magnetized) atmosphere. This is motivated because the spatial resolution of the HMI/SDO instrument is not high enough to fully resolve small-scale magnetic structures in the photosphere.

The two-component model solves for an atmosphere that includes in each pixel both a magnetic component and a non-magnetic component, with the magnetic fraction (filling factor) specified by $\alpha_{\rm mag}$. In order to keep the model complexity as simple as possible, we have chosen to use the same non-magnetic model parameters ($\eta_{0}$, $a$, $\Delta\lambda_{D}$, ${\rm V}_{\rm los}$, $\rm S_{0}$, $\rm S_{1}$) for both the magnetic and non-magnetic contributions. Additionally, the magnetic parameters for the non-magnetic component are set to zero. This way only one additional model parameter is required compared to the standard VFISV model: the filling factor $\alpha_{\rm mag}$.

By means of the Bayesian approach, we have demonstrated that the application of the two-component ME model to weak- and medium-strength magnetic features provides a significantly better inference of the magnetic field vector. Of particular interest are medium field strength features where the linear polarization signals are just strong enough to break the degeneracy of the inferred magnetic field strength, inclination, and filling factor, i.e. outside of strong-field active regions and outside of the quiet-Sun). In this regime, the two-component model is strongly preferable (as deduced from the Bayes factor between both models). The new inversion leads to a different inferred magnetic field vector that is more consistent with the spectropolarimetric observations. We also find that for quiet-Sun pixels with low polarization signals, even if the two-component model is marginally preferable, the additional computational cost together with the fact that the inferred magnetic field vector components are weakly constrained (the actual constraining parameter is the magnetic flux density) discourages the usage of the two-component model. In other words, if this method is to be systematically applied, setting a limit on the total polarimetric signal for which pixels it is useful would inform a smart hybrid systematic inversion of the HMI spectropolarimetric data.

We have shown that the application of this new version of VFISV mitigates the problem of the hemisphere bias in east-west component of magnetic field in regions with intermediate magnetic fields. This is because a change in the orientation of the magnetic field vector in the line-of-sight reference frame produces a change in the orientation of the vector in the local one. This change also alters $B_{\phi}$, improving the value of this component with the strategy shown in this paper. This hemispheric bias may also affect other components, including $B_{r}$. For a nominal full-disk vector magnetogram, roughly 22% of pixels with decent polarization signal might have the wrong sign in $B_{p}$. The computing time using a lower polarization threshold (0.25 % in total polarization) increases by about 55%. Searching for the optimal threshold of polarization and speeding up the code are under investigation.