11institutetext: LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne
Université, Université de Paris
5, Place Jules Janssen, 92190 Meudon, France

# Magnetometry: solving the fundamental ambiguity from line pair inversion

Véronique Bommier<EMAIL_ADDRESS>
(Received … / Accepted …)

###### Abstract

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

Aims. We show that observations in two different lines belonging to the same
multiplet but having different absorption coefficients so that they are formed
at two different depths like Fe I 6302.5 Å and 6301.5 Å, enable the resolution
of the azimuth ambiguity remaining from the Zeeman signal interpretation.

Methods. What is measured by interpretation of the Zeeman effect is the
magnetic field $\@vec{H}$, and not the divergence-free magnetic induction
$\@vec{B}$. We analyze how the anisotropy of the photosphere, which is
strongly stratified due to gravity and density at the star surface, affects
$\mathop{\mathrm{d}iv}\@vec{H}$ and how the ambiguity resolution has to be
performed in these conditions.

Results. 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 $\mathop{\mathrm{r}ot}\@vec{H}=\@vec{J}$.
This reveals the horizontal component of the current density, which is
generally found markedly larger than the better known vertical one. We observe
some systematical trends, of which we present examples in the paper, like
circular currents wrapping spots clockwise about a positive polarity spot and
anticlockwise about a negative polarity spot and strong horizontal current
components crossing active region neutral lines.

Conclusions. We finally remark that the Na I D1 and D2 lines form another such
line pair and we successfully tested their application to full vector mapping
(magnetic field and current density) with THÉMIS telescope observations.
Therefore, we propose them as an access to the low chromosphere where they are
formed. Examples of such observations are also reported in the paper.

###### Key Words.:

Magnetic fields – Polarization – Sun: magnetic fields – Sun: photosphere –
Sun: chromosphere – Sun: magnetic topology

††offprints: V. Bommier,

##  1 Introduction

In solar atmosphere, magnetic field and electric currents are strongly coupled
via the Ampère’s circuital law. Moreover, and without entering in details, the
presence of strong electric currents, for instance under the form of a current
sheet, could be the place where an eruptive event could occur. Thus, the
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
$\mathop{\mathrm{r}ot}\@vec{H}=\@vec{J}$. This makes them particularly
difficult to derive from the observations because this requires a good
accuracy for the primitive quantity (in the mathematical meaning), which is
the magnetic field. This requires also that the magnetic field vector is fully
known without any 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.
They are not discriminated at the spectropolarimetric data inversion step.
This ambiguity is often referred to as the azimuth ambiguity, or
$180{{}^{\circ}}$ ambiguity, or fundamental ambiguity.

The derivation of the electric current full vector requires that the ambiguity
is 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 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 in the spot. The same method was applied to HINODE/SOT/SP data by
Puschmann et al. (2010). The quality of HINODE data permitted to show a
coherence 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 that
harbor strong upflows and radial outflows (the intraspines).

Indeed, the derivation of the three coordinates of the current density vector
requires the knowledge of the depth variations of the magnetic field, which is
rarely achieved. In general, a single line is observed, which is formed at a
given depth. Except when the SIR code or a code of similar technique is
applied, the magnetic field is derived at this single depth and only the
vertical component $J_{z}$ of the current density can be finally derived.
Thus, full vector current density determinations remain rare. We present below
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
$\mathop{\mathrm{d}iv}\@vec{H}$ from the observations, which enables the
ambiguity resolution, the method of which is described in section 3. In
particular, we analyze how the anisotropy of the photosphere, which is
strongly stratified due to gravity and density at the star surface, affects
$\mathop{\mathrm{d}iv}\@vec{H}$ and how the ambiguity resolution has to be
performed in these conditions. We present example results in section 4. As
Puschmann et al. (2010), we obtain that the horizontal component of the
current density vector is much stronger than the vertical one. The vertical
component is the most frequently published and discussed. The horizontal
component remained generally unknown, which is yet much larger. We discuss
examples of typical structures we derived for the horizontal component, in
sunspot penumbrae and across some neutral lines.

We tested our multiline ambiguity solution method on different spectral lines
observed with the THÉMIS telescope operating in its new multislit mode. The
two splitted beams enter the spectrograph by two parallel slits. As a
consequence, the different spectral windows (Fe I 6301/6302, Na I D1, Na I D2)
are not simultaneously observed but sequentially. As a result, we obtained
that the pair of lines must belong to the same multiplet in order to be formed
in the same manner, but with different absorption coefficients. This is 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. We observed that the Na I D1 and D2 line pair
obeys the same conditions. But their depth formation is higher in the solar
atmosphere. It is located near the temperature minimum in the low
chromosphere. Thus, they could provide the disambiguated magnetic field at
this altitude. With the THÉMIS telescope, we observed also this line pair.
After the magnetic inversion of the spectropolarimetric observations, we
successfully resolved the ambiguity with this line pair. In sections 4.3 and
4.4, we present same regions observed with the Fe I 6301/6302 line pair
(photospheric) and with the Na I D1/D2 line pair (chromospheric). We present
the magnetic field maps observed in a sunspot penumbra and above a neutral
line at the two altitude ranges, photosphere and low chromosphere
respectively. In Sect. 5, we present a verification of our method, by applying
available codes to the same solar data and by comparing the results.

It has to be recalled that what is measured by interpretation of the Zeeman
effect is the magnetic field $\@vec{H}$, and not the magnetic induction
$\@vec{B}$ (Bommier 2020). The magnetic induction $\@vec{B}$ is divergence-
free, when the magnetic field $\@vec{H}$ is not. They are related by
$\@vec{B}=\mu_{0}(\@vec{H}+\@vec{M})$, where $\@vec{M}$ is the magnetization.
However, in usual solar models, the reported electron density is such that the
modulus of $\@vec{M}$ is very small, with respect to the modulus of
$\@vec{H}$, which results in $\@vec{B}\approx\mu_{0}\@vec{H}$, which makes
$\@vec{H}$ quite divergence-free also. In this respect, minimizing
$\mathop{\mathrm{d}iv}\@vec{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, and in addition 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, $\mathop{\mathrm{d}iv}\@vec{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 will apply Gauss units for describing the magnetic field. They
are in fact units for $\mu_{0}\@vec{H}$.

##  2 Data Analysis First Step: Milne-Eddington Inversion

We applied the Milne-Eddington inversion code UNNOFIT. Based on the Unno-
Rachkovsky solution of the transfer equation for the radiation Stokes
parameters of a spectral 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 fitting
the observed profile. Eight parameters were entered in the algorithm, which
are namely: the line center frequency $\lambda_{0}$, the Doppler width
$\Delta\lambda_{\mathrm{D}}$, $\eta_{0}$, the three magnetic field
coordinates, the Voigt parameter $a$, and the ratio $B_{1}/B_{0}$ of the two
parameters characterizing the Milne-Eddington atmosphere. This was the
original code developed by Landolfi et al. (1984).

However, since the pioneering work by Stenflo (1973), it is well known that
the atmosphere is permeated by unresolved magnetic structures, which can be
roughly modelled with a magnetic filling factor $\alpha$. Stenflo’s analysis
is based on magnetic field measurements from the two Fe I 5247.1 Å and 5250.2
Å lines, which are two lines of the same multiplet but with identical
absorption coefficients, so that they are formed at the same depth. They have
different Landé factors and then different magnetic sensitivities. This
enabled the unresolved magnetic structures discovery. We introduced the
magnetic filling factor $\alpha$ as the ninth parameter of the UNNOFIT
inversion, in the Levenberg-Marquardt algorithm. In other words, denoting by
‘$\mathrm{m}$’ and ‘$\mathrm{nm}$’ the magnetic and non-magnetic
contributions, one has for the four Stokes profiles, to which the algorithm is
applied

$\left\\{\begin{array}[]{l}I=(1-\alpha)I_{\mathrm{nm}}+\alpha
I_{\mathrm{m}}\\\ Q=\alpha Q_{\mathrm{m}}\\\ U=\alpha U_{\mathrm{m}}\\\
V=\alpha V_{\mathrm{m}}\end{array}\right.\;.$ (1)

The magnetic and non-magnetic contributions are computed from the Unno-
Rachkovsky solution with all their physical parameters taken as equal except
for the presence or absence of the three magnetic field coordinates. In other
words, the dependencies are
$I_{\mathrm{m}}(\lambda_{0},\Delta\lambda_{\mathrm{D}},\eta_{0},H,\theta_{H},\varphi_{H},a,B_{1}/B_{0})$
and
$I_{\mathrm{nm}}(\lambda_{0},\Delta\lambda_{\mathrm{D}},\eta_{0},a,B_{1}/B_{0})$,
where $(H,\theta_{H},\varphi_{H})$ are the three spherical coordinates of the
magnetic field vector: strength, inclination, and azimuth. In the present
work, the same weight of unity was given to the four Stokes parameters for the
chi-square calculation because $I$ plays an important role for the
determination of $\alpha$ (usually, $I$ is less weighted).

We tested the code and validated it (Bommier et al. 2007), but we obtained
that the unresolved magnetic field strength $H$ and the magnetic filling
factor $\alpha$ cannot be separately determined (see Fig. 4 of that paper).
$H$ is here the modulus of the magnetic field vector, which is different from
its longitudinal component. Only their product $\alpha H$, which is the local
average magnetic field strength, is 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 central depth, and three Stokes
parameters, which enable the calculation of the three magnetic field
components. This results in only 6 independent parameters. The Voigt $a$
parameters characterizes 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 stregth $\alpha H$ is finally well determined. Later
on, a tenth parameter was introduced in the inversion, which is the velocity
gradient along the line formation depth, which depends on an additional
observed parameter, which is the line asymmetry (Molodij et al. 2011).

The validation was achieved by creating test profiles from the Unno-Rachkovsky
solution, which we submitted to the inversion code. We thus compared the ouput
magnetic field values with the input ones. We analyzed 183,600 theoretical
profiles. The details of their definition are given in Bommier et al. (2007),
Sect. 2.2. The theoretical initial profiles were noised 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). This informed us also about the accuracy of the inversion method by
looking at the widths of the histograms of the difference between input and
output, which we plotted in Fig. 5 of Bommier et al. (2007). We thus obtained
accuracies of 15 G on the field strength (averaged with the magnetic filling
factor, $\alpha H$), 5∘ on the field inclination and azimuth for
$\alpha\mu_{0}H>45$ G, and 20 G on the field strength, 10∘ on the field
inclination, and 15∘ on the field azimuth for $\alpha\mu_{0}H<45$ G, for the
$1.5\times 10^{-3}$ level of polarimetric noise, for Fe I 6302.5 Å. For
6301.5, the values are similar for $\alpha\mu_{0}H>45$ G and a bit larger for
$\alpha\mu_{0}H<45$ G, due to the lower magnetic sensitivity, which results
from the lower effective Landé factor of this line (see Fig. 7 of the same
paper).

The inversion is performed in the four Stokes profiles simultaneously. It is
well-known that in weak fields the linear polarization Stokes parameters
$(Q,U)$ depend quadratically on the transverse magnetic field, whereas the
circular polarization Stokes parameter $V$ depends linearly 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 accuracies. We plotted the difference histograms for the cartesian
components $\mu_{0}H_{x},\mu_{0}H_{y},\mu_{0}H_{z}$ in the line-of-sight
reference frame assumed to be at disk-center for the same test data as above.
We obtained quite comparable widths or inaccuracy of 15 G for $\mu_{0}H_{x}$
and $\mu_{0}H_{y}$, and 10 G for $\mu_{0}H_{z}$, for $\alpha\mu_{0}H>10$ G.

Although the inversion finally provides only the $\alpha\mu_{0}H$ 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). We were thus able to
treat any kind of line.

##  3 Data Analysis Second Step: Azimuth Ambiguity Resolution

We applied a modification to the ME0 code developed by Metcalf, Crouch,
Barnes, & Leka and now available on the web111http://www.cora.nwra.com/AMBIG/
(Leka et al. 2009). This code applies the “Minimum Energy” Method initially
described by Metcalf (1994), which consists in searching for the field vector
orientation that minimizes
$\left|\mathop{\mathrm{d}iv}\@vec{H}\right|+\lambda\left|J_{Z}\right|$, where
$\@vec{J}$ is the current density vector, $J_{Z}$ its vertical component, and
$\lambda$ a positive weight parameter usually fixed at unity (with $\@vec{H}$
expressed in Gauss). Minimizing $\left|\mathop{\mathrm{d}iv}\@vec{B}\right|$
is a natural requirement imposed by Maxwell’s equations. As explained in the
Introduction, we will extent this requirement to what is really measured,
i.e., $\left|\mathop{\mathrm{d}iv}\@vec{H}\right|$. 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 in 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, back and forth transformations have to be performed
between these two systems of coordinates.

More precisely, $\mathop{\mathrm{d}iv}\@vec{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
$\mathop{\mathrm{d}iv}\@vec{H}$ in the line-of-sight coordinates instead,
which was possible from our 2-line observations. We accordingly modified the
corresponding subroutine of ME0. We calculated the expression of
$\mathop{\mathrm{d}iv}\@vec{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 by Crouch et al. (2009); Crouch (2013, 2015). We have first to
introduce the various reference frames entering the calculation.

###  3.1 Reference frames

Our approach is in agreement with Gary & Hagyard (1990) in the limit $P=0$. In
the following, we denote as $Oxyz$ the line-of-sight (l.o.s.) reference frame
and $OXYZ$ the heliographic one at the observed region location. The l.o.s.
reference frame is defined by $Oz$ being the l.o.s. itself oriented towards
the observer, and $Oy$ being 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 solar north oriented, and $OZ$ along the solar radius oriented from
the Sun’s center towards outside. We denote as $\@vec{R}$ the unit vector
along the $OZ$ axis of the heliographic reference frame, which is also the
local solar radius. We denote with indexes $(x,y,z)$ the vector coordinates in
the l.o.s. reference frame, and $(X,Y,Z)$ 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 $\@vec{R}$ coordinates are

$\left\\{\begin{array}[]{l}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{array}\right.\ ,$ (2)

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 $\@vec{I}$ along $OX$ and
$\@vec{K}$ along $OY$, with

$\left\\{\begin{array}[]{l}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{array}\right.$ (3)

and

$\left\\{\begin{array}[]{l}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{array}\right.$ (4)

in the l.o.s. reference frame. These formulae are obtained by applying 3
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
$\@vec{i}$, $\@vec{k}$ and $\@vec{\ell}$ (l.o.s. vector) in the heliographic
reference frame are

$\left\\{\begin{array}[]{l}i_{X}=\cos(L-L_{c})\\\ i_{Y}=-\sin(L-L_{c})\sin
b\\\ i_{Z}=\sin(L-L_{c})\cos b\end{array}\right.$ (5)

and

$\left\\{\begin{array}[]{l}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{array}\right.$ (6)

and

$\left\\{\begin{array}[]{l}\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{array}\right.\ .$ (7)

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

$\left\\{\begin{array}[]{l}X=xi_{X}+yk_{X}+z\ell_{X}\\\
Y=xi_{Y}+yk_{Y}+z\ell_{Y}\\\ Z=xi_{Z}+yk_{Z}+z\ell_{Z}\end{array}\right.$ (8)

and conversely

$\left\\{\begin{array}[]{l}x=XI_{x}+YK_{x}+ZR_{x}\\\ y=XI_{y}+YK_{y}+ZR_{y}\\\
z=XI_{z}+YK_{z}+ZR_{z}\end{array}\right.\ .$ (9)

###  3.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 $\Delta x\neq\Delta y$ in general, 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 meshs are no more rectangular in general,
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}$.

###  3.3 The quantity to be minimized

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

As stated above, we applied the “simulated annealing” algorithm to globally
minimize
$\left|\mathop{\mathrm{d}iv}\@vec{H}\right|+\lambda\left|J_{Z}\right|$. The
parameter $\lambda$ was taken as unity (with $\@vec{H}$ expressed in Gauss).
Indeed, in his pioneering paper, Metcalf (1994) recommends to minimize the
alternative quantity
$\left|\mathop{\mathrm{d}iv}\@vec{H}\right|+\lambda\left\|\@vec{J}\right\|$ 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 2-line observations, but it did not converge correctly. On
the contrary, we found that the algorithm sucessfully worked out, when we
instead minimized
$\left|\mathop{\mathrm{d}iv}\@vec{H}\right|+\lambda\left|J_{Z}\right|$ 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 articifial defaults
like 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
$\left|\mathop{\mathrm{d}iv}\@vec{H}\right|+\lambda\left|J_{Z}\right|$ was
later performed also by Metcalf’s team (Metcalf et al. 2006).

However, the map under study is usually not located at 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$. We denote as “map formation plane” the heliographic
horizontal plane, i.e., the plane perpendicular to the local solar radius. The
map results from interpretation of polarization in spectral lines. As
explained below in Sect. 3.4, the radiation received along the line-of-sight
in a given spectral line comes essentially from a certain height in the solar
atmosphere (this is the Eddington-Barbier approximation), which defines an
heliographic horizontal plane, which 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
disk center.

But the divergence has to be computed in a rectangular reference frame. For
doing this, 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. Or, as for us, 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 line-of-sight or $H^{(l)}$-type magnetic field vector
coordinates, as derived in Appendix A.

In Appendix A we derive the expressions of $\mathop{\mathrm{d}iv}\@vec{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 $\@vec{R}$ of l.o.s.
components $R_{x,y,z\text{ }}$given in Eq. (2). We obtain

$\displaystyle\mathop{\mathrm{d}iv}\@vec{H}$ $\displaystyle=$
$\displaystyle\frac{\Delta_{x}H_{x}^{(m)}}{\Delta
x}+\frac{\Delta_{y}H_{y}^{(m)}}{\Delta y}+\frac{\Delta_{z}H_{z}^{(m)}}{\Delta
z}$ $\displaystyle+\frac{R_{x}}{R_{z}}\frac{\Delta_{z}H_{x}^{(m)}}{\Delta
z}+\frac{R_{y}}{R_{z}}\frac{\Delta_{z}H_{y}^{(m)}}{\Delta z}\ .$

We denote as $\@vec{H}^{(m)}$ 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. $\Delta_{x}H_{x}^{(m)}$ and $\Delta_{y}H_{y}^{(m)}$ are
respectively the difference of the measured component $H_{x}^{(m)}$or
$H_{y}^{(m)}$ between two neighboring pixels in $x$ or $y$ direction. $\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 of
calculating 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. (3.3) simplifies in

$\displaystyle\mathop{\mathrm{d}iv}\@vec{H}$ $\displaystyle=$
$\displaystyle\frac{\Delta_{x}H_{x}^{(m)}}{\Delta
x}+\frac{\Delta_{y}H_{y}^{(m)}}{\Delta y}+\frac{\Delta_{z}H_{z}^{(m)}}{\Delta
z}$ $\displaystyle+\tan\theta\frac{\Delta_{z}H_{x}^{(m)}}{\Delta z}$

where $\theta$ is the heliocentric angle (positive for the West side of the
central meridian and negative for the East 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 one considers two neighboring pixels, referred to as
$P\left(i,j\right)$ and $P\left(i+1,j\right)$, separated by the length $\Delta
x$ in the plane perpendicular to the line-of-sight, and if one considers a
quantity $A$ measured via the same line at the two places, namely $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$
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
disk center, i.e., 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}$, one has $z\left[P(i+1,j)\right]\neq
z\left[P(i,j)\right],$ 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. This problem was
addressed in a different way by Crouch et al. (2009); Crouch (2013, 2015).
However, the final results are in agreement, when observation data are treated
as described below in Sect. 5.

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

$\displaystyle\Delta_{x}H_{x}^{(m)}$ $\displaystyle=$
$\displaystyle(H_{x}^{(1)}(i+1,j)-H_{x}^{(1)}(i,j)$
$\displaystyle+H_{x}^{(1)}(i+1,j+1)-H_{x}^{(1)}(i,j+1)$
$\displaystyle+H_{x}^{(2)}(i+1,j)-H_{x}^{(2)}(i,j)$
$\displaystyle+H_{x}^{(2)}(i+1,j+1)-H_{x}^{(2)}(i,j+1))/4\ ,$
$\displaystyle\Delta_{y}H_{y}^{(m)}$ $\displaystyle=$
$\displaystyle(H_{y}^{(1)}(i,j+1)-H_{y}^{(1)}(i,j)$
$\displaystyle+H_{y}^{(1)}(i+1,j+1)-H_{y}^{(1)}(i+1,j)$
$\displaystyle+H_{y}^{(2)}(i,j+1)-H_{y}^{(2)}(i,j)$
$\displaystyle+H_{y}^{(2)}(i+1,j+1)-H_{y}^{(2)}(i+1,j))/4\ ,$
$\displaystyle\Delta_{z}H_{z}^{(m)}$ $\displaystyle=$
$\displaystyle(H_{z}^{(2)}(i,j)-H_{z}^{(1)}(i,j)$
$\displaystyle+H_{z}^{(2)}(i+1,j)-H_{z}^{(1)}(i+1,j)$
$\displaystyle+H_{z}^{(2)}(i,j+1)-H_{z}^{(1)}(i,j+1)$
$\displaystyle+H_{z}^{(2)}(i+1,j+1)-H_{z}^{(1)}(i+1,j+1))/4\ ,$

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)

$\displaystyle J_{Z}$ $\displaystyle=$ $\displaystyle
R_{z}\left[\frac{\Delta_{x}H_{y}^{(m)}}{\Delta
x}-\frac{\Delta_{y}H_{x}^{(m)}}{\Delta y}\right]$
$\displaystyle+R_{x}\frac{\Delta_{y}H_{z}^{(m)}}{\Delta
y}-R_{y}\frac{\Delta_{x}H_{z}^{(m)}}{\Delta x}\ .$

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)

$\displaystyle J_{x}$ $\displaystyle=$
$\displaystyle\left[\frac{\Delta_{y}H_{z}^{(m)}}{\Delta
y}-\frac{\Delta_{z}H_{y}^{(m)}}{\Delta
z}\right]+\frac{R_{y}}{R_{z}}\frac{\Delta_{z}H_{z}^{(m)}}{\Delta z}\ ,$ (16)
$\displaystyle J_{y}$ $\displaystyle=$
$\displaystyle\left[\frac{\Delta_{z}H_{x}^{(m)}}{\Delta
z}-\frac{\Delta_{x}H_{z}^{(m)}}{\Delta
x}\right]-\frac{R_{x}}{R_{z}}\frac{\Delta_{z}H_{z}^{(m)}}{\Delta z}\ .$ (17)

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

$J_{z}=\left(J_{Z}-J_{x}R_{x}-J_{y}R_{y}\right)/R_{z}\ ,$ (18)

which can be reduced into

$\displaystyle J_{z}$ $\displaystyle=$
$\displaystyle\left[\frac{\Delta_{x}H_{y}^{(m)}}{\Delta
x}-\frac{\Delta_{y}H_{x}^{(m)}}{\Delta y}\right]$
$\displaystyle+\frac{R_{x}}{R_{z}}\frac{\Delta_{z}H_{y}^{(m)}}{\Delta
z}-\frac{R_{y}}{R_{z}}\frac{\Delta_{z}H_{x}^{(m)}}{\Delta z}\ .$

The second line of the formula accounts for the inclination of the line
formation plane with respect to the line-of-sight. $J_{X}$ and $J_{Y}$ can
then be derived by applying Eq. (8) to $J_{x},J_{y}$ and $J_{z}$.

In the following, we have plotted the heliographic reference frame components
$J_{X,Y,Z}$ 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. (8).

Crouch et al. (2009); 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. But their formula given in Eqs. (4-5) of Crouch et al. (2009);
Crouch (2013, 2015) is different from our Eq. (3.3) because the corresponding
reference frames are not the same. The divergence of Crouch et al. (2009);
Crouch (2013, 2015) applies 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 applies spatial
derivatives with respect to all the line-of-sight $(x,y,z)$ reference frame
coordinates, including the fact that the observed lines are formed along
inclined planes. A possible method of comparison is to run codes based on each
formula on the same data. This is done in Sect. 5, where it can be seen that
the results obtained by both methods are in excellent agreement in similar
application conditions. Both approaches are different ways to treat the
problem of the inclined line-of-sight.

###  3.4 The two line formation depth difference

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

The quantity to be determined is the difference in 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 behave in a very parallel manner.
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
why we discarded 5250.2, though it is more sensitive to the magnetic field
than 6301.5. Such a parallelism probably originates in the fact that the two
lines 6302.5 and 6301.5 belong to the same multiplet n. 816 and have different
$gf$ values, which 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 3 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 Na I D1
and D2 lines form another such pair. They obey $gf=1.309$ for Na I D2 and
$gf=0.655$ for Na I D1 (from the Kurucz data-basis), which are in a factor of
2 ratio. Thus Na I D2 forms higher than Na I D1.

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), who obtained the observed value $63.2\pm 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
Uitenbroek’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 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 line center, this part
of the profile is the most contributing to the field vector determination, so
that it it is quite natural to assign the line center formation height to this
field value. For instance, it is found that the value of $1/2\arctan U/Q$
computed at 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 well 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 has an accuracy of the order of $\pm$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. Fig. 4 of Khomenko & Collados (2007) also shows
that the formation depth difference remains on the order of $\Delta z=$ 70 km
in the quiet sun, whatever the formation depth is in itself.

Finally, in Fig. 4 of Khomenko & Collados (2007) we observed that the depth
difference is larger in active regions than in quiet ones. We thus applied
instead the depth difference of 98 km at 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.

As for the Na I D lines, we derived their 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 IAC). The continuum
absorption coefficient was evaluated as in the MALIP code of Landi
Degl’Innocenti (1976), i.e., 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), i.e., where $\tau/\mu=1$, where $\tau$ is the line center
optical depth along the vertical and $\mu$ the cosine of the heliocentric
angle $\theta$ (here supposed to be 0). As it can be 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 533 km for Na I D1
and 604 km for Na I D2 above the $\tau_{5000}=1$ level, which are located in
the low chromosphere close to the temperature minimum and with a difference of
71 km. We investigated also the Fe I 6302.5 and 6301.5 quiet sun line center
formation heights. We obtained 262 km for 6302.5 and 328 km for 6301.5,
leading to a difference of 66 km in excellent agreement with the measurements
described above, which validates our computation method.

###  3.5 The aspect ratio of the strongly stratified atmosphere

We were indeed confronted with the problem that the vertical gradient of the
magnetic field $\partial H_{z}/\partial z$ was found on the order of 3 G/km in
spot umbrae, whereas the horizontal gradient $\partial H_{x}/\partial
x+\partial H_{y}/\partial y$ was only on the order of 0.3 G/km. This leads to
a non-vanishing value of $\mathop{\mathrm{d}iv}\@vec{H}$, whatever the signs
would be. We then investigated the literature and found 15 references, which
fully confirm these values provided by different instruments (groundbased as
well as spaceborn), different inversion mesthods (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) and Bommier (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 $\mathop{\mathrm{d}iv}\@vec{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 be also zero, within the
noise level. The question arised to know if the observed value of
$\mathop{\mathrm{d}iv}\@vec{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 a filtering. The mathematical demonstration
that the filtered divergence is the divergence of the filtered quantity is
given in Bommier (2013) and Bommier (2014). Bommier (2014) includes also an
easier 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 to a zero divergence computed
by finite differences, when the local divergence is zero. To 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 $\mathop{\mathrm{d}iv}\@vec{H}$ to be
zero. We argue that the existence of magnetic monopoles is not the only
possibility for a non-zero $\mathop{\mathrm{d}iv}\@vec{H}$. After a discussion
about the measurement noise level effect, we present below another possibility
we 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 disk center. This results in an inaccuracy of 0.2 G/km for
the vertical field gradient observed at 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 disk center. The total inaccuracy on the observed
$\mathop{\mathrm{d}iv}\@vec{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
$\mathop{\mathrm{d}iv}\@vec{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. 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 20 (Bommier
2013, 2014) in the quiet sun photosphere. We proposed a physical process able
to explain how the plasma strong stratification influences the magnetic field.
The magnetic field is created by moving charges and we showed that the Debye
screening applies also to moving charges. The anisotropy of the screening
finally leads to possible non-vanishing of $\mathop{\mathrm{d}iv}\@vec{H}$.
For a usual, spherical, Debye sphere, one has
$\mathop{\mathrm{d}iv}\@vec{H}=0$ (within negligible magnetization). The
strong stratification is responsible for anistropy of the velocities, which
results in a flattened Debye sphere and then to the possibility of
$\mathop{\mathrm{d}iv}\@vec{H}\neq 0$. However, orders of magnitude concluded
to a too weak resulting $\mathop{\mathrm{d}iv}\@vec{H}$ with respect to the
observations. We then turned to the interior electron thermal escape (Bommier
2020) as the final explanation of the observations. Nevertheless, we finally
showed in our observations that once the inverse of the aspect ratio is
applied to scale the different magnetic field components, the scaled
$\mathop{\mathrm{d}iv}\@vec{H}$ vanishes, which enables the ambiguity
resolution. An example of this result in the case of NOAA 10808 observed with
THÉMIS on 13 September 2005, is visible 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 different in sunpots
or plages or active regions, which is the case of the regions we treated. This
is not so well-known and sometimes we adjusted the ratio (kept constant in the
whole map) until the ambiguity resolution is consistent in the whole map. With
the word “consistent” here, we mean 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 sunpot
penumbrae as discussed in the following subsection. In Sect. 5.2 we obtain the
same result by applying the inverse scaling before submitting the observation
data to the available ambiguity resolution code AMBIG2.

It is also reported in Bommier (2015) that in the case of the photosphere
flattened Debye sphere the current density remains given by
$\mathop{\mathrm{c}url}\@vec{H}=\@vec{J}$. We applied this law for the current
density derivation from the magnetic field vector in different regions
reported in the following.

###  3.6 Verification of the ambiguity resolution results

We thus 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
regions222http://lesia.obspm.fr/perso/veronique-bommier/. 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 is 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 as 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.

###  3.7 Test of the method on theoretical data

Figure 3: Result of test of the ambiguity resolution by 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. (3.3) 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. (3.3).

Our method for taking the departure from disk center into account, which leads
to additional terms in the divergence formula Eq. (3.3), 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 2-line
observations instead from extrapolation as described at the beginning of Sect.
3. The 2-line observation had been simulated with theoretical results at two
different depths separated by 70km along the line-of-sight. The pixel sizes
are 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∘ from the central meridian in longitude
and at 45∘ 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.
(3.3) only, is displayed in the middle map. The minimized quantity is in fact
$\left|\mathop{\mathrm{d}iv}\@vec{H}\right|+\lambda\left|J_{Z}\right|$ with
$\lambda=1$ as in our method. It can be seen that the amiguity resolution
widely fails over the $\delta$-spot. When the full Eq. (3.3) is applied
instead, the obtained result is displayed in the bottom map of Fig. 3, which
is in total agreement with the reference field. It is thus confirmed by
numerical test that the second line of Eq. (3.3) is necessary to resolve the
ambiguity by applying a two-line analysis to maps observed out of disk center
as discussed after Eq. (3.3) and in Fig. 1.

###  3.8 Test of the resolution code on theoretical data

pixel | line | $\mathcal{M}_{\mathrm{area}}$ | $\mathcal{M}_{H_{\bot}>100\,\mathrm{G}}$ | $\mathcal{M}_{H_{\bot}>500\,\mathrm{G}}$
---|---|---|---|---
size |  |  |  |
0.3” | lower | 0.99 | 1.00 | 1.00
0.3” | upper | 0.99 | 1.00 | 1.00
0.9” | lower | 0.99 | 1.00 | 1.00
0.9” | upper | 0.99 | 1.00 | 1.00
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). The data were provided and the metrics were computed by K.D. Leka’s courtesy. noise | $\mathcal{M}_{\mathrm{area}}$ | $\mathcal{M}_{H_{\bot}>100\,\mathrm{G}}$ | $\mathcal{M}_{H_{\bot}>500\,\mathrm{G}}$
---|---|---|---
level |  |  |
no | 1.00 | 1.00 | 1.00
low | 0.85 | 0.95 | 1.00
high | 0.73 | 0.83 | 0.90
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. 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.

We tested our ambiguity resolution code DIVB2 built on the above described
method, 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); Crouch (2013, 2015), where
magnetic field data at two different heights are also used, were also
submitted to these test data. We used two series of these data.

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 in p. 93 of Leka et al.
(2009) and also in 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 disk center,
i.e., 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 $\mathcal{M}_{\mathrm{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 visible 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 in 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. (3.3) 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. (3.3). The success of our resolution of theoretical
data again validates our Eq. (3.3) and, in particular, its second line in the
case of lines formed along horizontal planes and observed out of disk center.

##  4 Results: Examples of Vector Maps

Figure 4: Active Region NOAA 11420 observed by HINODE/SOT/SP on 18 February
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/m2. Right: Lorentz force vector, vertical
component between $-448$ and $+448$ , horizontal component typical arrow
length 800 mN/m2. The spatial resolution was reduced by a factor $5\times 5$
for paper file size purposes. (color figure on-line)

Figure 5: Active Region NOAA 11494 observed by HINODE/SOT/SP on 7 June 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/m2. Right: Lorentz force vector, vertical
component between $-682$ and $+682$ , horizontal component typical arrow
length 800 mN/m2. The spatial resolution was reduced by a factor $5\times 5$
for paper file size purposes. (color figure on-line)

Figure 6: Active Region NOAA 11476 observed by HINODE/SOT/SP on 12 May 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$, horizontal component
typical arrow length 8000 mA/m2. Right: Lorentz force vector, vertical
component between $-1040$ and $+1040$ , horizontal component typical arrow
length 800 mN/m2. The spatial resolution was reduced by a factor $5\times 5$
for paper file size purposes. (color figure on-line)

Figure 7: Active Region NOAA 11429 observed by HINODE/SOT/SP on 6 March 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/m2. Right: Lorentz force vector, vertical
component between $-1180$ and $+1180$ , horizontal component typical arrow
length 800 mN/m2. The spatial resolution was reduced by a factor $5\times 5$
for paper file size purposes. (color figure on-line)

Figure 8: Active Region NOAA 11857 observed with THEMIS on 7 October 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$, horizontal component maximum arrow
length 1070 G. Middle: current density vector, vertical component between
$-68$ and $+68$, horizontal component maximum arrow length 746 mA/m2. Right:
Lorentz force vector, vertical component between $-56$ and $+56$ , horizontal
component maximum arrow length 133 mN/m2. (color figure on-line)

Figure 9: Active Region NOAA 11857 observed with THEMIS on 7 October 2013
between 09:24and 10:42 UT, in the chromospheric lines Na I D1/D2. The
ambiguity is resolved by our method. Left: magnetic field vector, vertical
component between $-772$ and $+772$, horizontal component maximum arrow length
1254 G. Middle: current density vector, vertical component between $-188$ and
$+188$, horizontal component maximum arrow length 1495 mA/m2. Right: Lorentz
force vector, vertical component between $-155$ and $+155$ , horizontal
component maximum arrow length 54 mN/m2. (color figure on-line)

Figure 10: Active Region NOAA 11865 observed with THEMIS on 11 October 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$, horizontal component maximum arrow
length 1320 G. Middle: current density vector, vertical component between
$-96$ and $+96$, horizontal component maximum arrow length 655 mA/m2. Right:
Lorentz force vector, vertical component between $-83$ and $+83$ , horizontal
component maximum arrow length 43 mN/m2. (color figure on-line)

Figure 11: Active Region NOAA 11865 observed with THEMIS on 7 October 2013
between 10:14 and 11:19 UT, in the chromospheric lines Na I D1/D2. The
ambiguity is resolved by our method. Left: magnetic field vector, vertical
component between $-930$ and $+930$, horizontal component maximum arrow length
979 G. Middle: current density vector, vertical component between $-148$ and
$+148$, horizontal component maximum arrow length 1127 mA/m2. Right: Lorentz
force vector, vertical component between $-79$ and $+79$ , horizontal
component maximum arrow length 59 mN/m2. (color figure on-line)

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 page333http://lesia.obspm.fr/perso/veronique-bommier. 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-11 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 $\Delta
x=862$ km and $\Delta y=928$ km. This was also done for the tests described in
the following Sect. 5. The THÉMIS data pixel size was $\Delta x=794$ km and
$\Delta y=725$ km, the binning was not performed.

###  4.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 anticlockwise
about a negative polarity spot. As a consequence, the Lorentz force vector is
systematically found 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 18 February 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 7 June 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
12 December 2006, 7 September 2011, 1 February 2012, 1 May 2012, and 13 June
2012.

###  4.2 HINODE data: examples of a neutral line

We often observe a strong horizontal current component, which crossesthe
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 12 May 2012
between 02:38 and 03:11 UT, which is represented in Fig. 6, and NOAA 11429
observed by HINODE/SOT/SP on 6 March 2012 between 22:10 and 22:43, which is
represented in Fig. 7. One 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 7 September 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.

###  4.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 7 October 2013, between 08:57 and 09:22 UT in Fe I 6301/6302,
between 09:24 and 09:59 UT in Na I D2 and between 10:10 and 10:42 UT in Na I
D1. In the case of the photosphere, which is observed in Fe I 6301/6302 and is
represented in Fig. 8, we find again the current and Lorentz force
characteristics described above, namely circular current about the spot and
centripetal Lorentz force. Higher, at the low chromosphere level, which is
observed in Na I D1/D2 and is represented in Fig. 9, the field and current
organizations are found more complex. We observe, however, a good agreement
between the photospheric and chromospheric magnetic fields, which leads us to
think that the ambiguity was correctly resolved at the chromospheric level
observed with the Na I D lines. NOAA 11857 was observed with THÉMIS in both
iron and sodium lines on 5 October, 7 October (two times, morning and
afternoon), 8 October, and 9 October. 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.

###  4.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 11 October 2013, between 09:39 and 09:59 UT in Fe I 6301/6302,
between 10:14 and 10:43 UT in Na I D2, and between 10:50 and 11:19 UT in Na I
D1. Again, in the photosphere, which is observed in Fe I 6301/6302 and is
represented in Fig. 10, we observe a strong horizontal current component
across the neutral line. The structure of the currents is more involved at the
low chromosphere level, which is observed in Na I D1/D2 and is represented in
figure 11.

##  5 Verification: comparison with results obtained by applying available
codes

Figure 12: Current vector maps of the four active regions of Figs. 4-7, but
the ambiguity is resolved by another method. The ambiguity is resolved
independently for each line map by applying the Metcalf ME0 code (Leka et al.
2009). The current vector is then plotted by assuming the depth difference 66
km (see Sect. 3.4), corrected for the disk center departure, between the two
line maps.

Figure 13: Same as Fig. 4, but the ambiguity is resolved by applying the two
line code by Crouch (2013) in the conditions corresponding to Eq. (8) of that
paper. The ambiguity resolution solution is not in agreement with the spot
polarity in a zone South-West (i.e., bottom-right) of the spot center. (color
figure on-line)

Figure 14: Same as Fig. 13, but our proposed scaling has been applied without
modifying the ambiguity resolution code, by artificially dividing by 10 the
$x$ and $y$ pixel size in the image plane. The ambiguity resolution solution
is now in agreement with the spot polarity in the whole map. (color figure on-
line)

In order to check our method, we applied different available codes for
ambiguity resolution to our data.

###  5.1 Case of the 1-line Metcalf ME0 code applied to each line

First, we applied the Metcalf ME0 code described by Leka et al. (2009) to each
line map. This code is able to resolve the ambiguity by using only one line
data. For computing the magnetic field divergence for minimization, the
horizontal derivatives are computed from the data, whereas the vertical
derivative is derived from a current-free extrapolation from the data. The
departure from disk center is accounted for as explained in Crouch et al.
(2009). This code does not apply the scaling by the aspect ratio that we
propose in Sect. 3.5 and we did not modify this point. This is not necessary
because the extrapolation step is performed within the divergence-free
condition. We did not at all modify the code.

We applied the code independently to each of the two Fe I 6301.5 Å and 6302.5
Å lines of our HINODE/SOT/SP observations of Figs. 4-7. We then plotted the
current density vector by computing it from the two line magnetic field maps,
by assuming a depth difference of 66 km between the two lines following Sect.
3.4. The results are given in Fig. 12, where it can be seen that the direction
of the current density vector well matches the one we obtained with our method
and given in Figs. 4-7. In particular, the direction of the spot wrapping
current is fully recovered, together with the way how the current density
vector crosses the neutral line. We then infer that our ambiguity resolution
was correct. However, it can be seen that the result is more noisy with the
ME0 resolution method thus applied.

###  5.2 Case of the two-line AMBIG2 code proposed by Crouch (2013)

We then applied the two-line code AMBIG2 proposed by Crouch
(2013)444http://www.cora.nwra. com/~ash/ambig2.tar.gz, who presented runs on
theoretical data. We applied it to our real data represented in Figs. 4-7. For
a better comparison, we ran the code in the conditions of Eq. (8) of Crouch
(2013), which are the closest conditions to those we applied as described in
Sect. 3.3. As in our resolution code, we assigned a value of 98 km to the two
line depth difference, which is a bit higher than the 66 km value described
above in Sect. 3.4 because this last value was obtained in quiet Sun
conditions, whereas our data concern a sunspot. It can be seen in Figure 4 of
Khomenko & Collados (2007) that the two Fe I 6301.5 Å and 6302.5 Å formation
depth difference is a bit higher in active regions. The depth difference value
was corrected from the disk center departure.

As mentioned in Sect. 3.6, the ambiguity solution is intuitively known in
sunspots from their polarity. This enables a verification procedure. We
applied the code to our HINODE/SOT/SP data of Fig. 4. The results for
comparison are displayed in Fig. 13, where it can be seen that the ambiguity
solution does not correspond to the spot polarity inside a zone located South-
West (i.e., bottom-right) from the spot center.

For comparison, we then simulated our aspect ratio scaling described in Sect.
3.5 without modifying the code, by artificially dividing by 10 the $x$ and $y$
pixel sizes. The spot was located at longitude $-8{{}^{\circ}}$ from the
central meridian and latitude $23{{}^{\circ}}$. As a result, the heliospheric
angle was $24{{}^{\circ}}$ with cosine $\mu=0.91$. The spot was not so far
from disk center. Therefore, the above described simulation of the aspect
ratio scaling was sufficient without any other modification and in particular
not any code modification. The results are plotted in Fig. 14, where it can be
seen that now the ambiguity solution well matches the spot polarity in the
whole map. The success of the ambiguity resolution when the scaling by the
aspect ratio is applied as in Fig. 14, whereas the resolution is unsuccessful
when the scaling by the aspect ratio is not applied as in Fig. 13, confirms
our conviction about the necessity of applying this aspect ratio scaling to
vanish the magnetic field divergence in the data following our analysis in
Bommier (2013, 2014) and recalled in Sect. 3.5.

The current density vector directions obtained are also in full agreement with
the directions we obtained with our method. This is in particular the case of
the direction of the spot wrapping current. This also validates our ambiguity
resolution.

The same analysis was applied with the same results to the data corresponding
to Figs. 5-7. No map was successfully resolved without the scaling. For the
data of Fig. 5 and Fig. 7, the ambiguity resolution was successful when the
applied scaling factor was 20 instead of 10, which exactly matches the
theoretical aspect ratio we derived for the photosphere in Bommier (2013,
2014). For the data of Fig. 6, it was not possible to get a correct ambiguity
resolution whatever the applied scaling factor is. These analyses on several
cases confirm our intuition that the theoretical value of 20 for the scaling
factor has to be considered more as an order of magnitude than as a well-
defined value, as previously stated in Sect. 3.5.

One of the main results of the present paper is the necessity to apply the
scaling factor to the data before disambiguation. This can be done in an
external way by artificially dividing the pixel sizes before sumbission to the
AMBIG2 code, provided that the observed region is located close to disk
center. Our disambiguation code is able to also disambiguate data of regions
located out of disk center because departure from disk center is fully
integrated into Eq. (3.3).

##  6 Conclusion

We have shown that observations in two different lines, which belong to the
same multiplet but have different absorption coefficients so that they are
formed at two different depths, like Fe I 6302.5 Å and 6301.5 Å, enable the
resolution of the azimuth ambiguity, which remain from the Zeeman signal
interpretation. The anisotropy of the strongly stratified plasma of the
photosphere has also to be accounted for following Bommier (2013) and Bommier
(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 $\mathop{\mathrm{r}ot}\@vec{H}=\@vec{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. (2010) from HINODE/SOT/SP data. We observe some systematical
trends, like circular currents wrapping spots clockwise about a positive
polarity spot and anticlockwise 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. (2010), who applies completely
different methods. As a result, the Lorentz force may be computed. It is found
to be centripetal in sunspots.

We successfully compared the results of our method to the results of other
codes based on other methods and applied to the same observation data. In
particular, we applied the publicly available AMBIG2 code by Crouch (2013) to
our real solar data. We obtained 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) and Bommier (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.

We finally remarked that the Na I D1 and D2 lines form another such line pair.
We successfully tested their application to full vector mapping of the
magnetic field and current density with THÉMIS telescope observation.
Observing this line pair would open access to the low chromosphere where they
are formed.

###### Acknowledgements.

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 HINODE/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 (2018) Balthasar, H. 2018, Sol. Phys., 293, 120
  * Bommier (2013) Bommier, V. 2013, Physics Research International, 2013, 195403
  * Bommier (2014) Bommier, V. 2014, Comptes Rendus Physique, 15, 430
  * Bommier (2015) Bommier, V. 2015, in IAU Symposium, Vol. 305, Polarimetry, ed. K. N. Nagendra, S. Bagnulo, R. Centeno, & M. Jesús Martínez González, 28–34
  * Bommier (2020) Bommier, V. 2020, A&A, 634, A40
  * Bommier et al. (2007) Bommier, V., Landi Degl’Innocenti, E., Landolfi, M., & Molodij, G. 2007, A&A, 464, 323
  * Bommier et al. (1991) Bommier, V., Landi Degl’Innocenti, E., & Sahal-Bréchot, S. 1991, A&A, 244, 383
  * Brethouwer et al. (2007) Brethouwer, G., Billant, P., Lindborg, E., & Chomaz, J. M. 2007, Journal of Fluid Mechanics, 585, 343
  * Bruls et al. (1991) 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–456
  * Crouch (2013) Crouch, A. D. 2013, Sol. Phys., 282, 107
  * Crouch (2015) Crouch, A. D. 2015, Sol. Phys., 290, 2677
  * Crouch et al. (2009) Crouch, A. D., Barnes, G., & Leka, K. D. 2009, Sol. Phys., 260, 271
  * Fan & Gibson (2004) Fan, Y. & Gibson, S. E. 2004, ApJ, 609, 1123
  * Faurobert et al. (2009) Faurobert, M., Aime, C., Périni, C., et al. 2009, A&A, 507, L29
  * Gary & Hagyard (1990) Gary, G. A. & Hagyard, M. J. 1990, Sol. Phys., 126, 21
  * Grec et al. (2010) Grec, C., Uitenbroek, H., Faurobert, M., & Aime, C. 2010, A&A, 514, A91
  * Khomenko & Collados (2007) Khomenko, E. & Collados, M. 2007, ApJ, 659, 1726
  * Landi Degl’Innocenti (1976) Landi Degl’Innocenti, E. 1976, A&AS, 25, 379
  * Landi Degl’Innocenti et al. (1990) Landi Degl’Innocenti, E., Bommier, V., & Sahal-Bréchot, S. 1990, A&A, 235, 459
  * Landolfi & Landi Degl’Innocenti (1982) Landolfi, M. & Landi Degl’Innocenti, E. 1982, Sol. Phys., 78, 355
  * Landolfi et al. (1984) Landolfi, M., Landi Degl’Innocenti, E., & Arena, P. 1984, Sol. Phys., 93, 269
  * Leka et al. (2009) Leka, K. D., Barnes, G., Crouch, A. D., et al. 2009, Sol. Phys., 260, 83
  * Leka et al. (2005) Leka, K. D., Fan, Y., & Barnes, G. 2005, ApJ, 626, 1091
  * Leka et al. (2022) Leka, K. D., Wagner, E. L., Griñón-Marín, A. B., Bommier, V., & Higgins, R. E. L. 2022, Sol. Phys., 297, 121
  * Maltby et al. (1986) Maltby, P., Avrett, E. H., Carlsson, M., et al. 1986, ApJ, 306, 284
  * Metcalf (1994) Metcalf, T. R. 1994, Sol. Phys., 155, 235
  * Metcalf et al. (2006) Metcalf, T. R., Leka, K. D., Barnes, G., et al. 2006, Sol. Phys., 237, 267
  * Molodij et al. (2011) Molodij, G., Bommier, V., & Rayrole, J. 2011, A&A, 531, A139
  * Puschmann et al. (2010) Puschmann, K. G., Ruiz Cobo, B., & Martínez Pillet, V. 2010, ApJ, 721, L58
  * Ruiz Cobo & del Toro Iniesta (1992) Ruiz Cobo, B. & del Toro Iniesta, J. C. 1992, ApJ, 398, 375
  * Socas-Navarro (2005) Socas-Navarro, H. 2005, ApJ, 633, L57
  * Stenflo (1973) Stenflo, J. O. 1973, Sol. Phys., 32, 41
  * Vernazza et al. (1973) 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 $\@vec{R}$ in $O$ is perpendicular to the the line formation
plane, which is tangent to the solar surface. The $\@vec{R}$ coordinates in
the l.o.s. reference frame are given in Eq. (2). In the l.o.s. reference
frame, the line formation plane, which is perpendicular to $\@vec{R}$ of
coordinates $(R_{x},R_{y},R_{z})$, has then for general equation

$R_{x}x+R_{y}y+R_{z}z=0\ .$ (20)

Let’s denote as $\@vec{H}^{(1)}$ the magnetic field vector measured as a
function of $x$ and $y$ along the formation plane of the line number (1).
$\@vec{H}^{(1)}$ depends on the following variables

$\@vec{H}^{(1)}(x,y,z;R_{x}x+R_{y}y+R_{z}z=0)\ ,$ (21)

where the semicolon means “such as”. A similar law holds for the magnetic
field vector $\@vec{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)\ ,$ (22)

and analogously

$\Delta_{y}H_{y}^{(1)}=H_{y}^{(1)}(i,j+1)-H_{y}^{(1)}(i,j)\ .$ (23)

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)\ .$ (24)

###  A.1 Calculation of the divergence

The mathematical expression of $\mathop{\mathrm{d}iv}\@vec{H}$ is frame-
independent

$\mathop{\mathrm{d}iv}\@vec{H}=\frac{\partial H_{x}}{\partial
x}+\frac{\partial H_{y}}{\partial y}+\frac{\partial H_{z}}{\partial z}\ .$
(25)

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

The present Appendix is devoted to account 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 apply
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 consider
first 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

Figure 15: 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.

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. 15. 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. 15, $\@vec{S}$,
$\@vec{S}^{\prime}$, and $\@vec{S}^{\prime\prime}$. One has to calculate the
flux of the magnetic field through each elementary surface. Consider first the
surface element associated to the $\@vec{S}$ vector. This surface element is
part of the formation plane of line number (1), whereas $\@vec{S}^{\prime}$
and $\@vec{S}^{\prime\prime}$ join the two line formation planes. In the
l.o.s. reference frame $Oxyz$, the $\@vec{S}$ unit vector has the following
coordinates

$\@vec{S}\left\\{\begin{array}[]{l}S_{x}=\sin\theta\\\ S_{y}=0\\\
S_{z}=\cos\theta\end{array}\right.\ .$ (26)

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

$\@vec{H}^{(1)}\cdot\@vec{S}=H_{z}^{(1)}\cos\theta+H_{x}^{(1)}\sin\theta\ .$
(27)

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

$\Delta\ell=\frac{\Delta x}{\cos\theta}\ .$ (28)

The flux through this elementary surface is then finally

$\@vec{H}^{(1)}\cdot\@vec{S}\ \Delta\ell\Delta
y=(H_{z}^{(1)}+H_{x}^{(1)}\tan\theta)\ \Delta x\Delta y\ .$ (29)

The flux through the elementary surface associated to $\@vec{S}^{\prime}$ is

$\@vec{H}^{(m)}\cdot\@vec{S}^{\prime}\ \Delta z\Delta y=H_{x}^{(m)}\Delta
z\Delta y\ ,$ (30)

where the superscript $(m)$ denotes the average between the two lines $(1)$
and $(2)$. Finally, the flux through the third elementary surface
$\@vec{S}^{\prime\prime}$ is

$\@vec{H}^{(m)}\cdot\@vec{S}^{\prime\prime}\ \Delta\ell\Delta
z\cos\theta=H_{y}^{(m)}\Delta x\Delta z\ .$ (31)

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

$\displaystyle\mathop{\mathrm{d}iv}\@vec{H}\ \Delta x\Delta y\Delta z$
$\displaystyle=$
$\displaystyle\left[(H_{z}^{(1)}-H_{z}^{(2)})+(H_{x}^{(1)}-H_{x}^{(2)})\tan\theta\right]\
\Delta x\Delta y$ $\displaystyle+\left[H_{x}^{(m)}(x+\Delta
x,y)-H_{x}^{(m)}(x,y)\right]\ \Delta z\Delta y$
$\displaystyle+\left[H_{y}^{(m)}(x,y+\Delta y)-H_{y}^{(m)}(x,y)\right]\ \Delta
x\Delta z\ ,$

which results in Eq. (3.3) 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

Figure 16: Reference frame for the divergence evaluation: general case.

One considers now a point $O$ on the solar surface with longitude $L$ and
latitude $b$. The solar radius in $O$ denoted as $\@vec{R}$ is perpendicular
to the line formation plane. The $\@vec{R}$ coordinates in terms of $(L,b)$
are given in Eq. (2). There is a common perpendicular to this radius $R$ and
to the l.o.s. $Oz$. We denote as $Oy^{\prime}=OY$ this common perpendicular,
which is represented in the left part of Fig. 16. 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^{\prime}$ is
also different from the $Oy$ axis of the l.o.s. reference frame, therefore we
use a different notation. In the $Ox^{\prime}y^{\prime}z$ reference frame the
geometry is the same as in the equatorial case studied above, so that its
final result may be applied

$\displaystyle\mathop{\mathrm{d}iv}\@vec{H}$ $\displaystyle=$
$\displaystyle\frac{\Delta_{x^{\prime}}H_{x^{\prime}}^{(m)}}{\Delta
x^{\prime}}+\frac{\Delta_{y^{\prime}}H_{y^{\prime}}^{(m)}}{\Delta
y^{\prime}}+\frac{\Delta_{z}H_{z}^{(m)}}{\Delta z}$
$\displaystyle+\tan\theta\frac{\Delta_{z}H_{x^{\prime}}^{(m)}}{\Delta z}\ ,$

where the $\Delta$’s 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. 16. In this reference frame, the
$\@vec{i}_{x^{\prime}}$ unit vector along $Ox^{\prime}$ has the following
coordinates

$\@vec{i}_{x^{\prime}}\left\\{\begin{array}[]{l}\cos\varphi\\\ \sin\varphi\\\
0\end{array}\right.\ .$ (34)

Accordingly,

$H_{x^{\prime}}^{(m)}=H_{x}^{(m)}\cos\varphi+H_{y}^{(m)}\sin\varphi\ .$ (35)

One has by rotational invariance

$\frac{\Delta_{x^{\prime}}H_{x^{\prime}}^{(m)}}{\Delta
x^{\prime}}+\frac{\Delta_{y^{\prime}}H_{y^{\prime}}^{(m)}}{\Delta
y^{\prime}}=\frac{\Delta_{x}H_{x}^{(m)}}{\Delta
x}+\frac{\Delta_{y}H_{y}^{(m)}}{\Delta y}\ ,$ (36)

so that finally

$\displaystyle\mathop{\mathrm{d}iv}\@vec{H}$ $\displaystyle=$
$\displaystyle\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}$
$\displaystyle\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]\ .$

The $\@vec{R}$ vector coordinates can also be expressed as

$\@vec{R}\left\\{\begin{array}[]{l}\sin\theta\cos\varphi\\\
\sin\theta\sin\varphi\\\ \cos\theta\end{array}\right.\ .$ (38)

Accordingly, the above Eq. (A.1.2) simply results in Eq. (3.3) 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

Figure 17: 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 $\@vec{R}$.

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

$\@vec{i}_{p}\left\\{\begin{array}[]{c}1\\\ 0\\\
-\dfrac{R_{x}}{R_{z}}\end{array}\right.\ .$ (39)

Analogously for the $\@vec{k}$ and $\@vec{k}_{p}$ vectors, the $BB^{\prime}$
line of Fig. 17 is defined by the intersection of the two planes of equations
$x=0$ and $y=1$ respectively and the coordinates of $\@vec{k}_{p}$ are
obtained by intersection with the line formation plane given in Eq. (20)

$\@vec{k}_{p}\left\\{\begin{array}[]{c}0\\\ 1\\\
-\dfrac{R_{y}}{R_{z}}\end{array}\right.\ .$ (40)

We now 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.
17, the surface of the projected $\Delta x\Delta y$ is $\@vec{i}_{p}\ \Delta
x\times\@vec{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 $\@vec{i}_{p}$ or $\@vec{k}_{p}$ with $\@vec{\ell}$,
which is the unit vector along the l.o.s.). The divergence is then obtained as

$\displaystyle\mathop{\mathrm{d}iv}\@vec{H}\ \Delta x\Delta y\Delta z$
$\displaystyle=$ $\displaystyle\left[\@vec{H}^{(m)}(z+\Delta
z)-\@vec{H}^{(m)}(z)\right]$ (41) $\displaystyle\cdot\left[\@vec{i}_{p}\
\Delta x\times\@vec{k}_{p}\ \Delta y\right]$
$\displaystyle+\left[\@vec{H}^{(m)}(x+\Delta x)-\@vec{H}^{(m)}(x)\right]$
$\displaystyle\cdot\left[\@vec{k}_{p}\ \Delta y\times\@vec{\ell}\ \Delta
z\right]$ $\displaystyle+\left[\@vec{H}^{(m)}(y+\Delta
y)-\@vec{H}^{(m)}(y)\right]$ $\displaystyle\cdot\left[\@vec{\ell}\ \Delta
z\times\@vec{i}_{p}\ \Delta x\right]$

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

###  A.2 Calculation of the curl

The curl is similarly obtained by applying the 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

Figure 18: 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. 15.
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.

We consider the calculation of $J_{Z}$ from the circulation of $\@vec{H}$
along the oriented frontier of the elementary surface perpendicular to $OZ$.
We have represented this elementary surface in Fig. 18, which is another view
of Fig. 15. The circulation of $\@vec{H}$ about the oriented surface is given
by the scalar products

$\displaystyle\left[\@vec{H}^{(m)}(y)-\@vec{H}^{(m)}(y+\Delta
y)\right]\cdot\@vec{I}\Delta\ell$ (42)
$\displaystyle+\left[\@vec{H}^{(m)}(x+\Delta
x)-\@vec{H}^{(m)}(x)\right]\cdot\@vec{K}\Delta y\ .$

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

$\@vec{I}\left\\{\begin{array}[]{l}\cos\theta\\\ 0\\\
-\sin\theta\end{array}\right.\ ,$ (43)

whereas $\@vec{K}$ is

$\@vec{K}\left\\{\begin{array}[]{l}0\\\ 1\\\ 0\end{array}\right.\ .$ (44)

As a consequence, the circulation of $\@vec{H}$ is

$\displaystyle\left[H_{x}^{(m)}(y)\cos\theta-
H_{z}^{(m)}(y)\sin\theta\right]\Delta\ell$
$\displaystyle+\left[H_{y}^{(m)}(x+\Delta x)\right]\Delta y$ (45)
$\displaystyle-\left[H_{x}^{(m)}(y+\Delta y)\cos\theta-H_{z}^{(m)}(y+\Delta
y)\sin\theta\right]\Delta\ell$
$\displaystyle-\left[H_{y}^{(m)}(x)\right]\Delta y\ .$

This circulation equals the flux of $\@vec{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. (28). This can be written as

$\displaystyle J_{Z}$ $\displaystyle=$
$\displaystyle\cos\theta\left[\frac{\Delta_{x}H_{y}^{(m)}}{\Delta
x}-\frac{\Delta_{y}H_{x}^{(m)}}{\Delta y}\right]$ (46)
$\displaystyle+\sin\theta\frac{\Delta_{y}H_{x}^{(m)}}{\Delta y}\ ,$

which is Eq. (3.3) of the paper in the case of the surface located at equator,
where the coordinates of the $\@vec{R}$ vector are given by Eq. (38) with
$\varphi=0$.

####  A.2.2 General case

Referring now to Fig. 17, the sides of the elementary volume are respectively
$\@vec{i}_{p}\Delta x$, $\@vec{k}_{p}\Delta y$ and $\@vec{\ell}\ \Delta z$.
The circulation of $\@vec{H}$ is given by its scalar product with each of
these elementary vectors, eventually reversed following the surface
orientation. This leads to

$\displaystyle J_{Z}\ \frac{\Delta x\Delta y}{R_{z}}$ $\displaystyle=$
$\displaystyle\@vec{J}\cdot\left[\@vec{i}_{p}\ \Delta x\times\@vec{k}_{p}\
\Delta y\right]$ $\displaystyle=$
$\displaystyle\left[\@vec{H}^{(m)}(y)-\@vec{H}^{(m)}(y+\Delta
y)\right]\cdot\@vec{i}_{p}\ \Delta x$
$\displaystyle+\left[\@vec{H}^{(m)}(x+\Delta
x)-\@vec{H}^{(m)}(x)\right]\cdot\@vec{k}_{p}\ \Delta y$

which results in the expression given in Eq. (3.3). 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

$\displaystyle J_{y}\ \Delta x\Delta z\ $ $\displaystyle=$
$\displaystyle\@vec{J}\cdot\left[\@vec{\ell}\ dz\times\@vec{i}_{p}\ \Delta
x\right]$ $\displaystyle=$ $\displaystyle\left[\@vec{H}^{(m)}(z+\Delta
z)-\@vec{H}^{(m)}(z)\right]\cdot\@vec{i}_{p}\ \Delta x$
$\displaystyle+\left[\@vec{H}^{(m)}(x)-\@vec{H}^{(m)}(x+\Delta
x)\right]\cdot\@vec{\ell}\ \Delta z$

and

$\displaystyle J_{x}\ \Delta y\Delta z\ $ $\displaystyle=$
$\displaystyle\@vec{J}\cdot\left[\@vec{k}_{p}\ \Delta y\times\@vec{\ell}\
\Delta z\right]$ $\displaystyle=$
$\displaystyle\left[\@vec{H}^{(m)}(z)-\@vec{H}^{(m)}(z+\Delta
z)\right]\cdot\@vec{k}_{p}\ \Delta y$
$\displaystyle+\left[\@vec{H}^{(m)}(y+\Delta
y)-\@vec{H}^{(m)}(y)\right]\cdot\@vec{\ell}\ \Delta z\ .$

This results in the expressions given in Eqs. (16-17). Again, for clarity, the
repeated and unchanged indices have been omitted in $\@vec{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.
3.3, in which all the coordinates of the current density vector are finally
derived in both reference frames.