© ESO 2020

1. Introduction

The concept of plage in solar physics was coined more than 100 years ago from the analysis of low spatial-resolution Hα observations. These regions showed enhanced chromospheric line emission in many diagnostics such as the Hα line (Deslandres 1893). Our understanding of the physical processes that are at work in plage is still not complete. For instance, the origin of the large values of microturbulence that are required to model chromospheric observations remains unclear (Shine & Linsky 1974; Carlsson et al. 2015).

In the photosphere, reconstructions of the magnetic field vector indicate a strongly vertical orientation. On the other hand, the balance between magnetic pressure and gas pressure in the chromosphere is dominated by the former and the magnetic field is expected to expand horizontally, forming a canopy above the photosphere where the magnetic concentrations are confined in the intergranular lanes. The chromospheric canopy is hot and it leaves strong imprints in many chromospheric lines (Shine & Linsky 1974; Solanki et al. 1991; de la Cruz Rodríguez et al. 2013a; Carlsson et al. 2015). The first direct observational evidence of the expansion of the magnetic field was reported by Sanchez Almeida & Martinez Pillet (1994). More recent studies with high spatial resolution observations have reconstructed this expansion in the photosphere, where the reconstructed magnetic field concentrations become slightly more space-filling in the upper photosphere (Buehler et al. 2015).

The chromospheres of plage have been of particular interest for researchers because strong radiative cooling is predicted through the Ca II H&K lines, the Ca II infrared triplet lines, and the Mg II h&k lines (Vernazza et al. 1981; Carlsson & Leenaarts 2012) and that energy must be provided efficiently by an energy transport mechanism. Physical processes related to the presence of strong magnetic fields are obvious candidates to deliver this energy (see, e.g., Hasan & van Ballegooijen 2008; Judge et al. 2010; van Ballegooijen et al. 2011; Priest et al. 2018). Therefore, understanding the magnetic field topology and its connectivity between the photosphere and the chromosphere could also help to understand the energy balance in these regions. Interestingly, the gradient of the magnetic field is large at the lower edge of the magnetic canopy and, therefore, a large electric current density is predicted by Ohm’s law (j = ∇ × B/μ₀). The presence of those electric currents can lead to energy dissipation at the lower boundary of the chromosphere. The methods presented in this study can greatly contribute to understanding the energetic contribution of these currents to the chromospheric energy budget, which is our ultimate goal.

Although inversion methods have allowed retrieving plage magnetic fields in the photosphere in great detail, in the chromosphere the situation has been different. The limited selection of spectral lines with sufficient magnetic sensitivity combined with complex non local/non equilibrium physical processes that must be accounted for in order to model most of these lines have hindered enormously this line of research (see, e.g., de la Cruz Rodríguez & van Noort 2017; Carlsson et al. 2019, and references therein). Inversion codes that allow including non local thermodynamical equilibrium effects (NLTE) exist in different implementations but in all cases they require a large amount of computational resources in order to process a large field-of-view (Socas-Navarro et al. 2015; Milić & van Noort 2018; de la Cruz Rodríguez et al. 2019). Another limitation of these codes is that they assume a smooth depth-stratified atmosphere that might not always be able to describe complex stratifications in all physical parameters, and that can influence their ability to properly reconstruct the magnetic field.

The weak-field approximation (WFA; Landi Degl’Innocenti & Landi Degl’Innocenti 1973) provides an alternative method to quickly obtain magnetic field estimates whenever certain assumptions are satisfied (see Sect. 2). The WFA does not require knowledge of atom level populations to derive the emerging intensities in Stokes Q, U and V. This method has been used extensively to analyze observations in the Ca II 8542 Å and the Mg Ib₁ lines (e.g., Pietarila et al. 2007; Martínez González & Bellot Rubio 2009; Kleint 2017; Robustini et al. 2018; Centeno 2018). However, the WFA is equally affected by low signal-to-noise in observations as inversions are.

A new generation of inversion codes have made use of the idea of spatial coherency to improve the fidelity of their algorithms. A study by van Noort (2012) proposed to couple the solution of neighboring pixels using the telescope point-spread-function. A different approach has also been suggested by Asensio Ramos & de la Cruz Rodríguez (2015), using transformations to a different space where their model parameters are sparse and a decreased number of free parameters could be used to fit their data. In a recent paper, de la Cruz Rodríguez (2019) studied the effect of Tikhonov regularization in the inversion of spatially resolved spectra using a non linear Milne-Eddington model atmosphere in combination with a Levenberg-Marquardt algorithm. They showed that the quality of the inferred model parameters is greatly improved when spatial constraints are present in the fitting algorithm.

In this paper, we propose to use a spatially-regularized WFA algorithm, based on the imposition of Tikhonov regularization (Tikhonov & Arsenin 1977) to infer the magnetic field vector from spectropolarimetric (chromospheric) observations. The underlying idea is to find magnetic field values that are not only compatible with the spectrum of the local pixel, but that also keep some spatial coherence with the nearest neighboring pixels. The latter is achieved by penalizing strong horizontal gradients in the reconstructed magnetic field. Random noise solution with large amplitudes are therefore not compatible with these spatial constraints. Here we show that this modification of the WFA least-squares equations allows improving the fidelity of the reconstruction, while also extending the usability of this technique to data with a lower signal-to-noise ratio.

We have used this new method with high spatial resolution observations acquired in the Mg I 5173 Å, the Na I 5896 Å and the Ca II 8542 Å lines. The combination of these lines allows sampling the solar atmosphere in the photosphere, the upper photosphere and the lower chromosphere with some redundancy and depth overlap. We study the expansion of the canopy as a function of height and reconstruct a vertical slice of the line-of-sight magnetic field stratification.

2. Spatially-regularized weak-field approximation

The weak-field approximation is an approximate solution of the radiative transfer equation that allows predicting Stokes Q, U and V as a function of Stokes I and its derivatives as a function of wavelength (Landi Degl’Innocenti & Landi Degl’Innocenti 1973). The assumption underlying the WFA is that the magnetic field does not change with depth and that the splitting induced by the Zeeman effect is significantly smaller than the Doppler width of the line. In strong chromospheric lines, like the Ca II 8542 Å line, not having a constant magnetic field with depth would mean that the wings, sensitive to the photosphere, would be affected by a different value of the magnetic field than the core of the line, which is sensitive to the chromosphere. The latter condition, (Δλ_B ≪ Δλ_D) defines very well the regime of applicability of the approximation (Asensio Ramos 2011):

$$\begin{array}{c}\hfill B<\frac{4\pi m_{\mathrm{e}}c}{\overline{g}{\lambda }_0{e}_0}\sqrt{\frac{2kT}{M}+v_{\mathrm{mic}}^2},\end{array}$$(1)

where $\overline{g}$ is the effective Landé factor (that in the case of 8542 Å is equal to 1.10), λ₀ is the central wavelength of the line, m_e and e₀ are respectively the electron mass and charge, c is the speed of light, k is the Boltzmann constant, T is the temperature, M is the mass of the atomic element corresponding to the line and v_mic the microturbulent velocity.

In strong magnetic concentrations like sunspots, photospheric lines frequently exceed this limit (e.g., van Noort et al. 2013), but in the chromosphere the magnetic field is typically weaker (e.g., Joshi et al. 2016). In addition, the spectral lines that sample chromospheric conditions normally have a larger Doppler width than photospheric ones as well as relatively low Landé factors. So, the WFA is well-suited for chromospheric lines in the Zeeman regime, for example Ca II 8542 Å. In fact, the WFA has been employed regularly to extract the magnetic field from observations in this line (e.g., Pietarila et al. 2007; de la Cruz Rodríguez et al. 2013b; Centeno 2018; Kleint 2017; Kuridze et al. 2019).

For Stokes Q and U an additional constraint is imposed and the approximation only works if there are no (strong) velocity gradients in the atmosphere, which can lead to opacity jumps that break the validity of the approximation.

In this study we have followed the derivation presented in Landi Degl’Innocenti & Landolfi (2004) to obtain the differential equations that associate the different Stokes parameters (Q, U, V) to the derivative of Stokes I with respect to wavelength. In these equations, the two components of the magnetic field appear: the line-of-sight component connects Stokes V to the derivative of Stokes I (Sect. 2.1), while the horizontal one connects Stokes Q and U to the derivative of Stokes I (Sect. 2.2). We have also calculated the azimuth of B in the plane of the sky (Sect. 2.3). A similar derivation, along with uncertainties estimates, is provided in Martínez González et al. (2012).

Here, we have derived our mathematical expressions imposing spatial constraints from the nearest-neighbor pixels by minimizing the difference between their values with a certain weight. The latter is done in the form of spatial Tikhonov regularization (Tikhonov & Arsenin 1977).

In principle, given that our observations are critically sampled spatially and the magnetic field is expected to be relatively smooth in the chromosphere, the inferred magnetic field cannot be extremely different between pixels that are close to each other. In particular we have considered the four nearest neighbors: the pixel on the right, the left, above and below. The value of the magnetic field (longitudinal and transverse) in a certain pixel is then also determined by the solution of the surrounding pixels and it cannot be computed individually anymore, like in the case of the standard WFA. As a result, the problem becomes global and the magnetic field must be solved for the entire field of view simultaneously. Taking into account the values of nearby pixels is expected to be a decisive improvement of the fidelity of the WFA solution as it should limit the impact of noise.

2.1. Line-of-sight component

In this section, we focus on the line of sight, or longitudinal, component of B. The differential equation of the first order gives the relation between Stokes V and the derivative of Stokes I:

$$\begin{array}{c}\hfill V(\lambda )=-\mathrm{\Delta }{\lambda }_B\overline{g}cos\theta \frac{\mathrm{d}I}{\mathrm{d}\lambda }\end{array}$$(2)

where Δλ_B is given by:

$$\begin{array}{c}\hfill \mathrm{\Delta }{\lambda }_B=4.6686\times {10}^{-13}{\lambda }_0^{2}B\end{array}$$

Equation (2) can be written in a more compact way:

$$\begin{array}{c}\hfill V(\lambda )=-C_1{B}_{\Vert }\frac{\mathrm{d}I}{\mathrm{d}\lambda }\end{array}$$(3)

where the constant C₁ includes all the constant terms that appear in Eq. (2) and B_|| is the line-of-sight component of B in which is included cos θ, with θ the angle between the observer’s line-of-sight and the normal to the solar surface.

A priori, the standard WFA treats each pixel in the image individually, ignoring the rich and redundant information that is present in the spatial dimensions.

In the standard case, the WFA is applied to observations by defining a merit penalty function and trying to minimize it, for example, for longitudinal component of the magnetic field:

$$\begin{array}{cc}& {\chi }^2=\frac{1}{N}{\displaystyle \sum _{i={\lambda }_0}^{{\lambda }_n}{\left(\frac{V_{\mathrm{obs}}^i-V_{\mathrm{synt}}^i}{{\sigma }_i}\right)}^2}\hfill \\ \hfill & =\frac{1}{N}{\displaystyle \sum _{i={\lambda }_0}^{{\lambda }_n}\frac{1}{{\sigma }_i^2}{(V_{\mathrm{obs}}^i-(-C_1\frac{\mathrm{d}{I}_i}{\mathrm{d}\lambda }{B}_{\Vert }))}^2.}\hfill \end{array}$$(4)

Hereafter, we assume that each data point is divided by the corresponding $\sqrt{N}{\sigma }_i$, dropping those terms from the notation. Since both Stokes V and I are observed quantities, we can derive the longitudinal component of the magnetic field by simply minimizing the value of χ², which in this case is trivial because the dependence on the model parameters is linear:

$$\begin{array}{c}\hfill B_{\Vert }=\frac{{\sum }_i{V}_i\frac{\mathrm{d}{I}_i}{\mathrm{d}\lambda }}{C_1{\sum }_i(\frac{\mathrm{d}{I}_i}{\mathrm{d}\lambda }{)}^2}.\end{array}$$(5)

The least-squares formula in Eq. (5) has been used widely in solar and stellar applications in the past.

With the following definition the solution is forced to be compatible also with the surrounding pixels, imposing spatial constraints from the nearby pixels by reducing the difference with a certain weight, α:

$$\begin{array}{cc}& {\chi }^2={\displaystyle \sum _{i={\lambda }_0}^{{\lambda }_n}(V_i-C_1\frac{\mathrm{d}{I}_i}{\mathrm{d}\lambda }{B}_{\Vert }^{(x,y)}{)}^2+\alpha [{(B_{\Vert }^{(x,y)}-B_{\Vert }^{(x,y-1)})}^2}\hfill \\ \hfill & +{(B_{\Vert }^{(x,y)}-B_{\Vert }^{(x,y+1)})}^2+{(B_{\Vert }^{(x,y)}-B_{\Vert }^{(x-1,y)})}^2\hfill \\ \hfill & +{(B_{\Vert }^{(x,y)}-B_{\Vert }^{(x+1,y)})}^2].\hfill \end{array}$$(6)

Here we have added two superscripts to B_|| that correspond to its (x,y) location in the field of view.

Taking the derivative of the previous formula with respect to the line-of-sight components of the magnetic field in a specific pixel yields:

$$\begin{array}{cc}& [C_1^2{\displaystyle \sum _i(\frac{\mathrm{d}{I}_i}{\mathrm{d}\lambda }{)}^2+4\alpha ]B_{\Vert }^{(x,y)}-\alpha B_{\Vert }^{(x,y-1)}-\alpha B_{\Vert }^{(x,y+1)}-\alpha B_{\Vert }^{(x-1,y)}}\hfill \\ \hfill & -\alpha B_{\Vert }^{(x-1,y)}=C_1{\displaystyle \sum _i{V}_i\frac{\mathrm{d}{I}_i}{\mathrm{d}\lambda }.}\hfill \end{array}$$(7)

This equation refers to the pixel (x, y), so, if the image has dimensions (n_x, n_y), we obtain (n_x ⋅ n_y) equations. The factor 4, in this case, that multiplies α with the number of subtracted terms in the left part of the equation depends on the position of the pixel in the image. Pixels in the corners and on the edges of the field-of-view (FOV) have a different factor and number of terms. It is also clear that every equation is bound to the others through the values of the neighboring pixels, so it is not possible anymore to obtain the longitudinal component of the magnetic field for every pixel separately. Instead, to get the solution of the problem, we have to solve a sparse linear system of equations of the form Ax = b.

The sparse matrix A contains all the terms that multiply the values of the magnetic field: the diagonal term that depends on the derivative of Stokes I with respect to the wavelength plus a contribution of the local magnetic field to the regularization term (4α). Additionally we get four non diagonal terms that are proportional to −α and that account for the dependence on the solution for neighboring pixels. The matrix A has dimensions (n_x ⋅ n_y, n_x ⋅ n_y). Each row of the matrix maps the solution of our model into one observed pixel. In each row, only five elements (or fewer) are non zero, so we can store this matrix in a compact way, discarding all zeros. Figure 1 visualizes how the regularization terms in the A matrix are structured for an example field-of-view of (6 × 6) pixels: the dimensions are (36, 36) and the values on the diagonal change according to the number of the other non zero elements. The matrix x includes all the values of the line-of-sight component of the magnetic field, which is the quantity we want to derive. The vector b contains all the right-hand side terms in Eq. (7).

Fig. 1. Regularization terms of the left-hand-side matrix for a FOV of 6  ×  6 pixels. Each row of the matrix represents how one pixel is coupled to its neighbors. The image is scaled between ±4α. Red elements are positive and blue elements are negative.

We use the biconjugate gradient stabilized method (BiCGSTAB; van der Vorst 1992) to solve the sparse linear system of equations. This is an iterative method to solve symmetric and non symmetric linear systems and it provides very good convergence properties among other methods used to solve these types of linear systems (Saad 2003).

2.2. Horizontal component

The perturbation analysis of the polarized radiative transfer equation presented by Landi Degl’Innocenti & Landi Degl’Innocenti (1973) allows obtaining a second-order solution that expresses Stokes Q as a function of the derivative of Stokes I. The dependence is different in the wings (λ ≫ Δλ_D) and in the core (λ ≪ Δλ_D):

$$\begin{array}{cc}\hfill Q({\lambda }_0)& =-\frac{1}{4}\mathrm{\Delta }{\lambda }_B^2\overline{G}{sin}^2\theta cos2{\varphi }_B(\frac{{\mathrm{d}}^{2}I}{\mathrm{d}{\lambda }^2})\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill Q({\lambda }_w)& =\frac{3}{4}\mathrm{\Delta }{\lambda }_B^2\overline{G}{sin}^2\theta cos2{\varphi }_B\frac{1}{{\lambda }_w-{\lambda }_0}(\frac{\mathrm{d}I}{\mathrm{d}\lambda })\hfill \end{array}$$(9)

where $\overline{G}$ evaluates the sensitivity of linear polarization to B, which just depends on the type of the transition. Equations (8) and (9) are used in the center and in the wing of the line, respectively. In this study we only use the line wing Eq. (9), since the applicability of Eq. (8) is restricted to the line core, while common observational samplings are likely to cover more wavelength points outside Δλ_D. In both equations there is a term depending on ϕ_B, which is the azimuth angle of B with respect to a reference direction (Martínez González et al. 2012).

Similarly, the line-wing equation for Stokes U can be defined as:

$$\begin{array}{c}\hfill U({\lambda }_w)=\frac{3}{4}\mathrm{\Delta }{\lambda }_B^2\overline{G}{sin}^2\theta sin2{\varphi }_B\frac{1}{{\lambda }_w-{\lambda }_0}(\frac{\mathrm{d}I}{\mathrm{d}\lambda }).\end{array}$$(10)

To simplify we can define two different components of the horizontal magnetic field, one with the contribution of Stokes Q and the other one with the contribution of Stokes U:

$$\begin{array}{cc}& B_{{\perp }_Q}=Bsin\theta cos{\varphi }_B\hfill \end{array}$$(11)

$$\begin{array}{cc}& B_{{\perp }_U}=Bsin\theta sin{\varphi }_B\hfill \end{array}$$(12)

As in Sect. 2.1, we obtain the traditional WFA for the transverse magnetic field component by defining a merit function and minimizing it. In this case, we have obtained it for both contributions of Q and U to B_⊥:

$$\begin{array}{cc}\hfill {\left(B_{{\perp }_Q}^2\right)}^2& ={\left(\frac{{\sum }_i{Q}_i(\frac{\mathrm{d}{\overline{I}}_i}{\mathrm{d}\lambda })}{C_2{\sum }_i(\frac{\mathrm{d}{\overline{I}}_i}{\mathrm{d}\lambda }{)}^2}\right)}^2,\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill {\left(B_{{\perp }_U}^2\right)}^2& ={\left(\frac{{\sum }_i{U}_i(\frac{\mathrm{d}{\overline{I}}_i}{\mathrm{d}\lambda })}{C_2{\sum }_i(\frac{\mathrm{d}{\overline{I}}_i}{\mathrm{d}\lambda }{)}^2}\right)}^2,\hfill \end{array}$$(14)

where C₂ includes all the constant terms in Eqs. (9) and (10) and it is different from C₁ due to the different Landé factors for Stokes Q and U. On the other hand, the term $(\frac{1}{{\lambda }_w-{\lambda }_0})$ is now included in the derivatives of the intensity and we have changed the notation to $\overline{I}$ in order to reflect this contribution of the wavelength term. The total horizontal component of the magnetic field then follows as:

$$\begin{array}{c}\hfill B_{\perp }^2=\sqrt{{\left(B_{{\perp }_Q}^2\right)}^2+{\left(B_{{\perp }_U}^2\right)}^2}\end{array}$$(15)

(see also, e.g., Martínez González et al. 2012). In contrast, many studies adopt an alternative formulation that makes use of $L=\sqrt{Q^2+U^2}$ instead of using separate estimates of B_⊥ from Q and U (e.g., Asensio Ramos 2011; Jennerholm Hammar 2014; Kuridze et al. 2018; Centeno 2018). When faced with noisy data this leads to issues because it always predicts a non zero value of L even if the Q and U observations are pure noise. As a result, when applied to very noisy data, the reconstructed |B_⊥| maps predict strong magnetic fields everywhere. In our case, Eqs. (13) and (14) keep the sign of Q and U in the numerator and, in noisy cases, the summation should average to a number close to zero if Q and U are pure noise and the denominator is non zero.

We can define a modified merit function for the transverse field that accounts for regularization in exactly the same way as for the parallel component of the magnetic field in Sect. 2.1, adding a term dependent on the regularization parameter α. Like in the case of B_∥, α is a constant that has to be calibrated according to the observational data that are being processed. Since there are two different contributions to B_⊥, the solutions for Stokes Q and U need to be calculated separately with two different linear systems of equations: A_Qx_Q = b_Q and A_Ux_U = b_U. Since Stokes Q and U are linearly dependent on $B_{\perp }^2$, we make a variable change ${\overline{B}}_{\perp }=B_{\perp }^2$ to simplify the notation. The equation per pixel for the Stokes Q component of ${\overline{B}}_{\perp }$ is:

$$\begin{array}{cc}& [C_2^2{\displaystyle \sum _i(\frac{\mathrm{d}{\overline{I}}_i}{\mathrm{d}\lambda }{)}^2+4\alpha ]{\left({\overline{B}}_{{\perp }_Q}^{(x,y)}\right)}^2-\alpha {\left({\overline{B}}_{{\perp }_Q}^{(x,y-1)}\right)}^2-\alpha {\left({\overline{B}}_{{\perp }_Q}^{(x,y+1)}\right)}^2}\hfill \\ \hfill & -\alpha {\left({\overline{B}}_{{\perp }_Q}^{(x-1,y)}\right)}^2-\alpha {\left({\overline{B}}_{{\perp }_Q}^{(x-1,y)}\right)}^2=C_2\sum _i{Q}_i\frac{\mathrm{d}{\overline{I}}_i}{\mathrm{d}\lambda },\hfill \end{array}$$(16)

while the one for the Stokes U component of ${\overline{B}}_{\perp }$ is:

$$\begin{array}{cc}& [C_2^2{\displaystyle \sum _i(\frac{\mathrm{d}{\overline{I}}_i}{\mathrm{d}\lambda }{)}^2+4\alpha ]{\left({\overline{B}}_{{\perp }_U}^{(x,y)}\right)}^2-\alpha {\left({\overline{B}}_{{\perp }_U}^{(x,y-1)}\right)}^2-\alpha {\left({\overline{B}}_{{\perp }_U}^{(x,y+1)}\right)}^2}\hfill \\ \hfill & -\alpha {\left({\overline{B}}_{{\perp }_U}^{(x-1,y)}\right)}^2-\alpha {\left({\overline{B}}_{{\perp }_U}^{(x-1,y)}\right)}^2=C_2\sum _i{U}_i\frac{\mathrm{d}{\overline{I}}_i}{\mathrm{d}\lambda }.\hfill \end{array}$$(17)

As mentioned previously the number of terms with α and the factor 4, that multiplies α in square brackets, in the left part of the equations depend on the position of the pixel in the field of view.

Each component has been calculated separately using the BiCGSTAB method, and combined afterwards using Eq. (15).

2.3. Azimuth

The azimuth has already been introduced in Eqs. (8)–(10). By the definition of Q and U, dividing Eq. (10) by (9), we obtain:

$$\begin{array}{c}\hfill tan2{\varphi }_B=\frac{U(\lambda )}{Q(\lambda )}\end{array}$$(18)

By replacing the previously derived formulas, we can obtain the standard WFA estimate of the azimuthal angle of the magnetic field:

$$\begin{array}{c}\hfill {\varphi }_B=\frac{1}{2}arctan\frac{{\sum }_i{U}_i(\frac{\mathrm{d}{I}_i}{\mathrm{d}\lambda })}{{\sum }_i(\frac{\mathrm{d}{I}_i}{\mathrm{d}\lambda }{)}^2}\frac{{\sum }_i(\frac{\mathrm{d}{I}_i}{\mathrm{d}\lambda }{)}^2}{{\sum }_i{Q}_i(\frac{\mathrm{d}{I}_i}{\mathrm{d}\lambda })}=\frac{1}{2}arctan\frac{{\sum }_i{U}_i(\frac{\mathrm{d}{I}_i}{\mathrm{d}\lambda })}{{\sum }_i{Q}_i(\frac{\mathrm{d}{I}_i}{\mathrm{d}\lambda })}.\end{array}$$(19)

Once we introduce the α parameter and the magnetic field of the nearby pixels it is no longer possible to simplify the term $\sum _i{\left(\frac{\mathrm{d}{I}_i}{\mathrm{d}\lambda }\right)}^2$ of Eq. (19), because it is precisely through that term that we impose the spatial constraints. Instead, we have to compute ${\overline{B}}_{{\perp }_Q}$ and ${\overline{B}}_{{\perp }_U}$ individually using Eq. (16) and (17) and then calculate the azimuth directly:

$$\begin{array}{c}\hfill {\varphi }_B=\frac{1}{2}arctan\frac{{\overline{B}}_{{\perp }_U}}{{\overline{B}}_{{\perp }_Q}}.\end{array}$$(20)

2.4. Scaling of the problem

The terms in Eq. (6) mix quantities of a very different nature. The leftmost term is the usual definition of χ² and if the noise is correctly estimated and the model can realistically represent the observed data, then that term should be of the order of one. The rightmost term depends on a scaling factor (α) and the square of differences between values of the magnetic field, in Gauss units in our case. In order to make that term also close to one, we can add a constant that scales those terms so that the term itself is close to one. Our proposal is to use:

$$\begin{array}{c}\hfill {\alpha }^{\prime }=\frac{\alpha }{4B_{\mathrm{norm}}^2},\end{array}$$(21)

where B_norm is a typical norm for those magnetic field differences (≈100 G). If the problem is scaled this way, both terms should be close to unity and the optimal value of α can be usually found in the range (0.1 − 10), depending on the target and how good the estimate of the norm is. For B_⊥, the norm must be squared twice as the WFA is linear with the square of the magnetic field.

If the noise estimate is not included in the WFA least-squares problem, the two terms are unbalanced, usually requiring very small values of α. In the following section, we do not include this scaling (so we assume σ = 1 and α′=α) in order to be able to compare absolute numbers in our results, but with real data this scaling of the problem makes the search for the optimal α parameter trivial.

3. Tests with an 3D rMHD simulation

In order to assess the performance of the spatially-regularized WFA, we have performed a simulated observation in the Ca II 8542 Å line using a realistic 3D magnetohydrodynamics (MHD) simulation. We have taken one snapshot from a publicly available enhanced network simulation (Carlsson et al. 2016) performed with the BIFROST code (Gudiksen et al. 2011). Snapshots from this simulation have been extensively used in previous studies and have been the workhorse of polarized line formation studies in recent years using the same spectral line (e.g., Štěpán & Trujillo Bueno 2016; Quintero Noda et al. 2017; Jurčák et al. 2018; de la Cruz Rodríguez et al. 2019).

Since the height sensitivity of the 8542 Å line is target-dependent and corrugated, the magnetic field inferred using the WFA might not correspond to a unique geometrical height in the simulation in all pixels, making the comparison with the input magnetic field vector non trivial. Since our goal is not to assess the effect of magnetic field gradients on the WFA, but the improvement delivered by our new method when compared to the traditional WFA, we have set the vertical stratification of the magnetic field in each pixel to a constant value, which corresponds to the vertical average of the magnetic field at the same pixel from 1000 ≤ z ≤ 1500 km from the mean continuum formation layer. This situation is perhaps more idealized than what we would find in a normal observation where magnetic field gradients are expected as a function of height, however, many studies have overcome the effects derived from opacity jumps by restricting the wavelength range that is used to calculate the WFA (e.g., de la Cruz Rodríguez et al. 2013b; del Pino Alemán et al. 2016; Kleint 2017; Esteban Pozuelo et al. 2019).

We computed our synthetic observations using a modified version of the RH 1.5D code (Pereira & Uitenbroek 2015; Uitenbroek 2001 that has been updated with a cubic DELO-Bezier solver of the polarized radiative transfer equation (de la Cruz Rodríguez & Piskunov 2013). Figure 2 presents monochromatic images in the four Stokes parameters at 80 mÅ from line center. The panels show a clearly chromospheric landscape with elongated magnetic features connecting two opposite-polarity patches.

Fig. 2. Synthetic Ca II 8542 Å Stokes components at line center from the Bifrost snapshot. Left to right: Stokes I, Q, V and U.

The effect of noise in Stokes I. Before discussing our results, we would like to note that the WFA model relies on having a good estimate of ∂I/∂λ, which is also derived from an observed quantity. If the observations are very noisy also in Stokes I, the estimate of ∂I/∂λ has a large uncertainty and that has an impact on the reconstructions, regardless of the noise that is present in Q, U and V. We have performed a simple test in which we assume a large noise level of σ = 10⁻¹ in I only. We have reconstructed B_∥ twice using the traditional WFA algorithm, once with the noise added I and a second time without adding noise. As Fig. 3 shows, the effect is clear and when the noise level in Stokes I is high, the reconstructed magnetic fields have lower amplitudes than in the noiseless case (cf. left column). The reason is probably that the correlation between V and ∂I/∂λ starts to break and the numerator in Eq. (5) becomes smaller. We have repeated the same test but this time with the same amount of noise in Stokes V and we get identical results when we use the spatially-regularized algorithm (right column in Fig. 3).

Fig. 3. Effect of enhanced noise levels in Stokes I. Top two rows: maps of B∥ assuming σI = 10−1 (first row), σI = 0 (second row), σV = 0 (left column), and σV = 10−1 (right column). Bottom row: density plots comparing the results obtained within each column. Each row of the 2D histograms has been normalized by the value of its integral.

The line-of-sight component. We have applied four different levels of noise σ to the simulated data (relative to the continuum intensity), including σ = 0. These values have been added to the data through a random Gaussian distribution. For each component of the magnetic field we have explored a range of α values that illustrate the effect of our method for each noise level.

Depending on the choice of noise level σ, the signal might get drowned completely (in particular for σ = 10⁻¹) and we need to verify that our Stokes V values are non zero. However, even for high noise levels we find non negligible signal in Stokes V at wavelengths around line center, as Fig. 4 shows for Stokes V at Δλ = 118 mÅ from line center. Here the signal is very weak, but with the spatially-regularized method we are nonetheless able to recover a non random map for B_∥.

Fig. 4. Monochromatic Ca II 8542 Å Stokes V image at Δλ = 118 mÅ when the noise level is σ = 10−1Ic. The signal is barely visible above the noise in some regions.

Figure 5 presents the retrieval of the longitudinal component of B for the entire FOV obtained by imposing spatial constraints on the WFA. The grid shows all combinations of chosen regularization parameter values α = [0, 10⁻¹⁰, 10⁻⁸, 10⁻⁷] and noise levels σ = [0, 10⁻³, 5 × 10⁻², 10⁻¹]. For comparison, in the upper row of the grid we have plotted the constant magnetic field stratification that we used to synthesize the profiles.

Fig. 5. Inferred B∥ from the spatially-regularized WFA of Bifrost synthetic data for different values of noise σ and regularization parameter α. The top row shows the same B∥ map from the MHD model for reference. The remainder of the grid presents results for increased regularization parameter (from left to right) – where the first column is the equivalent of the classical WFA – and increased noise level σ (from top to bottom).

When the noise levels are low (bottom rows), the amplitude and shape of the Stokes profiles is well-determined and very little or no regularization is required (lower-left part of the grid). As we increase the noise level, the solution becomes more affected by noise, but using spatial constraints allows diminishing the impact of noise in the maps, while keeping the amplitudes of the reconstructed magnetic field patches. In the extreme case of σ = 10⁻¹ (second row) there is a hint of the original structure of the magnetic field in the non regularized WFA but the amplitude of the noise is as large as those signals. Increasing the regularization parameter reduces the noise level and a coherent structure can be recovered. We note however that the reconstruction is not perfect and the regions with very weak magnetic field are lost as they did not leave an imprint above the noise in the Stokes profiles. This imprint does not need to be large at all, and even a very small hint of signal above the noise level seems to be enough (see Fig. 4). Overall, the optimal α parameter for each noise level seems to be located around the counter-diagonal of the grid.

Figure 6 presents 2D histograms of the reconstructed magnetic field amplitudes between the MHD simulation values and each case of the grid, allowing a quantitative comparison. Overall, these scatter plots show that when the noise is low and the α parameter is not overestimated (lower left six panels), both the traditional and the spatially-regularized reconstructions follow the one-to-one line. When the α parameter is overestimated, the spatially-constrained method does not allow for a sufficiently rapid spatial variation in the FOV and the strongest field amplitudes are suppressed.

Fig. 6. Density maps comparing the values for B∥ from the MHD simulation (x-axis) with those from the spatially-regularized WFA (y-axis). The panels are arranged in the same order as the ones in the grid of Fig. 5. Each column in each plot has been normalized by the sum of the values in that column.

When the noise level is large and comparable to the signal amplitudes, the traditional method performs poorly and only a noisy image is reconstructed (top left four panels). In this case the spatially-regularized method clearly outperforms the traditional method by using information from surrounding pixels. The scatter plots corresponding to σ = 10⁻¹ show that even for an α value that was not overestimated for lower noise levels (third column), the reconstructed magnetic field has lower amplitudes than the original. We have investigated this effect and found that it partly originates from the fact that Stokes I is very noisy and therefore its derivative is also poorly determined observationally (see Fig. 3 and discussion thereof). We have conducted a test where we did not add the noise to Stokes I and the reconstructed amplitudes were much closer to the one-to-one reconstruction. Having a noise level of σ = 10⁻¹ is very unlikely in a real observation.

The transverse component of the magnetic field. Since the amplitudes of the signals in Q and U are smaller than in V, the range of α parameter values has been adjusted accordingly when computing the transverse component of the magnetic field. Specifically, we have chosen α = [0, 10⁻⁸, 10⁻⁷, 10⁻⁶]. The noise levels in this case are σ = [0, 10⁻³, 10⁻², 5 × 10⁻²]. Using the same regularization values as for B_∥, it would not have been possible to obtain B_⊥ for the largest noise levels.

Figure 7 presents the results obtained for B_⊥. The first row corresponds to the horizontal component of magnetic field from the MHD simulation. The effect of adding spatial constraints in this case yields similar results to the reconstructions of B_∥, but with some notable differences. When the noise level is low or absent (lower-left corner of the grid) we would expect the WFA to yield an almost perfect reconstruction, however this is not the case. There are several reasons why the result is actually showing more spatial structure (artifacts) than the top row:

• The velocity field present in the simulation is affected by gradients that can make the WFA not valid for the transverse component of B, in particular when those gradients are strong.

• Equation (9) assumes that we know the wavelength offset relative to the local frame of rest of the atmosphere where the signal is originating. So even if there are no velocity gradients, a constant velocity causes artifacts.

• The WFA for Stokes Q and U is a second-order term and its applicability is limited to very weak fields, which might not always be the case in these synthetic data.

Fig. 7. Inferred B⊥ for different values of noise σ and regularization parameter α. Format as for Fig. 5.

Although all these factors might be acting simultaneously, we can try applying the WFA, carefully assessing the quality of the results. Since in a real observation we have no prior knowledge of the velocity stratification, we have not attempted correcting the wavelength grid in each pixel for the chromospheric line-of-sight mean velocity. As a result, even for low noise levels, having a larger regularization parameter helps removing those artifacts. When the noise level is very large, the traditional WFA reconstruction yields only noise (second row, leftmost panel). As the regularization parameter is increased, the algorithm is capable of lowering the noise while recovering more structure from the data.

Figure 8 presents 2D histograms of original versus reconstructed B_⊥, in similar format as Fig. 6. We find that the unconstrained results overestimate the magnetic field amplitude due to the presence of small-scale artifacts. Therefore, for all cases using a larger α than what a low noise level would otherwise require seems to produce the correct amplitudes (see the rightmost column in Fig. 8). While we unfortunately have no way to know the real solution a priori from observations, a good approach would seem to be to choose the smallest value of α that yields a smooth reconstruction at the smallest scales.

Fig. 8. Density maps comparing the values for B⊥ from the MHD simulation (x-axis) with those from the spatially-regularized WFA (y-axis). Format as for Fig. 6.

The magnetic field azimuth. The same procedure as been repeated for ϕ using the formulae obtained in Sect. 2.3. The regularization and noise values are the same as in the case of the horizontal component of B, that is to say, α = [0, 10⁻⁸, 10⁻⁷, 10⁻⁶] and σ = [0, 10⁻³, 10⁻², 5 × 10⁻²]. Figure 9 presents the reconstructed grid for the azimuth for the entire FOV. For comparison we show the azimuth from the MHD simulation in the upper row. Again we find that applying WFA without any spatial constraints does not give reliable results if the level of noise is non zero. However, if we increase the weight that the nearby pixels have, we see that the map for the azimuth becomes much more similar to the predicted ones.

Fig. 9. Inferred ϕ for different values of noise σ and regularization α. Format as for Figs. 5 and 7.

Fig. 10. Density maps comparing the values for ϕ from the MHD simulation (x-axis) with those from the spatially-regularized WFA (y-axis). Format as for Figs. 6 and 8.

In the case of the azimuth it is clear that the Stokes Q and U signal outside the central part of the FOV is very weak, but given the very smooth nature of the azimuth we can choose a relatively large α value that tightly constrains the largest scales. Even for the largest noise level that we have considered, the reconstruction of the azimuth is very consistent spatially with the MHD map, although the errors in the outer part of the FOV are noticeably larger.

With real observations we can optimize the regularization parameter α by exploring a range of α values (which then is not a grid since the noise level is fixed) and again selecting the smallest α for which the reconstruction is smooth also at small scales.

In Appendix A we discuss other benefits of imposing regularization that we found when analyzing real observations.

4. Magnetic fields in active region plage

4.1. Data acquisition and reduction

Our observations were acquired at the Swedish 1-m Solar Telescope (SST; Scharmer et al. 2003) on 2018-06-19 at 07:30:46–07:43:43 UT, targeting a plage region at coordinates (X, Y) = (229″,61″), corresponding to viewing angle μ = 0.97. The CRisp Imaging Spectro-Polarimeter (CRISP Scharmer et al. 2008) obtained full-Stokes observations in the Ca II 8542 Å, Na I 5896 Å and Mg I 5173 Å lines.

Ca II 8542 Å was critically sampled at 21 equidistant positions around the line center with a step size of Δλ = 55 mÅ. Two extra points in the outer wings are located at ±0.88 Å from the line center. The Mg I line was sampled at 13 equidistant positions ±0.04 Å around line center and the two wing points located at ±0.24 Å from the center of the line. Finally, Na I 5896 Å line was sampled at 13 positions around line center with a step size of Δλ = 60 mÅ, with an extra wing point located at −0.6 Å from line center. The CRISP data were sampled with a pixel scale of 0$\stackrel{″}{.}$059 and the total cadence is 37 s.

The data were reduced with the CRISPRED/SSTRED data pipeline (de la Cruz Rodríguez et al. 2015; Löfdahl et al. 2018; Schnerr et al. 2011). Compensation of atmospheric seeing effects was performed using the Multi-Object Multi-Frame Blind Deconvolution method (MOMFBD; van Noort et al. 2005; Löfdahl et al. 2002). All three datasets were co aligned and corrected for residual rubber-sheet distortions using local cross-correlation of small patches in the image where the Ca II dataset was used as a reference.

4.2. Overview of the observations

We analyze a single scan taken at 07:33:14 UT, which had the best seeing conditions of the time-series. Figures 11–13 show an overview of the observations in Ca II 8542 Å, Mg I 5173 Å and Na I 5896 Å, respectively. The maps are not exactly at line center but at Δλ= 55 mÅ in the Ca II line, Δλ= 40 mÅ in the Mg I line and Δλ= 60 mÅ in the Na I line, where more signal can be expected in Stokes Q, U and V than at line center.

Fig. 11. Overview of the Ca II 8542 Ådata at Δλ = 0.036 Å from line center. Top: Stokes I, Q and U. Bottom: Stokes V, wavelength-averaged linear polarization map based on Stokes Q and U, and wavelength-averaged circular polarization map based on Stokes V.

Fig. 12. Overview of the Mg I 5173 Å data at Δλ = 0.040 Å from line center. Format as for Fig. 11.

Fig. 13. Overview of the Na I 5896 Å data at Δλ = 0.040 Å from line center. Format as for Figs. 11 and 12.

In the Stokes V maps (bottom left) two regions with opposite polarities are visible in the upper-left and lower-right corners of the image. These regions are connected by elongated, fibrilar features evident in the Ca II Stokes I panel (Fig. 11, top left). The signal, especially in Stokes Q and U (top middle and right), is stronger and above the noise level in the Mg I 5173 Å and Na I 5896 Å lines, while for Ca II 8542 Å Q and U are weak and at the detection and calibration limit. Stokes V exhibits signal above the noise level for all three lines. One way of improving the signal-to-noise would be by stacking several scans through time, however, in our case we did not consider this to be a good option because of the cadence of the observations (over half a minute) and the presence of waves in the solar atmosphere.

The bottom middle and right panels in Figs. 11–13 show wavelength-averaged linear and circular polarization maps, based respectively on Stokes Q, U and V. They were computed using only the central wavelength points – five for Ca II and three for Mg I and Na I – and saturated to increase the visibility of the structures present in the images. Finally, we note that the grid in the linear polarization map of Ca II (Fig. 11) is an artifact of the mosaicing technique used in the data reduction process.

4.3. Sensitivity of the observed spectral lines

Our observations include three strong spectral lines that sample a vast range of heights of the solar photosphere and chromosphere. Due to opacity variations, the outer wings of these lines are sensitive to the photosphere, whereas the inner cores are sensitive to the chromosphere.

The Mg I 5173 Å and Na I 5896 Å lines share some similarities. Both are strong lines that sample a fuzzy boundary between the upper photosphere and the lower chromosphere. They both correspond to minority species and therefore they are strongly affected by scattering close to line center. In fact the Stokes I images in Figs. 12 and 13 are strikingly similar in appearance. The Mg I line is more sensitive than the Na I to Zeeman-induced polarization (see Landé factors in Table 1). For more information about these lines we refer the reader to the publications by Rutten et al. (2011) and Quintero Noda et al. (2018), and references therein.

The core of the Ca II 8542 Å line samples the lower/middle chromosphere, and compared to the other two lines that we have discussed, it is more strongly coupled to the local physical conditions, which makes it also a better temperature diagnostic. The Stokes I panel in Fig. 11 is dominated by chromospheric elongated feature that connect bright magnetically active patches. This line is the least sensitive to Zeeman-induced polarization of our sample. For more information about the 8542 line we refer the reader to previous studies (e.g., Socas-Navarro et al. 2000; Pietarila et al. 2007; Cauzzi et al. 2008; Leenaarts et al. 2009; Manso Sainz & Trujillo Bueno 2010; de la Cruz Rodríguez et al. 2012; Štěpán & Trujillo Bueno 2016, and references therein).

In order to extract the stratification of the magnetic field, we can use the opacity variations within these lines to our advantage. We have defined four different spectral windows that sample four different height regions distributed throughout the photosphere and lower chromosphere.

In order to quantify to what heights these four windows correspond, we have extracted the spectra from a slit indicated in Fig. 14 and performed a NLTE inversion of the three lines included in our observations using the STiC code (de la Cruz Rodríguez et al. 2019). We have used the output model from that inversion to calculate response functions (RF; see review by del Toro Iniesta & Ruiz Cobo 2016), which allow for an analysis of the sensitivity of each line to perturbations applied at each depth-point of the model. By utilizing the RFs from the model and the z-scale that is calculated during the inversion (in hydrostatic equilibrium) we can easily assign an approximate geometrical formation height to each of our spectral windows.

Fig. 14. Inferred maps of B∥ from the four spectral windows summarized in Table 1. The panels are ordered according to their average formation height, as estimated using response functions. Panel a: photospheric estimate derived from the far wings of the Na I 5896 Å line. Panel b: upper photospheric estimate derived from the inner wing of the Na I 5896 Å line. Panel c: lower chromospheric estimate derived from the core of the Mg I 5173 Å line. Panel d: chromospheric estimate derived from the core of the Ca II 8542 Å line. The yellow contours indicate the area used to compute the spatial averages that are listed in the rightmost column in Table 1.

The formation height of each spectral window has been estimated using the center of gravity of the RFs within the wavelength window that was used to calculate each panel. To derive the zero point of the geometrical scale we have calculated the mean gas pressure along the slit at the height where τ₅₀₀ = 1, including the effect of magnetic pressure, and we have shifted the z-scale of each column so that there is horizontal gas pressure equilibrium at the continuum formation layer. This approach is only an approximation but it is the best we can currently do without including the Lorentz force and magnetic pressure in the inversion process, which is not a trivial task (Puschmann et al. 2010; Pastor Yabar et al. 2019; Borrero et al. 2019).

Some properties from the four spectral regions are summarized in Table 1, where we have indicated the wavelength range along with the spectral line properties for each window. The mean formation height of each region, from the deepest one to the highest one is: (157, 474, 760, 1168) km from the surface, although we have knowledge of the individual values for each column within the slit.

The Na I and Mg I lines included in this study reach similar heights at line center. However, our coarser sampling of the Na I line makes it less optimal than the Mg I line to use the line center as a chromospheric window due to strong opacity effects as a function of wavelength. Therefore the Na I was used to estimate the magnetic field in the photosphere at two different heights. The better sampling achieved in the Mg I line makes it better suited for an estimate using spectral positions close to line center while avoiding strong opacity effects, providing an estimate at the base of the chromosphere. The window that includes the Ca II line samples the highest height range.

5. Results

5.1. Reconstruction of B_∥

We have applied the spatially-regularized WFA to the data in the four spectral regions as defined in Table 1 and derived two-dimensional maps of B_∥. The α parameter was calibrated using a grid similar to those displayed in Sect. 3. The optimal value was selected in such a way that noise was reduced compared to the unregularized WFA while not degrading the spatial structures that are present in the image.

In the discussion of our results we focus on the positive polarity patch in the lower-right corner of the FOV. In the deepest photospheric panel (Fig. 14a) the magnetic field is spatially concentrated, interspersed with apparently field-free gaps (i.e., white in the maps). The line-of-sight field strength reaches kG values in most of these concentrations with a mean value (excluding the field-free gaps) of 719 G (within the yellow contours). In panel b the mean formation height is estimated to be 200 km higher up and the magnetic field has expanded and filled up a large fraction of the field-free gaps. The average strength of the field is now 692 G. In panel c the magnetic canopy has almost completely formed, filling most of the FOV. The mean magnetic field strength has decreased to 558 G. In this window we also estimate the horizontal component, which we discuss in more detail in Sect. 5.2 below. In the final window the canopy has completely formed filling all the plage patch with an average canopy value of 444 G. Most of the footpoint magnetic elements that were still visible in panel c are now hidden below the canopy. The clearest example of this effect can be found in the negative polarity patches located around (x, y) = (10″, 43″). The overall topology is smooth over the entire FOV.

5.2. Reconstruction of B_⊥

While we have applied our method to all spectral windows, only the results from spectral window (c) could be meaningfully interpreted as the data in the other spectral windows did not have sufficient signal-to-noise over the entire plage patch. Figure 15 presents the resulting maps of |B|, |B_⊥| and θ. The |B| map confirms the results from our inference of B_∥ and the entire region seems to be covered by a magnetic canopy with a mean value of 658 G. From the maps of B_⊥ it also becomes clear that the magnetic field in the center of the patch is more vertical and becomes stronger only closer to the edges. The strongest horizontal values that we get are around 607 G and they are anchored at the edge of the patch with an approximate inclination of θ = 83° relative to the line-of-sight. The mean inclination inside the plage patch is about 25° in the observer’s reference frame.

Fig. 15. Inferred |B|, |B⊥| and θ (inclination) maps for spectral region (c) (i.e., based on Mg I 5173 Å core spectral window).

5.3. Reconstruction of the depth stratification

Using the formation heights that we calculated in Sect. 4.3 we can reconstruct the stratification of B_∥ along the slit illustrated in Fig. 14. In each pixel we have four values of the magnetic field and an associated z-value. We have interpolated B_∥ in all pixels along the slit to a equidistant z-scale and the resulting stratification is displayed in Fig. 16. On the left vertical axis we have indicated the mean formation height of each spectral window using colored ticks.

Fig. 16. Vertical reconstruction of the canopy magnetic field from the observations. The vertical cut corresponds to the blue slit indicated in panel a of Fig. 14. The colored ticks indicate the mean formation region of each spectral region indicated in Table 1: (a) blue, (b) red, (c) yellow (d) black. The z-scale has been estimated for each pixel from its atmosphere derived from an NLTE inversion assuming hydrostatic equilibrium.

Assuming that our approximations are valid, this figure illustrates that the edge of the magnetic canopy must be located between 300 km and 600 km from the continuum formation layer. Table 1 summarizes some statistics of the reconstruction of B_∥. Although the mean magnetic field value does not change significantly as a function of height from 500 km to 1000 km, the large change in the standard deviation illustrates the effect of having a more homogeneous value in the upper layers and more confined and extreme values deeper down. The mean canopy magnetic field value at z = 1000 km is B_∥ = 449 G.

6. Conclusions

We have implemented a spatially-regularized weak-field approximation by imposing Tikhonov regularization. This type of ℓ-2 regularization has been commonly used in different stellar Doppler imaging applications (Piskunov 1990; Piskunov & Kochukhov 2002; Rosén et al. 2015; Kochukhov 2017), but to our knowledge this is the first time it is used in combination with the WFA in solar applications with spatially resolved data. This method exploits the sparsity of solar data to find solutions that are spatially coherent. Our implementation of the spatially-regularized WFA is publicly available¹.

Setting spatial constraints on the WFA reduces the noise present in the reconstructed maps and it improves the fidelity of the reconstruction by coupling the solution spatially. However, the spatial resolution that this method is capable of achieving is still set by the telescope diffraction limit, by the spatial and spectral sampling of the data and by the noise level. The latter is the dominant factor in the upper-right part of our grid (Figs. 5 and 7) as the higher spatial frequencies of the signal cannot be distinguished from noise any longer and we can only aim at retrieving larger scale details. In our study the parameter dependence is linear and therefore, the regularized WFA method only required a modification of the left-hand side term in the WFA equations.

The estimation of errors for the WFA becomes less obvious because now the magnetic field solution depends on more than one pixel and it would imply propagating the errors though a system of equations. The formulas provided by Martínez González et al. (2012) allow estimating errors based on a less-constrained solution and therefore they can be considered an upper limit of the error estimate.

We have applied the method to high spatial resolution observations of plage acquired in the Mg I 5173, Na I 5896 and Ca II 8542 Å lines. The WFA assumes a constant stratification of B. To overcome this limitation we have defined four narrow spectral windows from observations in three different lines, in order to reconstruct the line-of-sight stratification from the photosphere to the chromosphere. Within each of these narrow spectral windows we have assumed that the amplitude of the Stokes parameters can be described by a constant value of the magnetic field vector. Although a careful application of the methods presented in this paper could allow for the calculation of electrical current densities though Ohm’s law, our observations did not achieve the required S/N ratio to get a meaningful stratification of the transverse magnetic field component, which is also required for the calculation of j = ∇ × B/μ₀.

We have focused our analysis on plage regions, where we have studied the expansion of the field as a function of height and the formation of magnetic canopies in the chromosphere. We found that in the photosphere the field is concentrated and there are field-free gaps between magnetic patches. These magnetic concentrations expand horizontally closer to the chromosphere forming a quasi-homogeneous canopy above the photosphere. We found that the plage canopy has a mean total magnetic field of 658 G at approximately z = 760 km above the photosphere (leftmost panel Figs. 15 and 16). In the chromosphere our results suggest that the magnetic field must be strongly vertical. This could also explain the lack of strong Stokes Q and U signal in the 8542 Å line, in contrast to the 5173 and 5896 Å lines whose cores form closer to the canopy lower boundary where the field is more horizontal (due to the horizontal expansion). In addition, these lines are more sensitive to the presence of magnetic fields in the atmosphere so both effects are likely at play.

Acknowledgments

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (SUNMAG, Grant agreement 759548). JdlCR is supported by grants from the Swedish Research Council (2015-03994) and the Swedish National Space Agency (128/15). G.V. is supported by a grant from the Swedish Civil Contingencies Agency (MSB). The Swedish 1-m Solar Telescope is operated on the island of La Palma by the Institute for Solar Physics of Stockholm University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. The Institute for Solar Physics is supported by a Grant for research infrastructures of national importance from the Swedish Research Council (registration number 2017-00625). This research has made use of NASA’s Astrophysics Data System Bibliographic Services.

References

Appendix A: Other effects of the regularization terms

The Ca II 8542 Å line exhibits peculiar profiles in plage targets (de la Cruz Rodríguez et al. 2013a). Due to the stratification of the atmosphere and due to the presence of a hot magnetic canopy, the line source function is very shallow in the chromosphere and the resulting line profile has a flat line core. When applying the standard WFA, the denominator in Eq. (4) tends to be very small in plage regions. Combined with low S/N data, that equation might combine a numerator and denominator that are both very small and dominated by noise, resulting often in predictions of strong magnetic field. By adding regularization terms, we both construct a non diagonal linear system of equations and add a term to the diagonal (the local contribution). The effect of that term is not only a coupling one, but it also stabilizes the denominator when it is very small, by avoiding a division by zero.

We have illustrated this effect in Fig. A.1, where we have calculated B_⊥ using a list of increasing values of the α parameter. In the leftmost panel α = 0 and the prediction is unrealistically large. We could not see this effect in the simulated data because the simulation does not contain these type of profiles. We note that even when α is very small and the effect of spatial coupling is barely visible, regularization seems to be enough to remove those unrealistically large predictions (α = 3 × 10⁻⁶).

Fig. A.1. Inferred values of |B⊥| for the 8542 Å line using an increasing value of he α parameter from left to right. The scaling is set between 0 and 800 G in all panels. The tickmarks are given in arcsec.

A different type of regularization, which is local, can be applied instead to overcome this particular effect in the WFA. This type of regularization is also an ℓ − 2 method, usually referred to as low-norm regularization. It operates in a similar way as the diagonal damping λ parameter that is present in the Levenberg-Marquardt algorithm and it forces the solution to prefer low values of the magnetic field that are compatible with the observations.

The resulting least-squares formula for all components can be expressed as:

$$\begin{array}{c}\hfill B_x=\frac{C{\sum }_w{X}_w(I_w^{\prime })}{C^2{\sum }_w{(I_w^{\prime })}^2+\beta },\end{array}$$(A.1)

where β is the new regularization parameter, $I_w^{\prime }$ is the corresponding derivative term for a given line position and X_w is Stokes Q, U or V. While largely overestimating β does not lead to destruction of small scale structures, it suppresses large values of B_x, so a careful calibration of the β parameter is required when using Eq. (A.1).

All Tables

All Figures

	Fig. 1. Regularization terms of the left-hand-side matrix for a FOV of 6  ×  6 pixels. Each row of the matrix represents how one pixel is coupled to its neighbors. The image is scaled between ±4α. Red elements are positive and blue elements are negative.
In the text

	Fig. 2. Synthetic Ca II 8542 Å Stokes components at line center from the Bifrost snapshot. Left to right: Stokes I, Q, V and U.
In the text

	Fig. 3. Effect of enhanced noise levels in Stokes I. Top two rows: maps of B∥ assuming σI = 10−1 (first row), σI = 0 (second row), σV = 0 (left column), and σV = 10−1 (right column). Bottom row: density plots comparing the results obtained within each column. Each row of the 2D histograms has been normalized by the value of its integral.
In the text

	Fig. 4. Monochromatic Ca II 8542 Å Stokes V image at Δλ = 118 mÅ when the noise level is σ = 10−1Ic. The signal is barely visible above the noise in some regions.
In the text

Fig. 5. Inferred B∥ from the spatially-regularized WFA of Bifrost synthetic data for different values of noise σ and regularization parameter α. The top row shows the same B∥ map from the MHD model for reference. The remainder of the grid presents results for increased regularization parameter (from left to right) – where the first column is the equivalent of the classical WFA – and increased noise level σ (from top to bottom).

In the text

	Fig. 6. Density maps comparing the values for B∥ from the MHD simulation (x-axis) with those from the spatially-regularized WFA (y-axis). The panels are arranged in the same order as the ones in the grid of Fig. 5. Each column in each plot has been normalized by the sum of the values in that column.
In the text

	Fig. 7. Inferred B⊥ for different values of noise σ and regularization parameter α. Format as for Fig. 5.
In the text

	Fig. 8. Density maps comparing the values for B⊥ from the MHD simulation (x-axis) with those from the spatially-regularized WFA (y-axis). Format as for Fig. 6.
In the text

	Fig. 9. Inferred ϕ for different values of noise σ and regularization α. Format as for Figs. 5 and 7.
In the text

	Fig. 10. Density maps comparing the values for ϕ from the MHD simulation (x-axis) with those from the spatially-regularized WFA (y-axis). Format as for Figs. 6 and 8.
In the text

	Fig. 11. Overview of the Ca II 8542 Ådata at Δλ = 0.036 Å from line center. Top: Stokes I, Q and U. Bottom: Stokes V, wavelength-averaged linear polarization map based on Stokes Q and U, and wavelength-averaged circular polarization map based on Stokes V.
In the text

	Fig. 12. Overview of the Mg I 5173 Å data at Δλ = 0.040 Å from line center. Format as for Fig. 11.
In the text

	Fig. 13. Overview of the Na I 5896 Å data at Δλ = 0.040 Å from line center. Format as for Figs. 11 and 12.
In the text

Fig. 14. Inferred maps of B∥ from the four spectral windows summarized in Table 1. The panels are ordered according to their average formation height, as estimated using response functions. Panel a: photospheric estimate derived from the far wings of the Na I 5896 Å line. Panel b: upper photospheric estimate derived from the inner wing of the Na I 5896 Å line. Panel c: lower chromospheric estimate derived from the core of the Mg I 5173 Å line. Panel d: chromospheric estimate derived from the core of the Ca II 8542 Å line. The yellow contours indicate the area used to compute the spatial averages that are listed in the rightmost column in Table 1.

In the text

	Fig. 15. Inferred |B|, |B⊥| and θ (inclination) maps for spectral region (c) (i.e., based on Mg I 5173 Å core spectral window).
In the text

Fig. 16. Vertical reconstruction of the canopy magnetic field from the observations. The vertical cut corresponds to the blue slit indicated in panel a of Fig. 14. The colored ticks indicate the mean formation region of each spectral region indicated in Table 1: (a) blue, (b) red, (c) yellow (d) black. The z-scale has been estimated for each pixel from its atmosphere derived from an NLTE inversion assuming hydrostatic equilibrium.

In the text

	Fig. A.1. Inferred values of |B⊥| for the 8542 Å line using an increasing value of he α parameter from left to right. The scaling is set between 0 and 800 G in all panels. The tickmarks are given in arcsec.
In the text