# Magnetometry: Solving the fundamental ambiguity from line pair inversion

Véronique Bommier

► To cite this version:

Véronique Bommier. Magnetometry: Solving the fundamental ambiguity from line pair inversion. Astronomy & Astrophysics - A&A, 2025, 694, pp.A40. 10.1051/0004-6361/202348632. obspm-04959074

HAL Id: obspm-04959074

https://hal-obspm.ccsd.cnrs.fr/obspm-04959074v1

Submitted on 20 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Magnetometry: Solving the fundamental ambiguity from line pair inversion

Véronique Bommier*

LIRA, Observatoire de Paris, Université PSL, Sorbonne Université, Université Paris Cité, CY Cergy Paris Université, CNRS, 5 place Jules Janssen, 92190 Meudon, France

Received 16 November 2023 / Accepted 22 November 2024

## ABSTRACT

*Context.* All magnetic field vector measurements lead to ambiguous results.

*Aims.* We demonstrate that carrying out observations in two different lines belonging to the same multiplet, but with different absorption coefficients (and thus formed at two different depths, e.g., Fe I 6302.5 Å and 6301.5 Å) allow for the azimuth ambiguity remaining from the Zeeman signal interpretation to be resolved.

*Methods.* The measurement based on an interpretation of the Zeeman effect relates to the magnetic field, *H*, and not the divergence-free magnetic induction, *B*. We analyzed how the anisotropy of the photosphere, which is strongly stratified due to gravity and density at the star surface, affects the estimation of div *H*. We also considered how the ambiguity resolution ought to be performed in these conditions.

*Results.* We obtained two ambiguity-resolved field vector maps at two different, but nonetheless close altitudes. This allowed us to derive the current density full vector via curl *H* = *J*. The process revealed the horizontal component of the current density, which is generally found to be markedly larger than the better known vertical one. We observed some systematical trends and present some examples in this paper, such as circular currents wrapping spots clockwise about a positive polarity spot and counterclockwise about a negative polarity spot. We also detail the strong horizontal current components crossing active region neutral lines.

*Conclusions.* Finally, we remark that the Na I D1 and D2 lines form another such line pair. We propose them as a way to access the low chromosphere where they are formed. However, a better spatial resolution for both observations and analysis would likely be necessary in such a medium.

**Key words.** Sun: chromosphere – Sun: magnetic fields – Sun: photosphere – sunspots

## 1. Introduction

In the solar atmosphere, magnetic field and electric currents are strongly coupled via Ampère's circuital law. Moreover, and in brief, the presence of strong electric currents (e.g., taking the form of a current sheet) could play host to a potential eruptive event. Thus, knowledge of the electric current is essential for understanding the atmosphere stability and evolution for space weather purposes. However, the electric current vector components result from spatial derivatives of the magnetic field vector components via the Maxwell equation: curl *H* = *J*. This makes them particularly difficult to derive from the observations because this requires a good level of accuracy for the primitive quantity (in the mathematical sense) comprised by the magnetic field. This also requires the magnetic field vector to be fully known, with no ambiguity. The Zeeman effect observation is ambiguous in terms of magnetic field vector because two field vectors symmetrical with respect to the line of sight (i.e., which have the same longitudinal component but opposite transverse components) are responsible for the same polarization. Thus, they are not distinguishable at the spectropolarimetric data inversion step. This ambiguity is often referred to as the azimuth ambiguity (180° ambiguity or fundamental ambiguity).

The derivation of the electric current full vector requires for the ambiguity to be solved after the spectropolarimetric

data inversion. A review of different methods can be found in Metcalf et al. (2006). However, sunspots offer the opportunity of an intuitive solution of the ambiguity in their penumbra from the spot magnetic polarity. For a negative spot, the penumbra field vectors point towards the spot center; whereas for a positive spot, the penumbra field vectors diverge from the spot center. Some authors took advantage of this opportunity to obtain current density vectors in spot penumbrae. Socas-Navarro (2005) applied the Stokes Inversion based on Response functions (SIR) code (Ruiz Cobo & del Toro Iniesta 1992) for the spectropolarimetric data inversion. This code uses several depth nodes, where the magnetic field is independently determined, which enables the derivation of the field variation along depth. The ambiguity intuitive solution was further introduced in the spot penumbra. Then, the full current density vector was derived for the spot. The same method was applied to SpectroPolarimeter of the Solar Optical Telescope of the Japanese Solar Satellite HINODE (HINODE/SOT/SP) data by Puschmann et al. (2010a). The quality of HINODE data permitted a coherence to be seen between the filamentary structure of the penumbra and the current density vector lines. Significant currents are seen to wrap around the hotter, more elevated regions with lower and more horizontal magnetic fields, which harbor strong upflows and radial outflows (the intraspines).

Indeed, the derivation of the three coordinates of the current density vector requires knowledge of the depth variations of the magnetic field, which is rarely achieved with current

<footnote>* Corresponding author; V.Bommier@obspm.fr</footnote>
methods. In general, a single line has been observed, which is
formed at a given depth. We can derive the magnetic field at
this single depth and only the vertical component Jz of the cur-
rent density can be derived, except in the case where the SIR
code or a code of similar technique is applied. As a result, full
vector current density determinations are still rare. Below, we
present such determinations, which are not based on the SIR
code application, but on two line observations and single-line
Milne-Eddington inversions. The inversion method is described
in Section 2. The depth difference in line formation altitude was
previously directly derived from other HINODE/SOT/SP data
(Faurobert et al. 2009). This enables the complete derivation of
div **H** from the observations, which allows for the ambiguity
resolution to be determined according to the method described
in Section 4. In particular, we analyze how the anisotropy of the
photosphere, which is strongly stratified due to gravity and den-
sity at the star surface, affects div **H** and how the ambiguity reso-
lution ought to be performed in these conditions. We present our
example results in Section 5. As in the study by Puschmann et al.
(2010a), we find that the horizontal component of the current
density vector is much stronger than the vertical one. The ver-
tical component is the most frequently reported and discussed.
The horizontal component remained generally unknown, which
is much larger still. In this work, we discuss examples of the
typical structures we derived for the horizontal component in
sunspot penumbrae and across some neutral lines. Obviously,
applying the SIR code inversion (Ruiz Cobo & del Toro Iniesta
1992) would provide more refined vertical variations for evalu-
ating divergences and electric currents.

We tested our multiline ambiguity solution method on differ-
ent spectral lines observed with Télescope Héliographique pour
l'Étude du Magnétisme et des Instabilités Solaires (THÉMIS)
operating in its new multislit mode. The two splitted beams enter
the spectrograph by two parallel slits. As a result, we found that
the pair of lines must belong to the same multiplet in order to
have been formed in the same manner, but with different absorp-
tion coefficients. This is exemplified by the case of the Fe I
6301.5 Å and 6302.5 Å line pair. With these conditions, their
depth formation difference remains rather constant as visible in
the simulation by Khomenko & Collados (2007). This line pair
is interestingly the one observed by HINODE/SOT/SP. In Sec-
tions 5.3 and 5.4, we present the magnetic field maps observed
in a sunspot penumbra and above a neutral line.

We ought to recall here that what the measurement based
on the interpretation of the Zeeman effect is the magnetic field,
H, and not the magnetic induction, B (Bommier 2020). The
magnetic induction, B, is divergence-free, whereas the magnetic
field, H, is not. They are related via B = μ₀(H + M), where
M is the magnetization. However, in typical solar models, the
reported electron density gives the modulus of M as very small,
with respect to the modulus of H, which results in B ≈ μ₀H; this
also makes H appear rather divergence-free. In this respect, min-
imizing div H makes sense for solving the ambiguity. However,
the Sun's surface electron density could be much higher due to
electron thermal escape in the solar interior, where the electron
thermal velocity overpasses the star gravitation escape velocity
by a factor of 14 and the proton escape velocity by a factor of
6, as remarked by Bommier (2020). The model electron density
is, in fact, derived from ionization equilibria modelled from the
spectrum; in addition, this is considered within the local electric
neutrality hypothesis (Vernazza et al. 1973, Section V), which
has to be questioned in the presence of the star gravitation field,
as explained above. Thus, div H could be non-zero as observed
(see the review by Balthasar 2018), but its minimization could

be kept as a method for solving the fundamental ambiguity. The
results presented in the present paper confirm this approach. In
practice, we apply Gauss units to describe the magnetic field,
although they are, in fact, the units for μ₀H.

## 2. First step of the data analysis: Milne-Eddington inversion

We applied the Milne-Eddington inversion code known as the
fitting of the Unno solution of the radiative transfer equation
(UNNOFIT). Based on the Unno-Rachkovsky solution of the
transfer equation for the radiation Stokes parameters of a spec-
tral line, the code accounts for all the magneto-optical effects
(Landolfi & Landi Degl'Innocenti 1982). The code applies the
Levenberg-Marquardt algorithm to the theoretical profile to
make it fit the observed profile. Eight parameters were entered in
the algorithm, namely: the line center frequency, λ₀, the Doppler
width, ΔλD, η₀, the three magnetic field coordinates, the Voigt
parameter, a, and the ratio B₁/B₀ of the two parameters charac-
terizing the Milne-Eddington atmosphere. This was the original
code developed by Landolfi et al. (1984).

However, thanks to the pioneering work by Stenflo (1973),
it is known that the atmosphere is permeated by unresolved
magnetic structures, which can be roughly modelled with a
magnetic filling factor α. Following a wish by Egidio Landi
Degl'Innocenti, we introduced such a magnetic filling factor
in his UNNOFIT code. We tested the code and validated it
(Bommier et al. 2007), but we obtained that the unresolved mag-
netic field strength, *H*, and the magnetic filling factor α cannot
be separately determined (see Fig. 4 of that paper). Here, *H* is
here the modulus of the magnetic field vector, which is different
from its longitudinal component. Only their product *αH*, which
is the local average magnetic field strength, has been obtained
as final result. This occurs when the field is not strong enough
to well separate the Zeeman components with respect to the
Doppler width (typically 1000 G at least). This is the effect of
the larger number of parameters to be determined with respect
to the number of independent parameters provided by the line
observation. A line is characterized by a position, a width, a cen-
tral depth, and three Stokes parameters, which enable the cal-
culation of the three magnetic field components. This results in
only six independent parameters. The Voigt *a* parameter charac-
terizes the far wings and plays a specific role there. Thus, with
the filling factor, the number of searched for parameters exceeds
the number of available parameters. However, as shown by the
tests, the local average magnetic field stregh *αH* can be well
determined at last. Later on, a tenth parameter was introduced in
the inversion, which is the velocity gradient along the line for-
mation depth that depends on an additional observed parameter,
namely, the line asymmetry (Molodij et al. 2011).

The inversion was performed in the four Stokes profiles
simultaneously. The same weight of unity was given to the four
Stokes parameters for the chi-square calculation because *I* plays
an important role in the determination of the magnetic filling fac-
tor *α* (usually, *I* is less weighted). It is well-known that in weak
fields the linear polarization Stokes parameters (*Q*, *U*) quadrat-
ically depend on the transverse magnetic field, whereas the cir-
cular polarization Stokes parameter, *V*, exhibits a linear depen-
dence on the longitudinal magnetic field. As the chi-square of
the four Stokes parameter are simply added in the algorithm, the
longitudinal and transverse fields are simultaneously determined
in a unique procedure. This results in quite comparable accura-
cies. We plotted the difference histograms for the cartesian com-
ponents μ₀Hₓ, μ₀Hᵧ, μ₀Hz in the line of sight reference frame
assumed to be at the disk center for the same test data as above. We obtained quite comparable widths (or an inaccuracy of just 15 G) for $\mu_0H_x$ and $\mu_0H_y$, and 10 G for $\mu_0H_z$, for $\alpha\mu_0H > 10$ G, from noised profiles at the $1.5 \times 10^{-3}$ level for the polarimetric noise, which is not far from the HINODE/SOT/SP one that we evaluated at $1.2 \times 10^{-3}$ in profile far wings ($1.0 \times 10^{-3}$ during our last THÉMIS campaigns described in this paper).

Although the inversion finally provides only the $\alpha\mu_0H$ product, it is important to take $\alpha$ into account during the inversion for a good determination of the field inclination (Leka et al. 2022). Due to the sensitivity difference between the longitudinal and transverse fields, forcing $\alpha = 1$ (i.e., ignoring $\alpha$) would lead to bad inclinations in regions where $\alpha$ is far from unity.

Landolfi et al. (1984) developed the light UNNOFIT version for the normal Zeeman triplet and the heavier UNNOFIT2 version specific for Zeeman multiplets. More details are given in Bommier (2013). Thus, we were able to treat any kind of line.

### 3. Warning about the aspect ratio applied for the ambiguity resolution

In the present paper, we present an ambiguity resolution method, where we apply a scaling based on the aspect ratio of the stratified atmosphere (see below Sect. 4.5). Our intuition is based on the non-zero value generally observed for the magnetic field divergence, with vertical magnetic field gradients on the order of 3 G/km, when the horizontal gradient is solely on the order of 0.3 G/km (see the review by Balthasar 2018). However, it is important to advertise that the determination of these gradients depends on two things: first, the in-depth analysis, which requires appropriate inversion codes such as the SIR code (Ruiz Cobo & del Toro Iniesta 1992) for the vertical gradient; second, a refined spatial resolution, for horizontal gradient. Moreover, these gradients may also depend on the studied solar structure.

In the present paper, we did not apply the SIR or similar inversion code. We applied rougher methods instead and, in addition, the spatial resolution of the observation is not the optimal one. For instance, Puschmann et al. (2010b) conducted a detailed study of a sunspot penumbra with a better spatial resolution and they obtained smaller gradient values, on the order of 0.25 G/km for the vertical gradient (see their Fig. 8). However, their horizontal gradient is about 0.05 G/km only, which remains significantly different from the vertical gradient. They applied the SIR inversion code. Buehler et al. (2015) use observation with similar spatial resolution as ours, and find then a vertical gradient similar to ours. They apply the Stokes-Profiles INversion-O-Routines (SPINOR) inversion code (Frutiger et al. 2000). They do not strictly find a magnetic flux conservation, but they can finally impose the magnetic flux conservation to their results. When done, the gradient of their inversion agrees with the thin-tube approximation (see the end of their Sect. 4.2). They interestingly distinguish between core and canopy pixels (see their Fig. 3), which have different signs for their vertical magnetic field gradients. This underlines the effect of the spatial resolution.

Below, we present a non-exhaustive list of some effects that have been neglected when considering the vertical magnetic field gradient in sunspot umbrae in our work:

- The presence of velocity-temperature or magnetic field-temperature (or continuum intensity) correlations, which could be produced by, for instance, umbral dots at a spatial resolution below that of the observations;

- Blending of the atomic spectral lines used to measure H in umbrae by molecular lines (often unidentified or with poor molecular data);

- The effect of straylight from the penumbrae, from surrounding plage and from the quiet Sun. Straylight can play a major role in the umbra by affecting spectral lines significantly (it is not sufficient to model it via a simplistic filling factor, as we do);

- Departure from a purely hydrostatic equilibrium (e.g., by including the $J \times B$ force in the force-balance), which can affect the line formation height;

- The effect of a non-horizontal and corrugated $\tau = 1$ surface, which is strongly inclined in most parts of the Sun, including in the umbra. This can falsify the horizontal gradients of the field deduced from observations with finite spatial resolution;

- Non-LTE effects in spectral lines, which can also change their formation heights.

Consequently, as long as these and other effects have not been properly tested in realistic model atmospheres and with state-of-the-art radiative transfer computations, the scaling factor that we introduce below must be considered to be a purely empirical factor. This factor will depend on the quality of the data and of the inversion as well as on the solar feature being studied.

### 4. Data analysis second step: Azimuth ambiguity resolution

We applied a modification to the ME0 code (i.e. the Minimum Energy method) developed by Metcalf, Crouch, Barnes, & Leka and now available on the web<fnref n="1" /> (Leka et al. 2009). This code applies the “Minimum Energy” Method initially described by Metcalf (1994), which consists of searching for the field vector orientation that minimizes $|\text{div } H| + \lambda|J_z|$. In this approach, $J$ is the current density vector, $J_z$ is its vertical component, and $\lambda$ is a positive weight parameter usually fixed at unity (with $\mu_0H$ expressed in Gauss). Minimizing $|\text{div } B|$ is a natural requirement imposed by Maxwell’s equations. As explained in the introduction, we extended this requirement to what is actually being measured, namely, $|\text{div } H|$. On the other hand, minimizing the current density minimizes the maximum allowed free magnetic energy (Metcalf 1994). The minimization is performed globally on the whole map by applying the “simulated annealing” algorithm. Its application to this minimization problem is described by Crouch et al. (2009). The method is then complemented by propagating the solution via the acute angle method below a certain field strength threshold (Leka et al. 2009) presently taken at 400 G. The acute angle method consists of selecting from two ambiguous solutions symmetrical with respect to the line of sight, the solution that makes an acute angle (in the transverse plane) with the vector to be compared. In the original ME0 method, one single map is used and the vertical derivatives of the magnetic field are derived from a current-free reconstruction of the magnetic field. Reconstructions based on less restrictive hypotheses were later on introduced by this method authors (Metcalf et al. 2006; Leka et al. 2009). As the field vector ambiguity concerns its transverse component in the line of sight coordinates (whereas the reconstruction is performed in the heliographic coordinates), the back and forth transformations have to be performed between these two systems of coordinates.

<footnote n="1">¹ http://www.cora.nwra.com/AMBIG/</footnote>
More precisely, div **H** is calculated in the heliographic reference frame in the original ME0 code because $\partial H_z/\partial z$ is evaluated in this frame by the reconstruction. Our approach avoids the vertical reconstruction of the field by the introduction of two maps recorded in two lines formed at two different depths. In our method we kept the minimization procedure of ME0 but we calculated div **H** in the line of sight coordinates instead, which was possible from our two-line observations. Accordingly, we modified the corresponding subroutine of ME0. We calculated the expression of div **H** adapted to the case of a line formation plane inclined with respect to the line of sight in order to develop a method able to treat maps of any location on the solar disk. Below (and in Appendix A), we describe this calculation. A different calculation, which also accounts for the case of a line of sight inclined with respect to the local vertical, was developed in Crouch et al. (2009) and Crouch (2013, 2015). We first have to introduce the various reference frames entering the calculation.

### 4.1. Reference frames

Our approach is in agreement with Gary & Hagyard (1990) in the limit of $P = 0$. In the following, we denote the line of sight (l.o.s.) reference frame as Oxyz and OXYZ is the heliographic one at the observed region location. The l.o.s. reference frame is defined by Oz, given the l.o.s. is itself oriented towards the observer, and with Oy parallel to the disk central meridian oriented towards the solar north. The heliographic reference frame has OX aligned with the local parallel solar west oriented, OY aligned with the local meridian that is oriented to the solar north, while OZ goes along the solar radius oriented from the Sun's center and directed outward. We use **R** to denote the unit vector along the OZ axis of the heliographic reference frame, which is also the local solar radius. We used indexes $(x, y, z)$ as the vector coordinates in the l.o.s. reference frame and $(X, Y, Z)$ for those in the heliographic reference frame. Gary & Hagyard (1990) denote with the upper index, $l$, the vector coordinates in the l.o.s. reference frame, and $h$ those in the heliographic reference frame. In the l.o.s. reference frame, the **R** coordinates are:

$$ \begin{cases} R_x = \sin(L - L_c) \cos b, \\ R_y = \sin b \cos b_0 - \cos(L - L_c) \cos b \sin b_0, \\ R_z = \sin b \sin b_0 + \cos(L - L_c) \cos b \cos b_0, \end{cases} \quad (1) $$

where $L$ and $b$ are, respectively, the longitude and latitude of the center O of the observed region, $L_c$ is the disk central meridian longitude, and $b_0$ is the disk center latitude. The other unit vectors of the OXYZ heliographic reference frame are denoted as **I** along OX and **K** along OY, with

$$ \begin{cases} I_x = \cos(L - L_c), \\ I_y = \sin(L - L_c) \sin b_0, \\ I_z = -\sin(L - L_c) \cos b_0, \end{cases} \quad (2) $$

and

$$ \begin{cases} K_x = -\sin(L - L_c) \sin b, \\ K_y = \cos b \cos b_0 + \cos(L - L_c) \sin b \sin b_0, \\ K_z = \cos b \sin b_0 - \cos(L - L_c) \sin b \cos b_0, \end{cases} \quad (3) $$

in the l.o.s. reference frame. These formulae are obtained by applying three rotations successively to transform the heliographic reference frame into the l.o.s. one; namely: (1) rotation of $b$ about the OX axis; (2) rotation of $-(L - L_c)$ about the new OY axis; (3) rotation of $-b_0$ about the new Ox axis. Conversely,

the coordinates of the l.o.s. basic unit vectors **i**, **k** and **l** (l.o.s. vector) in the heliographic reference frame are

$$ \begin{cases} i_x = \cos(L - L_c), \\ i_y = -\sin(L - L_c) \sin b, \\ i_z = \sin(L - L_c) \cos b, \end{cases} \quad (4) $$

and

$$ \begin{cases} k_x = \sin(L - L_c) \sin b_0, \\ k_y = \cos b \cos b_0 + \cos(L - L_c) \sin b \sin b_0, \\ k_z = \sin b \cos b_0 - \cos(L - L_c) \cos b \sin b_0, \end{cases} \quad (5) $$

and

$$ \begin{cases} \ell_x = -\sin(L - L_c) \cos b_0, \\ \ell_y = \cos b \sin b_0 - \cos(L - L_c) \sin b \cos b_0, \\ \ell_z = \sin b \sin b_0 + \cos(L - L_c) \cos b \cos b_0, \end{cases} \quad (6) $$

The transformation from the coordinates $(x, y, z)$ of a given vector in the l.o.s. reference frame into its coordinates $(X, Y, Z)$ in the heliographic reference frame, can be written as

$$ \begin{cases} X = x_i x + y_k x + z \ell_x, \\ Y = x_i y + y_k y + z \ell_y, \\ Z = x_i z + y_k z + z \ell_z, \end{cases} \quad (7) $$

and conversely

$$ \begin{cases} x = XI_x + YK_x + ZR_x, \\ y = XI_y + YK_y + ZR_y, \\ z = XI_z + YK_z + ZR_z, \end{cases} \quad (8) $$

### 4.2. Deprojection

A map is obtained from a scan of the spectrograph slit along the solar image. The map is reconstructed by positioning side by side all the slit outputs. The mapped quantities are then given on a rectangular system of pixels in the l.o.s. reference frame, the $\Delta x$ pixel size along Ox being given by the scan step size, and the $\Delta y$ pixel size along Oy being given by the camera pixel size. Thus, we generally have $\Delta x \neq \Delta y$, which results in an anamorphosis of the map. Once the ambiguity is solved and a single magnetic field vector is obtained for each pixel, its coordinates may be transformed into the heliographic reference frame. The deprojection of the pixel array is a more complicated task because the meshes are no more rectangular overall, as described, for instance, in Gary & Hagyard (1990). In the present paper and related studies, we have roughly approximated a rectangular shape for the deprojected map, however, with pixel side sizes of $\Delta X = \Delta x/I_x$ and $\Delta Y = \Delta y/K_y$. This is exact when the map center is located on the solar equator and when $b_0 = 0$. This departs from exactitude when the latitudes of the disk and/or map centers depart from 0. The cosine of the heliocentric angle $\theta$ is $\mu = R_z$.

### 4.3. The quantity to be minimized

As stated above, we applied the “simulated annealing” algorithm to globally minimize $|\text{div } \mathbf{H}| + \lambda |J_z|$. The parameter $\lambda$ was taken as unity (with $\mu_0 \mathbf{H}$ expressed in Gauss). Indeed, in his pioneering paper, Metcalf (1994) recommends to minimize the alternative quantity $|\text{div } \mathbf{H}| + \lambda ||\mathbf{J}||$ as the physically meaningful quantity for “energy minimization”, with the full current vector and not only its Z-component. We first tried this minimization because we had also the transverse current component (within the ambiguity) from our two-line observations, but it did not converge
<figure><img src="image_2.png" /><figcaption>Fig. 1. $\Delta_x H_x^{(m)}$ is the difference between the $H_x$ components measured in pixels $P(i+1, j)$ and $P(i, j)$. It can be seen that the $z$ coordinate is not the same for the two pixels for an observation out of disk center, which is in contradiction with the definition of the partial derivative along $x$.</figcaption></figure>

correctly. On the contrary, we found that the algorithm successfully worked out when we instead minimized $|\text{div } \mathbf{H}| + \lambda |J_z|$ with $\lambda = 1$. We tried various $\lambda$ values and also found $\lambda = 1$ to be the best choice, which leads to maps at most free of artificial defaults such as strong current lines delineating zones of different ambiguous solution choice by the algorithm or with minimal checkered zones, where the ambiguity resolution failed. The minimization of $|\text{div } \mathbf{H}| + \lambda |J_z|$ was later performed also by Metcalf's team (Metcalf et al. 2006).

However, the map under study is typically not located at the disk center, which results in a certain inclination of the map formation plane with respect to the l.o.s. referred to by the heliocentric angle $\theta$ and its cosine $\mu = \cos \theta$. The “map formation plane” refers to the heliographic horizontal plane, namely, the plane perpendicular to the local solar radius. The map results from the interpretation of polarization in spectral lines. As explained below in Sect. 4.4, the radiation received along the line of sight in a given spectral line comes essentially from a certain height in the solar atmosphere (i.e., the Eddington-Barbier approximation), which defines an heliographic horizontal plane that contains the location of the points where the magnetic field is measured by spectropolarimetry. This plane is not perpendicular to the line of sight, when the observed region is not located at the disk center.

However, the divergence has to be computed in a rectangular reference frame. To do so, there are two possibilities. Either the coordinates of the two ambiguous magnetic field vectors obtained in the line of sight reference frame and symmetrical with respect to the line of sight, which are of the $H^{(l)}$ type, are transformed into heliographic coordinates of the $H^{(h)}$ type of Gary & Hagyard (1990). The ambiguity is then resolved in the heliographic reference frame, where the measurements are all located in the horizontal plane perpendicular to the $OZ$ axis.

This is the method applied by Leka et al. (2009) in their ME0 code described above. As for our work, we resolved the ambiguity in the line of sight coordinates of the field as obtained from the measurements, which are of the $H^{(l)}$ type of Gary & Hagyard (1990). In this case, the measurements are located along a plane inclined with respect to the line of sight, which is the $Oz$ axis. This inclination described below and in Fig. 1 is at the origin of balancing terms in the divergence and curl expressions in terms of the line of sight or $H^{(l)}$-type magnetic field vector coordinates (derived as described in Appendix A).

In Appendix A, we derive the expressions of div $\mathbf{H}$, $J_Z$ and $J_{x,y}$ as a function of the l.o.s. reference frame magnetic field coordinates, when the fields are measured along an inclined formation plane perpendicular to the local solar radius unit vector, $\mathbf{R}$, of l.o.s. components $R_{x,y,z}$ given in Eq. (1). We obtain

$$ \text{div } \mathbf{H} = \frac{\Delta_x H_x^{(m)}}{\Delta x} + \frac{\Delta_y H_y^{(m)}}{\Delta y} + \frac{\Delta_z H_z^{(m)}}{\Delta z} \\ + \frac{R_x \Delta_z H_x^{(m)}}{R_z \Delta z} + \frac{R_y \Delta_z H_y^{(m)}}{R_z \Delta z}. \qquad (9) $$

We use $\mathbf{H}^{(m)}$ to denote the vector measured along the inclined plane, where the line is formed. $Oxyz$ is the l.o.s. reference frame and $Oy$ is solar north oriented. Also, $\Delta_x H_x^{(m)}$ and $\Delta_y H_y^{(m)}$ are the difference of the measured component, $H_x^{(m)}$ or $H_y^{(m)}$, respectively, between two neighboring pixels in $x$ or $y$ direction. Then, $\Delta x$ or $\Delta y$ denotes the distance between the neighboring pixels in the “sky plane”, which is the plane perpendicular to the l.o.s. at the location of the observed region. The magnetic field may eventually be averaged between the two lines. If this average has to be performed before the ambiguity is solved, the acute angle method is applied between the two lines to select the ambiguous solutions to average. This is a reasonable approximation, which saves computation time in the simulated annealing. This corresponds to assume that the field lines do not twist so much. As for $\Delta_z H_z^{(m)}$, it is the difference between the values obtained from the two lines with a difference in depth formation $\Delta z$ along the l.o.s. Explicit definitions of all these quantities are given in Appendix A. The way to calculate these quantities is detailed at the beginning of the Appendix. When the map center is located on the solar equator and when $b_0 = 0$, the above Eq. (9) is simplified to

$$ \text{div } \mathbf{H} = \frac{\Delta_x H_x^{(m)}}{\Delta x} + \frac{\Delta_y H_y^{(m)}}{\Delta y} + \frac{\Delta_z H_z^{(m)}}{\Delta z} \\ + \tan \theta \frac{\Delta_z H_x^{(m)}}{\Delta z}, \qquad (10) $$

where $\theta$ is the heliocentric angle (positive for the western side of the central meridian and negative for the eastern side). The second line of these equations accounts for the $z$ variation hidden in $\Delta_x H_x^{(m)}/\Delta x$ and $\Delta_y H_y^{(m)}/\Delta y$, because $H^{(m)}$ is measured along the line formation plane, which is not perpendicular to $Oz$ in the general case.

In other words, if we consider two neighboring pixels, referred to as $P(i, j)$ and $P(i+1, j)$, separated by the length $\Delta x$ in the plane perpendicular to the line of sight, and if we also consider a quantity $A$ measured via the same line at the two places (i.e., $A_{i,j}$ and $A_{i+1,j}$) the ratio $(A_{i+1,j} - A_{i,j})/\Delta x = \Delta_x A/\Delta x$ is not an approximate value of the partial derivative $\partial A/\partial x$. This is because $\Delta_x A$ also involves a variation along the line of sight $z$ together with the variation in $x$, when the observation is not
performed at the disk center; namely, when the line formation plane is not perpendicular to the line of sight. From the definition, $\partial A/\partial x$ has to be evaluated at constant y and z. In the evaluation of $\Delta_x A = A_{i+1,j} - A_{i,j}$, we have $z[P(i+1, j)] \neq z[P(i, j)]$, as it is visible in Fig. 1, when the line formation plane is inclined with respect to the line of sight, which is also the Oz axis. As a consequence, z is not kept constant in the partial derivative numerical evaluation, when the observation is done out of disk center. This was the reason to perform the Appendix A calculations.

For the numerical computation, the divergence was computed at the center of a pixel, by averaging between the variations along each side of the pixel. We used two lines formed at two different depths. This gives us

$$ \Delta_x H_x^{(m)} = (H_x^{(1)}(i+1, j) - H_x^{(1)}(i, j) \\ + H_x^{(1)}(i+1, j+1) - H_x^{(1)}(i, j+1) \\ + H_x^{(2)}(i+1, j) - H_x^{(2)}(i, j) \\ + H_x^{(2)}(i+1, j+1) - H_x^{(2)}(i, j+1))/4, \quad (11) $$

$$ \Delta_y H_y^{(m)} = (H_y^{(1)}(i, j+1) - H_y^{(1)}(i, j) \\ + H_y^{(1)}(i+1, j+1) - H_y^{(1)}(i+1, j) \\ + H_y^{(2)}(i, j+1) - H_y^{(2)}(i, j) \\ + H_y^{(2)}(i+1, j+1) - H_y^{(2)}(i+1, j))/4, \quad (12) $$

$$ \Delta_z H_z^{(m)} = (H_z^{(2)}(i, j) - H_z^{(1)}(i, j) \\ + H_z^{(2)}(i+1, j) - H_z^{(1)}(i+1, j) \\ + H_z^{(2)}(i, j+1) - H_z^{(1)}(i, j+1) \\ + H_z^{(2)}(i+1, j+1) - H_z^{(1)}(i+1, j+1))/4, \quad (13) $$

where the indexes (1) and (2) correspond to “line 1” and “line 2”.

For the current density vector component, $J_z$, along the direction perpendicular to the inclined plane, we obtain (curl of the magnetic field):

$$ J_z = R_z \left[ \frac{\Delta_x H_y^{(m)}}{\Delta x} - \frac{\Delta_y H_x^{(m)}}{\Delta y} \right] \\ + R_x \frac{\Delta_y H_z^{(m)}}{\Delta y} - R_y \frac{\Delta_x H_z^{(m)}}{\Delta x}. \quad (14) $$

It is possible to obtain analogously the two components of the current density vector in the plane perpendicular to the line of sight (see Appendix A):

$$ J_x = \left[ \frac{\Delta_y H_z^{(m)}}{\Delta y} - \frac{\Delta_z H_y^{(m)}}{\Delta z} \right] + \frac{R_y \Delta_z H_z^{(m)}}{R_z \Delta z}, \quad (15) $$

$$ J_y = \left[ \frac{\Delta_z H_x^{(m)}}{\Delta z} - \frac{\Delta_x H_z^{(m)}}{\Delta x} \right] - \frac{R_x \Delta_z H_z^{(m)}}{R_z \Delta z}. \quad (16) $$

From these relations and by applying Eq. (7), the z component of the current density vector in the line of sight reference frame can be derived, which is

$$ J_z = (J_z - J_x R_x - J_y R_y) / R_z, \quad (17) $$

which can be reduced into

$$ J_z = \left[ \frac{\Delta_x H_y^{(m)}}{\Delta x} - \frac{\Delta_y H_x^{(m)}}{\Delta y} \right] \\ + \frac{R_x \Delta_z H_y^{(m)}}{R_z \Delta z} - \frac{R_y \Delta_z H_x^{(m)}}{R_z \Delta z}. \quad (18) $$

The second line of the formula accounts for the inclination of the line formation plane with respect to the line of sight. Then, $J_x$ and $J_y$ can be derived by applying Eq. (7) to $J_x$, $J_y$, and $J_z$.

In the following, we have plotted the heliographic reference frame components, $J_{X,YZ}$, of the current density vector. For our plots, we derived the current density vector coordinates in the heliographic reference frame from the magnetic field unique vector (after disambiguation) rotated into the heliographic reference frame by applying Eq. (7).

Crouch et al. (2009) and Crouch (2013, 2015) developed a different calculation for also accounting for the departure from disk center, which implies that the line of sight is not perpendicular to the line formation plane. They also obtained that additional terms have to be introduced into the usual divergence expression. However, their formula given in Eqs. (4)–(5) of Crouch et al. (2009) and Crouch (2013, 2015) is different from our Eq. (9) because the corresponding reference frames are not the same. The divergence applied in Crouch et al. (2009) and Crouch (2013, 2015) uses spatial derivatives with respect to the heliographic X and Y and line of sight z coordinates, which do not form a rectangular reference frame; whereas our divergence uses spatial derivatives with respect to all the line of sight (x, y, z) reference frame coordinates. Our approach takes into consideration the fact that the observed lines are formed along inclined planes. Both approaches are different ways to treat the problem of the inclined line of sight.

### 4.4. The two-line formation depth difference

The quantity to be determined is the difference in the line of sight formation depth between the two lines. As it can be seen in Fig. 4 of Khomenko & Collados (2007), the two lines Fe I 6302.5 Å and 6301.5 Å are particularly interesting because their formation depths exhibit a strikingly parallel behavior. As also visible in Fig. 4 of Khomenko & Collados (2007), such a parallelism is not the case for the pair of Fe I 6302.5 and 5250.2 lines. This was the reason we discarded 5250.2, even though it is more sensitive to the magnetic field than 6301.5. Such a parallelism probably originates from the fact that the two lines 6302.5 and 6301.5 belong to the same multiplet n. 816 and have different gf values; these are, respectively, gf = 0.180 for 6301.5 and gf = 0.0627 for 6302.5 from the Kurucz data-basis. Since differential non-LTE effects within multiplets are thought to be very small, as proven by detailed, multi-level, non-LTE computations, this implies that the absorption coefficient of 6301.5 is three times larger than the absorption coefficient of 6302.5. There is then no doubt that 6301.5 forms higher than 6302.5. In other words, the optical depth of 6301.5 is almost 3, where the optical depth of 6302.5 is unity.

The two-iron-line formation difference in the quiet sun atmosphere was recently determined by a phase-shift analysis of HINODE observations by Faurobert et al. (2009). These authors obtained the observed value 63.2±0.9 km, which is corroborated by the value of 69 km derived by the same phase-shift technique applied to theoretical profiles computed with the non-LTE Uiten-broek’s code (Grec et al. 2010).

However, it has to be remarked that the formation height varies along the line profile. The highest value, which is typically several hundreds of km, is reached at line center; whereas in the far wings, the formation height is nearly 0 km in the visible range. However, in the inversion algorithm, the mean squared difference between observed and theoretical profiles is computed over the whole profile. As the Zeeman effect is maximum around
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Contribution functions of the two lines Fe I 6302.5 Å and 6301.5 Å as a function of the height above $\tau_{5000} = 1$. The heights where $\tau = 1$, which is the line formation height in the Eddington-Barbier approximation, are indicated by arrows.</figcaption></figure>

line center, this part of the profile is the most contributing to the field vector determination, so it is quite natural to assign the line center formation height to this field value. For instance, it has been found that the value of $1/2 \arctan U/Q$ computed at the line center provides the field azimuth within a very good first approximation.

In addition, it has to be remarked that the formation depth has a certain thickness, which is aptly represented by the behavior of the contribution function. We have plotted the contribution functions in Fig. 2 at line center for the two iron lines. This figure has been obtained by applying a non-LTE polarized radiation transfer code (Landi Degl'Innocenti et al. 1990; Bommier et al. 1991) to the Fe I 6302.5 and 6301.5 line formation in the quiet sun reference model atmosphere of Maltby et al. (1986) and in the absence of a magnetic field. It can be seen in the figure that the formation depth of each line is characterized by a level of accuracy that is on the order of ±75 km, which is the halfwidth of the contribution function. However, the difference $\Delta z = 66$ km between the two lines remains clearly visible all along the profile of the contribution function. In Fig. 2, we have also indicated with arrows the height of $\tau = 1$, which we consider as the line formation depth following the Eddington-Barbier approximation. It is well visible that this height correctly represents the mean line formation depth. In Fig. 4 of Khomenko & Collados (2007), we can also see that the formation depth difference remains on the order of $\Delta z = 70$ km in the quiet sun, whatever the actual formation depth is.

Finally, in Fig. 4 of Khomenko & Collados (2007), we can observe that the depth difference is larger in active regions than in quiet ones. We thus applied instead the depth difference of 98 km at the disk center for active regions. This is also the value determined with the non-LTE Uitenbroek's code (Grec et al. 2010) in a previous approximation.

Finally, we derived the line formation depths from a model that we describe below. The temperature, electron pressure, and gas pressure were first taken from an atmospheric model. We used the Maltby et al. (1986) quiet sun photospheric reference model extrapolated downwards beyond –70 km to –450 km

below the $\tau_{5000} = 1$ level (courtesy of Instituto de Astrofísica de Canarias (IAC)). The continuum absorption coefficient was evaluated according to the MagnetoActive Fraunhofer LInes Stokes Parameters Profiles (MALIP) code of Landi Degl'Innocenti (1976). This was done, in particular, by including $H^-$ bound-free, $H^-$ free-free, neutral hydrogen atom opacity, Rayleigh scattering on H atoms, and Thompson scattering on free electrons. The line absorption coefficient was derived from the Boltzmann and Saha equilibrium laws, taking the two first ions of element into account. The atomic data were taken from Wiese or Moore and the partition functions from Wittmann. The iron abundance was assumed to be 7.60 and the sodium abundance 6.25 in the usual logarithmic scale, where the abundance of hydrogen is 12. A depth-independent microturbulent velocity field of 1 km/s was introduced. Finally, departures from LTE in the ionization equilibrium were simulated for layers above $\tau_{5000} = 0.1$ by applying Saha's law with a constant radiation temperature of 5100 K, instead of the electron temperature provided by the atmospheric model.

At final step, the line center optical depth grid was scaled to the continuum optical depth grid by applying the respective absorption coefficients. We used the continuum optical depth grid provided in the atmospheric model and the transfer equation was not explicitly solved again. The height of formation of the line center was then determined as follows. Once obtained the grid of line center optical depths, the height of formation of the line center was located where the optical depth along the line of sight is unity (Eddington-Barbier approximation); namely, where $\tau/\mu = 1$ and $\tau$ is the line center's optical depth along the vertical and $\mu$ is the cosine of the heliocentric angle $\theta$ (here supposed to be 0). As seen in Fig. 5 of Bruls et al. (1991) and in Fig. 2 of the present paper, this conveniently represents the depth where the contribution function has its maximum.

We thus obtained quiet sun line center formation heights of 262 km for Fe I 6302.5 and 328 km for Fe I 6301.5 above the $\tau_{5000} = 1$ level. This leads to a difference of 66 km in excellent agreement with the measurements described above, which validates our computation method.

## 4.5. The aspect ratio of the strongly stratified atmosphere

We were indeed confronted with the problem that the vertical gradient of the magnetic field, $\mu_0 \partial H_z / \partial z$, was found to be on the order of 3 G/km in the spot umbrae, whereas the horizontal gradient $\mu_0 \partial H_x / \partial x + \mu_0 \partial H_y / \partial y$ was only on the order of 0.3 G/km. This leads to a non-vanishing value of $\text{div } \mathbf{H}$, whatever the signs would be. We then investigated the literature and found 15 references, which fully confirm these values provided by different instruments (ground-based as well as spaceborne), different inversion methods (SIR or others), and different spectral lines. The detailed description of these references can be found in Bommier (2013). An observation review was also presented by Balthasar (2018). All concluded to the above cited values. In other words, a loss of magnetic flux is observed with increasing height, which is not compensated for by an increase of the horizontal flux. In Bommier (2013, 2014), we showed that the lack of spatial resolution in both transverse and along the line of sight directions cannot be held responsible for the seemingly non-vanishing observed $\text{div } \mathbf{H}$. This was based on mathematical study of the convolution procedure. It was shown that the divergence computed with averaged quantities is equal to the average of the local divergences. Accordingly, if the local divergence is zero, the divergence computed with averaged quantities should also be zero, within the noise level. The question arised to know
if the observed value of div **H** is an effect of different spatial resolution along the different space directions. The effect of the spatial resolution (horizontal as well as vertical) is filtering. The mathematical demonstration that the filtered divergence is the divergence of the filtered quantity is given in Bommier (2013, 2014). Bommier (2014) included also a more straightforward demonstration in the spatial Fourier space. The same demonstration applies to the spatial averaging and the effect of the magnetic filling factor $\alpha$. As explained in the previous section, we only know the “local average magnetic field”, $\alpha H$, from the measurements, but the average divergence is the divergence of the averaged field. In addition, Bommier (2014) provides results of numerical tests devoted to investigate an eventual effect of the limited spatial resolution (eventually anisotropic). The numerical tests are all negative. They conclude with a zero divergence computed by finite differences, when the local divergence is zero. In our opinion, the negation of this logical proposition is that the non-zero observed value for the divergence indicates a non-zero local value. In Bommier (2015), we present a discussion about the necessity for div **H** to be zero. We argue that the existence of magnetic monopoles is not the only possibility for a non-zero div **H**. Following a discussion about the measurement noise level effect, we present below another possibility we have investigated.

In the THÉMIS measurements described in Bommier et al. (2007), which is the validation paper of the UNNOFIT inversion method, the polarimetric noise is assumed to be $1.5 \times 10^{-3}$. As reported in Sect. 2, this results in an inaccuracy of 10 G for the longitudinal field and 15 G for the transverse field. The field difference observed between and from the two line 6301.5 and 6302.5 is about 300 G, when the difference in line formation height is about 100 km, close to the disk center. This results in an inaccuracy of 0.2 G/km for the vertical field gradient observed at the disk center. As for the transverse field, the typical THÉMIS pixel size is about 500 km (Bommier 2013), which results in an inaccuracy of less than 0.2 G/km for the horizontal field gradient observed at the disk center. The total inaccuracy on the observed div **H** value results in 0.4 G/km, to be compared to the non-vanishing value of 2.7 G/km, which results from the measurements reported in the literature. The non-vanishing value observed for div **H** is then markedly higher than the noise level.

We then investigated how the magnetic field is influenced by the plasma anisotropy due to the strong stratification due to the gravity and the density at the star surface, which is responsible for an “aspect ratio” between horizontal and vertical typical lengths, respectively denoted as $l_h$ and $l_v$. In the case of the solar photosphere, by applying strongly stratified fluid mechanics laws following Brethouwer et al. (2007), we evaluated this aspect ratio to be on the order of $l_h/l_v = 20$ (Bommier 2013, 2014) in the quiet sun photosphere. In Bommier (2020), we showed that what is measured by Zeeman effect is the magnetic field **H** and not the magnetic induction **B**. The magnetic induction **B** is divergence-free, which implies that div **H** = −div **M**, where **M** is the magnetization. The magnetization is linked to the matter. Accordingly, it scales following the typical lengths introduced above. It reveals that once the inverse of the aspect ratio is applied to scale the different magnetic field components, the scaled div **H** vanishes, which enables the ambiguity resolution that we report in this paper. An example of this result in the case of NOAA 10808 observed with THÉMIS on September 13, 2005 is shown in Fig. 3 of Bommier (2013). Although we evaluated the theoretical value of the aspect ratio on the order of 20 in the quiet sun photosphere, this ratio may eventually be differ-

ent in sunspots or plages or active regions, as in the case of the regions we treated. This is not so widely known and sometimes we continued to adjust the ratio (kept constant in the whole map) until the ambiguity resolution was consistent across the whole map. Here, the word “consistent” is meant to indicate that the resulting map does not show any artificial strong current line, which would delineate sharp azimuth change due to ambiguous solution selection change. The resulting map displays (at most) homogeneous field directions. We obtain that this inverse scaling is a necessity for a correct ambiguity resolution. The correctness is obtained by comparison with the intuitive ambiguity solution, which can be found in sunspot penumbrae as discussed in the following subsection. We applied the Maxwell law curl $H = J$ for the current density derivation from the magnetic field vector in different regions reported in the following subsections.

## 4.6. Verification of the ambiguity resolution results

In this way, we treated 60 maps observed in Fe I 6301/6302 and 26 maps observed in Na I D with THÉMIS (in 2010–2013) and 23 HINODE/SOT/SP observations of active regions<fnref n="2" />. For the sunspots the correctness of the solution can be verified because the spot umbra polarity is known from circular polarization. This is true even near the limb, where the solar vertical field at spot center is nearly transverse with respect to the l.o.s. The ambiguity solution is good if it is in agreement with the polarity. It has to be remarked that for a sunspot observed anywhere on the disk there is always a place in it where the ambiguity solution can be derived from the spot polarity. This place is either in the penumbra, when the spot is observed near the disk center, or in the umbra, when the spot is observed near the limb. We consider such cases to be possibilities of observational proofs of our disambiguation method and we obtained successful proofs even in a spot at $\mu = 0.42$, which is $\theta = 65^\circ$ from the disk center.

## 4.7. Test of the method on theoretical data

Our method for taking the departure from disk center into account, which leads to additional terms in the divergence formula Eq. (9), has been tested on theoretical data. Divergence-free data were prepared by simulating two close sunspots forming a $\delta$-spot and by applying the magnetic field model described by Eqs. (8)–(10) of Fan & Gibson (2004). Simulation results based on this model have been used to study $\delta$-spots (as in Fig. 3 of Leka et al. 2005) and also to test various ambiguity resolution methods, as in Metcalf et al. (2006). The ambiguity of the transverse magnetic field in the line of sight reference frame was simulated and submitted to the ME0 ambiguity resolution code of Leka et al. (2009), which was modified by us for taking the vertical magnetic field gradient $\partial H_z/\partial z$ from two-line observations instead from extrapolation as described at the beginning of Sect. 4. The two-line observation had been simulated with theoretical results at two different depths separated by 70 km along the line of sight. The pixel sizes were assumed to be $\Delta x = 1160$ km and $\Delta y = 581$ km, as in our THÉMIS observations described in Bommier (2013). The $\delta$-spot was assumed to be far from disk center and located at $30^\circ$ from the central meridian in longitude and at $45^\circ$ in latitude. The solar disk center latitude was assumed to be zero.

The test result is represented in Fig. 3, where the theoretical reference field is plotted in the top map. The result of the disambiguation when the usual divergence formula is applied, which is given by the first line of Eq. (9) only, is displayed in the

<footnote n="2">² http://lesia.obspm.fr/perso/veronique-bommier/</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Result of the ambiguity resolution test using our method applied to theoretical data. Top: Map of the input magnetic field vector. Middle: ambiguity resolution by applying the usual divergence formula (first line of Eq. (9) only) to these data, where the region is assumed to be out of disk center and where two lines formed at two different altitudes are supposed to be observed. Bottom: Ambiguity resolution by applying Eq. (9).</figcaption></figure>

<table><caption>**Table 1.** Performance metrics for the limited resolution case for our DIVB resolution algorithm applied to the test data of Leka et al. (2009) and Crouch (2013).</caption><thead><tr><th>Pixel size</th><th>Line</th><th>M<sub>area</sub></th><th>M<sub>H&perp;&gt;100 G</sub></th><th>M<sub>H&perp;&gt;500 G</sub></th></tr></thead><tbody><tr><td>0.3"</td><td>Lower</td><td>0.99</td><td>1.00</td><td>1.00</td></tr><tr><td>0.3"</td><td>Upper</td><td>0.99</td><td>1.00</td><td>1.00</td></tr><tr><td>0.9"</td><td>Lower</td><td>0.99</td><td>1.00</td><td>1.00</td></tr><tr><td>0.9"</td><td>Upper</td><td>0.99</td><td>1.00</td><td>1.00</td></tr></tbody></table>

**Notes.** The data were provided and the metrics were computed by K.D. Leka's courtesy.

middle map. The minimized quantity is in fact $|\text{div } \mathbf{H}| + \lambda|J_z|$ with $\lambda = 1$ as in our method. It can be seen that the ambiguity resolution widely fails over the $\delta$-spot. When the full Eq. (9) is applied instead, the obtained result is displayed in the bottom map of Fig. 3, which is in total agreement with the reference field. Thus, we confirm via numerical test that the second line of Eq. (9) is necessary to resolve the ambiguity by applying a two-line analysis to maps observed out of disk center as discussed after Eq. (10) and in Fig. 1.

## 4.8. Test of the resolution code on theoretical data

We tested our ambiguity resolution code DIVB2 (for: divergence of B) built on the method described above on the theoretical test data described in Leka et al. (2009). These data were used for testing various ambiguity resolution methods as in Leka et al. (2009). The methods by Crouch et al. (2009) and Crouch (2013, 2015) also used magnetic field data at two different heights, which were applied to the current test data. We used two series of these data.

<table><caption>**Table 2.** Performance metrics for the noise-added case for our DIVB resolution algorithm applied to the test data provided by K.D. Leka's courtesy.</caption><thead><tr><th>Noise level</th><th>M<sub>area</sub></th><th>M<sub>H&perp;&gt;100 G</sub></th><th>M<sub>H&perp;&gt;500 G</sub></th></tr></thead><tbody><tr><td>No</td><td>1.00</td><td>1.00</td><td>1.00</td></tr><tr><td>Low</td><td>0.85</td><td>0.95</td><td>1.00</td></tr><tr><td>High</td><td>0.73</td><td>0.83</td><td>0.90</td></tr></tbody></table>

**Notes.** These data were similar to that described in Crouch (2013), but with a two line formation height difference of 890 km along the l.o.s. instead. The metrics were computed by K.D. Leka's courtesy.

The first series was the “flower” data, where the sunspots are simulated with form of flowers, as visible in Fig. 4 of Leka et al. (2009). These data were prepared to test the code robustness against lack of spatial resolution. The preparation of the theoretical data is described on p. 93 of Leka et al. (2009) and also on p. 111 of Crouch (2013). We used the theoretical data with averaged pixel sizes of 0.3" and 0.9". These theoretical data are provided at two heights separated by 0.18". These data simulate location at the disk center, namely, the line of sight is assumed to be perpendicular to the line formation planes. For these data, we obtained the result metrics listed in Table 1. The metrics is given by the fraction of pixels with correct ambiguity resolution in the whole map, $M_{\text{area}}$, and the fraction of pixels with correct ambiguity resolution within the conditions of transverse field $H_{\perp}$ stronger than 100 G and 500 G, respectively. As seen in Table 1, our code obtains excellent results for the two resolution cases.

The second series of data was the simulation represented in Fig. 1 of Leka et al. (2009). These data were artificially noised by adding them a theoretical photon noise, at the level of polarimetric accuracy of $10^{-3}$ for the low noise case and $10^{-2}$ for the high noise case, as described on p. 89 of Leka et al. (2009). The assumed pixel size was 0.5". The case without any added noise was also treated. These data are available at two heights. We used data at two heights distant of 890 km along the l.o.s., whereas Crouch (2013) used data at two heights distant of the pixel size as described in p. 111 of Crouch (2013), which is 376 km only. For these data, we obtained the result metrics listed in Table 2. In this case also, our code obtains very good results. The results are perfect in the case without any noise. The highest noise level, which is $10^{-2}$ in polarization, is rather high and higher than the current observation polarization inaccuracies. In the case of this high level noise, we obtain nevertheless rather good results better than 80% correct for the high noise level and better than 95% correct for the low noise level. As visible in Table 2, our worst results locate in weak magnetic field regions.

This second series of data were assumed to be located out of disk center at latitude 9° South and longitude 36° East with respect to the central meridian. This leads to an heliocentric angle cosine $\mu = \cos\theta = 0.80$, which is rather far from disk center. It has to be remarked that we obtain perfect ambiguity resolution result with our code in the no noise case as visible in the first line of Table 2. In our code, the disambiguation is performed in the l.o.s. reference frame and we applied Eq. (9) to account for departure from disk center in the case of lines formed along horizontal planes, which is done in the second line of Eq. (9). The success of our resolution of theoretical data again validates our Eq. (9) and, in particular, its second line in the case of lines formed along horizontal planes and observed out of disk center.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Active region NOAA 11420 observed by HINODE/SOT/SP on February 18, 2012 between 11:08 and 12:04 UT. The ambiguity is resolved by our method. Left: Magnetic field vector, vertical component between -2040 and +2040, horizontal component maximum arrow length 1355 G. Middle: current density vector, vertical component between -64 and +64, horizontal component typical arrow length 5000 mA/m². Right: Lorentz force vector, vertical component between -448 and +448, horizontal component typical arrow length 800 mN/m². The spatial resolution was reduced by a factor 5 × 5 for paper file size purposes (color figure on-line).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Active region NOAA 11494 observed by HINODE/SOT/SP on June 7, 2012 between 14:15 and 15:05 UT. The ambiguity is resolved by our method. Left: magnetic field vector, vertical component between -2150 and +2150, horizontal component maximum arrow length 1743 G. Middle: Current density vector, vertical component between -219 and +219, horizontal component typical arrow length 5000 mA/m². Right: Lorentz force vector, vertical component between -682 and +682, horizontal component typical arrow length 800 mN/m². The spatial resolution was reduced by a factor 5 × 5 for paper file size purposes (color figure on-line).</figcaption></figure>

## 5. Results: Examples of vector maps

The results presented below were obtained from 86 maps observed with THÉ MIS and 23 maps observed by HINODE/SOT/SP, which are visible in V. Bommier's personal page<fnref n="3" />. The magnetic field final values are also available on-line there. Four examples of HINODE data and two examples of THÉ MIS data are included in Figs. 4-9 below. For accuracy reasons, the HINODE data pixels were two times binned two by two before the ambiguity resolution, which results in pixel sizes of $Δx = 862$ km and $Δy = 928$ km. The THÉ MIS data pixel size was $Δx = 794$ km and $Δy = 725$ km, the binning was not performed.

The plotted Lorentz force displays the correct direction of the force, but may not actually be the exact force strength. What is, in fact, plotted in the following figures is $J × μ_0H$, when the Lorentz force is in fact $J × B$. The magnetic induction $B$, the magnetic field $H$ and the magnetization $M$ are all parallel, because

$$M = -\frac{(n_e + n_i)k_B T}{B^2} \mathbf{B}, \qquad (19)$$

where $n_e$ and $n_i$ are the electron and ion densities respectively, $k_B$ is the Boltzmann constant, and $T$ is the temperature (Delcroix & Bers 1994). The minus sign implies that $M$ is

antiparallel to $B$, which is plasma diamagnetism. This can be also written as

$$\mu_0 M = -\frac{\beta_c}{2} \mathbf{B}, \qquad (20)$$

where $β_c$ is the plasma $β$, ratio of kinetic energy to magnetic energy, however, this is expressed in terms of the charged particle density

$$\beta_c = \frac{(n_e + n_i)k_B T}{B^2/2\mu_0}. \qquad (21)$$

In the usual photospheric electron and ion densities, $μ_0M$ remains negligible with respect to $μ_0H$ and $B$, which are then equal, and the following figures would display the full Lorentz force. However, if the electron density is higher, following the suggestion of interior electron thermal escape by Bommier (2020), $μ_0M$ can become non-negligible, so that the length of the Lorentz force vectors would not be exact in the following figures, even if the vector direction is correct.

### 5.1. HINODE data: Examples of a single regular spot

We systematically observe a strong circular horizontal current component, which wraps spots clockwise about a positive polarity spot and counterclockwise about a negative polarity spot. As a consequence, the Lorentz force vector is systematically found

<footnote n="3">³ http://lesia.obspm.fr/perso/veronique-bommier</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Active region NOAA 11476 observed by HINODE/SOT/SP on May 12, 2012 between 02:38 and 03:11 UT. The ambiguity is resolved by our method. Left: Magnetic field vector, vertical component between -2450 and +2450, horizontal component maximum arrow length 2470 G. Middle: current density vector, vertical component between -218 and +218, and horizontal component typical arrow length 8000 mA/m². Right: Lorentz force vector, vertical component between -1040 and +1040, and horizontal component typical arrow length 800 mN/m². The spatial resolution was reduced by a factor 5 × 5 for paper file size purposes (color figure online).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Active region NOAA 11429 observed by HINODE/SOT/SP on March 6, 2012 between 22:10 and 22:43 UT. The ambiguity is resolved by our method. Left: Magnetic field vector, vertical component between -2480 and +2480, horizontal component maximum arrow length 2441 G. Middle: Current density vector, vertical component between -185 and +185, horizontal component typical arrow length 5000 mA/m². Right: Lorentz force vector, vertical component between -1180 and +1180, and horizontal component typical arrow length 800 mN/m². The spatial resolution was reduced by a factor 5 × 5 for paper file size purposes (color figure online).</figcaption></figure>

centripetal with respect to the spot center. Numerous examples can be found along the web page. We can cite in particular the active regions: NOAA 11420 observed by HINODE/SOT/SP on February 18, 2012 between 11:08 and 12:04 UT (negative polarity spot), which is represented in Fig. 4, and NOAA 11494 observed by HINODE/SOT/SP on June 7, 2012 between 14:15 and 15:05 UT (positive polarity spot), which is represented in Fig. 5. These examples can be found in the web page with a better spatial resolution. Similar examples can be found in the web page for HINODE/SOT/SP observations of single spots on December 12, 2006, September 7, 2011, February 1, 2012, May 1, 2012, and June 13, 2012.

## 5.2. HINODE data: Examples of a neutral line

We often observe a strong horizontal current component, which crosses the neutral line. Two examples are presented in this paper, which can also be found in the web page: NOAA 11476 observed by HINODE/SOT/SP on May 12, 2012 between 02:38

and 03:11 UT, which is represented in Fig. 6, and NOAA 11429 observed by HINODE/SOT/SP on March 6, 2012 between 22:10 and 22:43, which is represented in Fig. 7. We can refer to the web page for a better spatial resolution. NOAA 11476, which is represented in Fig. 6, produced numerous C-class flares and NOAA 11429, which is represented in Fig. 7, produced an X-class flare followed by a CME. The example of September 7, 2011, which is visible in the web page, is also the case of an active region that produced an X-class flare. A strong horizontal current component is well visible across the neutral line.

## 5.3. THÉMIS data: Example of a single regular spot

The example presented in this paper is the case of NOAA 11857 observed with THÉMIS on October 7, 2013, between 08:57 and 09:22 UT in Fe I 6301/6302. In the case of the photosphere, which is observed in Fe I 6301/6302 and represented in Fig. 8, we find again the current and Lorentz force characteristics described above; namely, circular current about the spot and
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Active region NOAA 11857 observed with THEMIS on October 7, 2013 between 08:57 and 09:22 UT, in the photospheric lines Fe I 6301/6302. The ambiguity is resolved by our method. Left: Magnetic field vector, vertical component between −1790 and +1790, and horizontal component maximum arrow length 1070 G. Middle: Current density vector, vertical component between −68 and +68, and horizontal component maximum arrow length 746 mA/m². Right: Lorentz force vector, vertical component between −56 and +56, and horizontal component maximum arrow length 133 mN/m² (color figure online).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Active region NOAA 11865 observed with THEMIS on October 11, 2013 between 09:39 and 09:59 UT, in the photospheric lines Fe I 6301/6302. The ambiguity is resolved by our method. Left: Magnetic field vector, vertical component between −1450 and +1450, and horizontal component maximum arrow length 1320 G. Middle: Current density vector, vertical component between −96 and +96, and horizontal component maximum arrow length 655 mA/m². Right: Lorentz force vector, vertical component between −83 and +83, and horizontal component maximum arrow length 43 mN/m² (color figure on-line).</figcaption></figure>

centripetal Lorentz force. NOAA 11857 was observed with THÉMIS on October 5, October 7 (two times, morning and afternoon), October 8, and October 9. All these examples are visible in the web page. The magnetic field is more rotating about the spot center at the chromospheric level than at the photospheric level, where it is more radial.

## 5.4. THÉMIS data: Example of a neutral line

The example presented in this paper is the case of NOAA 11865 observed with THÉMIS on October 11, 2013, between 09:39 and 09:59 UT in Fe I 6301/6302. Again, in the photosphere, which is represented in Fig. 9, we observe a strong horizontal component across the neutral line.

## 6. Conclusion

We have shown that observations in two different lines, which belong to the same multiplet but have different absorption coefficients and thus formed at two different depths (e.g., Fe I 6302.5 Å and 6301.5 Å) enable the resolution of the azimuth ambiguity that is left over from the Zeeman signal interpretation. The anisotropy of the strongly stratified plasma of the

photosphere has also to be accounted for following Bommier (2013, 2014). As a consequence, two ambiguity-resolved field vector maps are obtained at two different but close altitudes, which enable the derivation of the current density full vector via curl $H = J$. This reveals the horizontal component of the current density, which is found markedly stronger than the better known vertical one, as already observed by Puschmann et al. (2010a) from HINODE/SOT/SP data. We observe some systematic trends, such as circular currents wrapping spots clockwise about a positive polarity spot and counterclockwise about a negative polarity spot and strong horizontal current components, which cross active region neutral lines. The wrapping direction with respect to the spot polarity is the same as the average one in Fig. 1 of Puschmann et al. (2010a), where completely different methods were used. As a result, the Lorentz force may be computed. It is found to be centripetal in sunspots. The sunspots are thus maintained by the force.

We found that the anisotropic scaling of the real solar data following the aspect ratio in the strongly stratified medium of the solar photosphere, suggested by Bommier (2013, 2014), has to be applied for a correct resolution of the ambiguity. The correctness can be established by comparison with the intuitive ambiguity solution, which can be derived in spot penumbrae
from the spot polarity. The necessity of applying the aspect ratio for solving the ambiguity is another main result of the present work.

This can be done by applying an anisotropic scaling factor to the data before disambiguation. For instance, by artificially dividing by 10 the pixel sizes before submission to the Azimuthal Ambiguity Resolution (AMBIG2) code proposed by Crouch (2013) for two lines. This code also requires that the observed region be located close to disk center. Our disambiguation code is able to also disambiguate data of regions located out of the disk center.

Finally, we can remark that the Na I D1 and D2 lines form another favorable line pair, because they are also two lines of the same multiplet, but with different absorption coefficients. We evaluated their quiet sun line center formation heights at 533 km for Na I D1 and 604 km for Na I D2 above the $\tau_{5000} = 1$ level. These lines are then located in the low chromosphere close to the temperature minimum and with a difference of 71 km. Observing this line pair would open access to the low chromosphere where they are formed. However, a better spatial resolution for both observations and analysis would probably be necessary in such a medium.

In the present analysis, the spatial resolution of both observations and analysis is very limited. The Stokes inversion was performed with average over the whole observed pixel. The presence of unresolved magnetic fields was however roughly accounted for by means of a magnetic filling factor $\alpha$, which was introduced and determined in the inversion. However, active regions such as those studied here, which include sunspot umbrae and penumbrae, have a structure akin to dots that are smaller than the present pixel sizes. This concerns the horizontal spatial resolution, but the vertical spatial resolution of our study is limited also, because it is limited to the difference in line formation depth of two lines, which is found not so far from the vertical height scale. As a rough first approach, we assume a constant formation depth difference for the two lines under interest, although our images include various different structures, where this difference may vary. A better vertical spatial resolution would be attained by applying for instance the SIR inversion code (Ruiz Cobo & del Toro Iniesta 1992).

A series of other effects able to perturb the sunspot magnetic field measurements is detailed in the above Sect. 3, including: velocity or magnetic field-temperature correlations, blending by molecular lines, straylight from penumbra, departure from a purely hydrostatic equilibrium, non-horizontal and corrugated $\tau = 1$ surface, and non-LTE effects in spectral lines. Consequently, as long as these and other effects have not been properly tested in realistic model atmospheres and with state-of-the-art radiative transfer computations, the scaling factor we introduce must be taken as a purely empirical factor.

**Acknowledgements.** The author is grateful to an anonymous referee for refined suggestions, reported in this paper. Some reported observations were made with the French-Italian telescope THÉMIS operated by the CNRS on the island of Tenerife in the Spanish Observatorio del Teide of the Instituto de Astrofísica de Canarias. We are grateful to B. Lites for the level-1 HIN-ODE/SOT/SP data. For the data inversion, this work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the *Investissements d'Avenir* program supervised by the *Agence Nationale pour la Recherche*.

## References

Balthasar, H. 2018, Sol. Phys., 293, 120

Bommier, V. 2013, Phys. Res. Int., 2013, 195403

Bommier, V. 2014, C.R. Phys., 15, 430

Bommier, V. 2015, in IAU Symposium, eds. K. N. Nagendra, S. Bagnulo, R. Centeno, & M. Jesús Martínez González, Polarimetry, 305, 28

Bommier, V. 2020, A&A, 634, A40

Bommier, V., Landi Degl'Innocenti, E., & Sahal-Bréchot, S. 1991, A&A, 244, 383

Bommier, V., Landi Degl'Innocenti, E., Landolfi, M., & Molodij, G. 2007, A&A, 464, 323

Brethouwer, G., Billant, P., Lindborg, E., & Chomaz, J. M. 2007, J. Fluid Mech., 585, 343

Bruls, J. H. M. J., Lites, B. W., & Murphy, G. A. 1991, in Solar Polarimetry, Proceedings of the 11st National Solar Observatory/Sacramento Peak Summer Workshop, Sunspot, New Mexico, 27-31 August 1990, ed. L. J. November (National Solar Observatory), 444

Buehler, D., Lagg, A., Solanki, S. K., & van Noort, M. 2015, A&A, 576, A27

Crouch, A. D. 2013, Sol. Phys., 282, 107

Crouch, A. D. 2015, Sol. Phys., 290, 2677

Crouch, A. D., Barnes, G., & Leka, K. D. 2009, Sol. Phys., 260, 271

Delcroix, J. L., & Bers, A. 1994, Physique des Plasmas (7, rue de l'Estrapade, 75005 Paris: InterÉditions/CNRS Éditions)

Fan, Y., & Gibson, S. E. 2004, ApJ, 609, 1123

Faurobert, M., Aime, C., Périni, C., et al. 2009, A&A, 507, L29

Frutiger, C., Solanki, S. K., Fligge, M., & Bruls, J. H. M. J. 2000, A&A, 358, 1109

Gary, G. A., & Hagyard, M. J. 1990, Sol. Phys., 126, 21

Grec, C., Uitenbroek, H., Faurobert, M., & Aime, C. 2010, A&A, 514, A91

Khomenko, E., & Collados, M. 2007, ApJ, 659, 1726

Landi Degl'Innocenti, E. 1976, A&AS, 25, 379

Landi Degl'Innocenti, E., Bommier, V., & Sahal-Bréchot, S. 1990, A&A, 235, 459

Landolfi, M., & Landi Degl'Innocenti, E. 1982, Sol. Phys., 78, 355

Landolfi, M., Landi Degl'Innocenti, E., & Arena, P. 1984, Sol. Phys., 93, 269

Leka, K. D., Fan, Y., & Barnes, G. 2005, ApJ, 626, 1091

Leka, K. D., Barnes, G., Crouch, A. D., et al. 2009, Sol. Phys., 260, 83

Leka, K. D., Wagner, E. L., Griñón-Marín, A. B., Bommier, V., & Higgins, R. E. L. 2022, Sol. Phys., 297, 121

Maltby, P., Avrett, E. H., Carlsson, M., et al. 1986, ApJ, 306, 284

Metcalf, T. R. 1994, Sol. Phys., 155, 235

Metcalf, T. R., Leka, K. D., Barnes, G., et al. 2006, Sol. Phys., 237, 267

Molodij, G., Bommier, V., & Rayrole, J. 2011, A&A, 531, A139

Puschmann, K. G., Ruiz Cobo, B., & Martínez Pillet, V. 2010a, ApJ, 721, L58

Puschmann, K. G., Ruiz Cobo, B., & Martínez Pillet, V. 2010b, ApJ, 720, 1417

Ruiz Cobo, B., & del Toro Iniesta, J. C. 1992, ApJ, 398, 375

Socas-Navarro, H. 2005, ApJ, 633, L57

Stenflo, J. O. 1973, Sol. Phys., 32, 41

Vernazza, J. E., Avrett, E. H., & Loeser, R. 1973, ApJ, 184, 605
# Appendix A: Expressions of the divergence and curl when the quantities are measured along an inclined line formation plane

We describe below how divergence and curl have to be computed, when the measurements are made in pixels, and when the spectral line used for the measurements is formed at a certain depth in the stellar atmosphere. This depth defines a formation plane for the line. In the following we neglect an eventual thickness of this plane, which we consider as infinitely sharp. We consider the general case where this plane is inclined with respect to the line of sight for an observation performed out of disk center. The case of disk center observation is however possible as a particular case.

We denote the heliographic reference frame as OXYZ and the l.o.s. reference frame as Oxyz with Oy solar north oriented. The point O where divergence and curl are computed is taken at longitude-latitude (L, b). The solar radius unit vector, **R**, in O is perpendicular to the line formation plane, which is tangent to the solar surface. The **R** coordinates in the l.o.s. reference frame are given in Eq. (1). In the l.o.s. reference frame, the line formation plane, which is perpendicular to **R** of coordinates ($R_x, R_y, R_z$), has then for general equation

$$R_x x + R_y y + R_z z = 0. \qquad (A.1)$$

We use $\mathbf{H}^{(1)}$ to denote the magnetic field vector measured as a function of x and y along the formation plane of the line number (1). $\mathbf{H}^{(1)}$ depends on the following variables

$$\mathbf{H}^{(1)}(x, y, z; R_x x + R_y y + R_z z = 0), \qquad (A.2)$$

where the semicolon means “such as”. A similar law holds for the magnetic field vector $\mathbf{H}^{(2)}$ measured with the second line number (2). In the l.o.s. reference frame Oxyz, we denote as $\Delta x$ and $\Delta y$ the distance between two neighboring pixels along the Ox and Oy axes respectively in the “sky plane”. Let’s denote as $\Delta_x H_x^{(1)}$ the variation of $H_x^{(1)}$ between two neighboring pixels in x direction (i + 1, j) and (i, j) distant of $\Delta x$ in the image

$$\Delta_x H_x^{(1)} = H_x^{(1)}(i+1, j) - H_x^{(1)}(i, j), \qquad (A.3)$$

and analogously

$$\Delta_y H_y^{(1)} = H_y^{(1)}(i, j+1) - H_y^{(1)}(i, j). \qquad (A.4)$$

For the second line one can similarly define $\Delta_x H_x^{(2)}$ and $\Delta_y H_y^{(2)}$. The magnetic field may be averaged between the two lines, with the acute angle method when the ambiguity is not resolved, and in this case we denote the average as $\Delta_x H_x^{(m)}$ and $\Delta_y H_y^{(m)}$. As for the z variation, it involves the two lines as

$$\Delta_z H_z^{(m)} = H_z^{(1)}(i, j) - H_z^{(2)}(i, j). \qquad (A.5)$$

## A.1. Calculation of the divergence

The mathematical expression of div $\mathbf{H}$ is frame-independent

$$\mathrm{div} \, \mathbf{H} = \frac{\partial H_x}{\partial x} + \frac{\partial H_y}{\partial y} + \frac{\partial H_z}{\partial z}. \qquad (A.6)$$

However, a problem arises when it is discretized in the above described system of line formation planes. $\Delta_x H_x^{(1)}/\Delta x$, which could a priori be considered as the $\partial H_x/\partial x$ contribution, includes also a variation of $H_x^{(1)}$ along the Oz axis because the magnetic

<figure><img src="image_11.png" /><figcaption>**Fig. A.1.** Reference frame for the divergence evaluation: case of an observation at the solar equator. The heliocentric $\theta$ angle is oriented by the Oy axis and is negative in the figure. The two formation planes corresponding to the two lines (1) and (2) are indicated by their respective number. The unit S vector is indeed the same as the unit R vector along the local solar radius.</figcaption></figure>

field is measured along the line formation plane. As a consequence, $\Delta_x H_x^{(1)}/\Delta x$ cannot finally be considered as an approximation of $\partial H_x/\partial x$ because $\partial H_x/\partial x$ is the variation of $H_x$ with x at constant y and z following the mathematical definition.

This appendix is devoted to accounting for this difficulty by applying the Ostrogradski’s theorem, which is that the volume integrated divergence equals the flux of the quantity through the surface that borders the volume. We applied this theorem to the small volume delineated by four neighboring pixels separated by $\Delta x$ and $\Delta y$ in the sky plane xOy, but located in the line formation plane, and the corresponding four other pixels for the second line formed in another formation plane separated from the first one by the distance $\Delta z$ along the line of sight. In the following, we detail three different derivations from the simplest to the most complex. We first consider an observation located at the solar equator without inclination of the solar rotation axis (Sect. A.1.1). We then generalize this demonstration to the general case (Sect. A.1.2) and we add a more elegant derivation leading to the same result as a third part (Sect. A.1.3).

### A.1.1. Case of a region located at solar equator

We first assume that O is located on solar equator and that the Sun rotation axis lies in the plane of the sky, i.e. that $b_0 = 0$. The heliocentric angle is denoted as $\theta$. A section of the small volume introduced above is represented in Fig. A.1. For each side of the volume, a unit vector perpendicular to the surface and pointing towards the exterior of the volume can be plotted. Three of these vectors are plotted in Fig. A.1, S, S', and S''. One has to calculate the flux of the magnetic field through each elementary surface. Consider first the surface element associated to the S vector. This surface element is part of the formation plane of line number (1), whereas S' and S'' join the two line formation planes. In the l.o.s. reference frame Oxyz, the S unit vector has the following coordinates:

$$S = \begin{cases} S_x = \sin \theta, \\ S_y = 0, \\ S_z = \cos \theta. \end{cases} \qquad (A.7)$$

As a consequence, the flux of the magnetic field through this surface per surface unit is

$$\mathbf{H}^{(1)} \cdot \mathbf{S} = H_z^{(1)} \cos \theta + H_x^{(1)} \sin \theta. \qquad (A.8)$$
<figure><img src="image_12.png" /><figcaption>Fig. A.2. Reference frame for the divergence evaluation: General case.</figcaption></figure>

As visible in Fig. A.1, the elementary surface value is $\Delta\ell\Delta y$, where

$$ \Delta\ell = \frac{\Delta x}{\cos \theta}. \qquad (A.9) $$

The flux through this elementary surface is then finally

$$ H^{(1)} \cdot S \Delta\ell \Delta y = (H_z^{(1)} + H_x^{(1)} \tan \theta) \Delta x \Delta y. \qquad (A.10) $$

The flux through the elementary surface associated to $S'$ is

$$ H^{(m)} \cdot S' \Delta z \Delta y = H_x^{(m)} \Delta z \Delta y, \qquad (A.11) $$

where the superscript (m) denotes the average between the two lines (1) and (2). Finally, the flux through the third elementary surface $S''$ is

$$ H^{(m)} \cdot S'' \Delta\ell \Delta z \cos \theta = H_y^{(m)} \Delta x \Delta z. \qquad (A.12) $$

As the elementary volume is $\Delta\ell\Delta y\Delta z \cos \theta = \Delta x\Delta y\Delta z$, one has finally in application of Ostrogradski's theorem

$$ \begin{aligned} \text{div } \mathbf{H} \Delta x \Delta y \Delta z &= \left[ (H_z^{(1)} - H_z^{(2)}) + (H_x^{(1)} - H_x^{(2)}) \tan \theta \right] \Delta x \Delta y \\ &\quad + \left[ H_x^{(m)}(x + \Delta x, y) - H_x^{(m)}(x, y) \right] \Delta z \Delta y \\ &\quad + \left[ H_y^{(m)}(x, y + \Delta y) - H_y^{(m)}(x, y) \right] \Delta x \Delta z, \end{aligned} \qquad (A.13) $$

which results in Eq. (10) of the paper. The term in $\tan \theta$ accounts for the line formation plane inclination with respect to the line of sight.

### A.1.2. Generalization

One considers now a point O on the solar surface with longitude L and latitude b. The solar radius in O denoted as R is perpendicular to the line formation plane. The R coordinates in terms of (L, b) are given in Eq. (1). There is a common perpendicular to this radius R and to the l.o.s. Oz. We denote as $Oy' = OY$ this common perpendicular, which is represented in the left part of Fig. A.2. In this section, the OXYZ reference frame is different from the rest of the paper, because the common perpendicular OY is not necessarily the solar meridian. $Oy'$ is also different from the Oy axis of the l.o.s. reference frame, therefore we use a different notation. In the $Ox'y'z$ reference frame the geometry is

the same as in the equatorial case studied above, so that its final result may be applied

$$ \begin{aligned} \text{div } \mathbf{H} = & \frac{\Delta_{x'} H_{x'}^{(m)}}{\Delta x'} + \frac{\Delta_{y'} H_{y'}^{(m)}}{\Delta y'} + \frac{\Delta_{z'} H_{z'}^{(m)}}{\Delta z'} \\ & + \tan \theta \frac{\Delta_{z'} H_{z'}^{(m)}}{\Delta z}, \end{aligned} \qquad (A.14) $$

where the $\Delta$ values are taken along each respective axis.

One has now to rotate this result into the l.o.s. reference frame represented in the right part of Fig. A.2. In this reference frame, the $i_{x'}$ unit vector along $Ox'$ has the following coordinates

$$ i_{x'} = \begin{cases} \cos \varphi \\ \sin \varphi \\ 0 \end{cases}. \qquad (A.15) $$

Accordingly,

$$ H_{x'}^{(m)} = H_x^{(m)} \cos \varphi + H_y^{(m)} \sin \varphi. \qquad (A.16) $$

One has by rotational invariance

$$ \frac{\Delta_{x'} H_{x'}^{(m)}}{\Delta x'} + \frac{\Delta_{y'} H_{y'}^{(m)}}{\Delta y'} = \frac{\Delta_x H_x^{(m)}}{\Delta x} + \frac{\Delta_y H_y^{(m)}}{\Delta y}, \qquad (A.17) $$

so that finally

$$ \begin{aligned} \text{div } \mathbf{H} = & \frac{\Delta_x H_x^{(m)}}{\Delta x} + \frac{\Delta_y H_y^{(m)}}{\Delta y} + \frac{\Delta_z H_z^{(m)}}{\Delta z} + \frac{1}{\cos \theta} \\ & \times \left[ \frac{\Delta_z H_x^{(m)}}{\Delta z} \sin \theta \cos \varphi + \frac{\Delta_z H_y^{(m)}}{\Delta z} \sin \theta \sin \varphi \right]. \end{aligned} \qquad (A.18) $$

The R vector coordinates can also be expressed as

$$ \mathbf{R} = \begin{cases} \sin \theta \cos \varphi \\ \sin \theta \sin \varphi \\ \cos \theta \end{cases}. \qquad (A.19) $$

Accordingly, the above Eq. (A.18) simply results in Eq. (9) of the paper. The second line accounts for the line formation plane inclination with respect to the line of sight

### A.1.3. General case

Let's introduce the l.o.s. basic unit vectors i, k and l (l.o.s. vector), and let's consider their projections $i_p$ and $k_p$ along Oz onto the line formation plane perpendicular to the solar radius unit vector R (see Fig. A.3). $i_p$ and $k_p$ are not unit vectors and have to be evaluated as follows. Let's denote by A the extremity of the i vector, and by A' the extremity of the $i_p$ vector. A' is the projection of A along Oz. The AA' line of Fig. A.3 is defined by the intersection of the two planes of equations $x = 1$ and $y = 0$, respectively. A' is the intersection of the AA' line with the line formation plane perpendicular to R whose equation is given above in Eq. (A.1). The coordinates of $i_p$ are then

$$ i_p = \begin{cases} 1 \\ 0 \\ -\frac{R_x}{R_z}. \end{cases} \qquad (A.20) $$

Analogously for the k and $k_p$ vectors, the BB' line of Fig. A.3 is defined by the intersection of the two planes of equations $x = 0$
<figure><img src="image_13.png" /><figcaption>**Fig. A.3.** Projection of the unit vectors of the line of sight reference frame onto the line formation plane, which is perpendicular to the solar radius unit vector, **R**.</figcaption></figure>

and $y = 1$, respectively, and the coordinates of $k_p$ are obtained by intersection with the line formation plane given in Eq. (A.1)

$$k_p \begin{cases} 0 \\ 1 \\ -\frac{R_y}{R_z} \end{cases}. \qquad (A.21)$$

We go on to apply the Ostrogradski’s theorem. The elementary volume is obtained by projecting the elementary surface $\Delta x \Delta y$ along $Oz$ onto the map plane, and by giving it a $\Delta z$ thickness along $Oz$. As visible in Fig. A.3, the surface of the projected $\Delta x \Delta y$ is $i_p \Delta x \times k_p \Delta y$, which has for value $\Delta x \Delta y / R_z = \Delta x \Delta y / \cos \theta$, where $\theta$ is the heliocentric angle. The surfaces of the other sides of the elementary volume are respectively $\Delta x \Delta z$ and $\Delta y \Delta z$ (they are obtained by taking the product of $i_p$ or $k_p$ with $\ell$, which is the unit vector along the l.o.s.). The divergence is then obtained as

$$\begin{aligned} \operatorname{div} \mathbf{H} \Delta x \Delta y \Delta z = & \left[ \mathbf{H}^{(m)}(z + \Delta z) - \mathbf{H}^{(m)}(z) \right] \\ & \cdot \left[ i_p \Delta x \times k_p \Delta y \right] \\ & + \left[ \mathbf{H}^{(m)}(x + \Delta x) - \mathbf{H}^{(m)}(x) \right] \\ & \cdot \left[ k_p \Delta y \times \ell \Delta z \right] \\ & + \left[ \mathbf{H}^{(m)}(y + \Delta y) - \mathbf{H}^{(m)}(y) \right] \\ & \cdot \left[ \ell \Delta z \times i_p \Delta x \right], \end{aligned} \qquad (A.22)$$

leading to the expression given in Eq. (9). For clarity, the repeated and unchanged indices have been omitted in $\mathbf{H}^{(m)}$.

<figure><img src="image_14.png" /><figcaption>**Fig. A.4.** Elementary surface for the calculation of $J_z$ by applying the Stokes theorem to the circulation of the magnetic field vector along the oriented frontier of the surface. This figure is another view of Fig. A.1. Accordingly, the heliocentric $\theta$ angle is oriented by the $Oy$ axis and is negative in the figure. The unit $S$ vector is indeed the same as the unit $R$ vector along the local solar radius.</figcaption></figure>

## A.2. Calculation of the curl

The curl is similarly obtained by applying Stokes’ theorem, which is that the flux of the curl through an elementary surface equals the “circulation” or work of the vector under study along the frontier of the surface, which is oriented as the surface itself. The elementary surfaces are the different faces of the elementary volume introduced above for the divergence. Analogously to the divergence derivation, we present first below a detailed calculation in the particular case of solar equator (Sect. A.2.1) and then the general case (Sect. A.2.2).

### A.2.1. Case of a region located at solar equator

We consider the calculation of $J_z$ from the circulation of $\mathbf{H}$ along the oriented frontier of the elementary surface perpendicular to $OZ$. We have represented this elementary surface in Fig. A.4, which is another view of Fig. A.1. The circulation of $\mathbf{H}$ about the oriented surface is given by the scalar products

$$\begin{aligned} & \left[ \mathbf{H}^{(m)}(y) - \mathbf{H}^{(m)}(y + \Delta y) \right] \cdot \mathbf{I} \Delta \ell \\ & + \left[ \mathbf{H}^{(m)}(x + \Delta x) - \mathbf{H}^{(m)}(x) \right] \cdot \mathbf{K} \Delta y. \end{aligned} \qquad (A.23)$$

For clarity, the repeated and unchanged indices have been omitted in $\mathbf{H}^{(m)}$. The unit vector $\mathbf{I}$ along $OX$ has the following coordinates in the l.o.s. reference frame $Oxyz$

$$\mathbf{I} \begin{cases} \cos \theta \\ 0 \\ -\sin \theta \end{cases}, \qquad (A.24)$$
whereas **K** is

$$ K \begin{cases} 0 \\ 1 \\ 0 \end{cases} . \qquad (A.25) $$

As a consequence, the circulation of **H** is

$$ \begin{aligned} & [H_x^{(m)}(y) \cos \theta - H_z^{(m)}(y) \sin \theta] \Delta \ell \\ & + [H_y^{(m)}(x + \Delta x)] \Delta y \qquad (A.26) \\ & - [H_x^{(m)}(y + \Delta y) \cos \theta - H_z^{(m)}(y + \Delta y) \sin \theta] \Delta \ell \\ & - [H_y^{(m)}(x)] \Delta y . \end{aligned} $$

This circulation equals the flux of **J** through the surface, which is $J_Z$ times the surface value, which is $\Delta\ell\Delta y = \Delta x\Delta y/\cos\theta$ following Eq. (A.9). This can be written as

$$ J_Z = \cos \theta \left[ \frac{\Delta_x H_y^{(m)}}{\Delta x} - \frac{\Delta_y H_x^{(m)}}{\Delta y} \right] \\ + \sin \theta \frac{\Delta_y H_x^{(m)}}{\Delta y}, \qquad (A.27) $$

which is Eq. (14) of the paper in the case of the surface located at equator, where the coordinates of the **R** vector are given by Eq. (A.19) with $\varphi = 0$.

### A.2.2. General case

Referring now to Fig. A.3, the sides of the elementary volume are respectively $i_p \Delta x$, $k_p \Delta y$ and $\ell \Delta z$. The circulation of **H** is given by its scalar product with each of these elementary vectors, eventually reversed following the surface orientation. This leads to

$$ \begin{aligned} J_Z \frac{\Delta x \Delta y}{R_z} &= \mathbf{J} \cdot [i_p \Delta x \times k_p \Delta y] \qquad (A.28) \\ &= [\mathbf{H}^{(m)}(y) - \mathbf{H}^{(m)}(y + \Delta y)] \cdot i_p \Delta x \\ &\quad + [\mathbf{H}^{(m)}(x + \Delta x) - \mathbf{H}^{(m)}(x)] \cdot k_p \Delta y \end{aligned} $$

which results in the expression given in Eq. (14). The flux through the two other faces, which are respectively perpendicular to the $Ox$ and $Oy$ of the line of sight reference frame, gives the current components along those axes

$$ \begin{aligned} J_y \Delta x \Delta z &= \mathbf{J} \cdot [\ell \, dz \times i_p \, \Delta x] \qquad (A.29) \\ &= [\mathbf{H}^{(m)}(z + \Delta z) - \mathbf{H}^{(m)}(z)] \cdot i_p \, \Delta x \\ &\quad + [\mathbf{H}^{(m)}(x) - \mathbf{H}^{(m)}(x + \Delta x)] \cdot \ell \, \Delta z \end{aligned} $$

and

$$ \begin{aligned} J_x \Delta y \Delta z &= \mathbf{J} \cdot [k_p \Delta y \times \ell \, \Delta z] \qquad (A.30) \\ &= [\mathbf{H}^{(m)}(z) - \mathbf{H}^{(m)}(z + \Delta z)] \cdot k_p \, \Delta y \\ &\quad + [\mathbf{H}^{(m)}(y + \Delta y) - \mathbf{H}^{(m)}(y)] \cdot \ell \, \Delta z . \end{aligned} $$

This results in the expressions given in Eqs. (15-16). Again, for clarity, the repeated and unchanged indices have been omitted in $\mathbf{H}^{(m)}$. These three derived current density vector coordinates are not all taken in the same reference frame. The frame rotation is fully accounted for in Sect. 4.3, where all the coordinates of the current density vector are finally derived in both reference frames.