© ESO, 2011

1. Introduction

Many inversion codes are available at present to derive magnetic field vectors from Zeeman Stokes profiles of solar lines. At the same time, these codes determine thermodynamical parameters of the solar atmosphere, which are coupled with magnetic fields in line formation. In general they employ least-squares inversion techniques, in which atmospheric parameters are modified by iterations until synthetic profiles match the observed ones (see for example the review by Bellot Rubio 2006). Some codes, SIR (Ruiz Cobo & del Toro Iniesta 1992) for example, use response functions to iteratively determine parameters at all altitudes across the atmosphere. NLTE extensions are possible (Socas-Navarro 2007). Some other codes such as UNNOFIT (Bommier et al. 2007) derive synthetic profiles with analytic expressions derived from the Milne-Eddington approximation and Unno-Rachkowsky solution (Harvey et al. 1972; Landolfi & Landi Degl’Innocenti 1982; Skumanich & Lites 1987). Multi-component atmospheres can also be investigated, with filling factors or micro-structures (MISMA code Sanchez Almeida 1997). Some other inversion codes are based on a direct comparison of observed profiles with a number of synthetic ones. They can be strongly accelerated by a principal component analysis (Rees et al. 2000; Socas-Navarro et al. 2001).

Central moments of unpolarized line profiles have been used in the past to determine velocity fields across the solar atmosphere (see for example Marmolino & Severino 1981). In the case of Zeeman sensitive lines, first order moments can be used to determine the center of gravity of the left- and right-hand circularly polarized profiles and estimate the longitudinal magnetic field (Semel 1967; Uitenbroek 2003). In this paper we expand on the use of moments, and propose the use of several quantities derived from central moments of the Stokes profiles I ± Q,U in addition to those of I ± V. These quantities turn out to be nearly linear funtions of the magnetic field components along and perpendicular to the line of sight, and are mostly insensitive to other atmospheric parameters. They can be used to determine vector magnetic fields with short polynomial expansions that can be properly calibrated with forward modeling in the context of a small number of one-dimensional model atmospheres.

In Sect. 2 we define the moments that we use in Sect. 3 to find the basic quantities to be used in the inversions. In Sect. 4 we present the radiative transfer code and the four model atmospheres used to compute synthetic profiles and inversion coefficients. In Sects. 5 and 6, polynomial expansions are specified and tested in the case of one model atmosphere. Noise effects are estimated in Sect. 7. To satisfy simultaneously a great number of solar conditions, an interpolation method is proposed and tested with different models and center-to-limb positions (Sect. 8). It is compared to UNNOFIT inversion of THEMIS/MTR data (Sect. 9). Because of the low sensitivity to thermodynamical parameters, the capabilities of a simpler method using one average model only are also investigated (Sect. 10). We end with conclusions in Sect. 11.

2. Moments of line profiles

We denote with s the Zeeman splitting corresponding to the longitudinal field B_∥, for an effective Lande factor $\overline{\mathit{g}}$ and a line-center wavelength λ₀: $\mathit{s}\mathrm{=}\mathit{k}{\mathit{B}}_{\mathrm{\parallel }}$(1)with $\mathit{k}\mathrm{=}\mathrm{4.67}\mathrm{\times }{\mathrm{10}}^{-14}\overline{\mathit{g}}\hspace{0.17em}{\mathit{\lambda }}_{\mathrm{0}}^{\mathrm{2}}\mathit{,}$(2)where wavelengths are expressed in nm and B_∥ in Gauss. The depth of a line profile P(λ) can be written $\mathit{R}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{-}\mathit{P}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathit{/}{\mathit{I}}_{\mathrm{c}}\mathit{,}$(3)where I_c is the continuum intensity. We note with ${\mathit{M}}_{\mathrm{1}}\mathrm{(}\mathit{P,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{1}}\mathrm{)}\mathrm{=}\frac{{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{R}}^{{\mathit{\beta }}_{\mathrm{1}}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathit{\lambda }\hspace{0.17em}\mathrm{d}\mathit{\lambda }}{\mathit{k}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{R}}^{{\mathit{\beta }}_{\mathrm{1}}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\lambda }}$(4)the first order moment of R, expressed in Gauss. The exponent β₁ is used to assign different weights to line wings and line center and to reduce noise effects (Sect. 8).

To characterize the widths of line profiles, we define a second moment of R similar to the classical second-order central moment: ${\mathit{M}}_{\mathrm{2}}\mathrm{(}\mathit{P,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{2}}\mathrm{)}\mathrm{=}\frac{{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{R}}^{{\mathit{\beta }}_{\mathrm{2}}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\hspace{0.17em}\mathrm{|}\mathit{\lambda }\mathrm{-}{\mathit{\lambda }}_{\mathrm{2}}\mathrm{|}\hspace{0.17em}\mathrm{d}\mathit{\lambda }}{\mathit{k}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{R}}^{{\mathit{\beta }}_{\mathrm{2}}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\lambda }}\mathit{,}$(5)where λ₂ is the wavelength corresponding to the first moment of R(λ) to the power β₂: ${\mathit{\lambda }}_{\mathrm{2}}\mathrm{=}\mathit{k}{\mathit{M}}_{\mathrm{1}}\mathrm{(}\mathit{P,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{2}}\mathrm{)}\mathit{.}$(6)We use the modulus of λ − λ₂ instead of (λ − λ₂)² to obtain the exact shifts of blue and red components in the case of strong transverse magnetic fields (see next section).

3. Quasi-linear functions A₁ and A₂ of magnetic field components, from weak to strong fields

For a weak line formed in a Milne-Eddington atmosphere without velocity and magnetic field, the line depth can be written as $\mathit{R}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}\mathit{c}{\mathrm{e}}^{\mathrm{-}\mathrm{(}\mathit{\lambda }\mathrm{-}{\mathit{\lambda }}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}\mathit{/}{\mathit{d}}^{\mathrm{2}}}$(7)by neglecting damping effects. The coefficient c and the mean Doppler-width d depend on the line and model atmosphere.

Let us now assume a constant vector magnetic field throughout the atmosphere, and denote with B_∥ and B_⊥ the longitudinal and the transverse components of that field. Below we show that one can express the moments of line depths corresponding to polarized profiles I ± S, (S = Q,U,V) as functions of B_∥ and B _⊥ . In this we restrict our calculations to the case of a Zeeman triplet.

3.1. B_∥ and the function A₁

Let us assume B_⊥ = 0. As before, we denote with s the wavelength shift corresponding to the longitudinal field B_∥ and λ₀ the wavelength at line center. The line depths of I − V and I + V are two Gaussian functions shifted by  ± s from λ₀: $\begin{array}{ccc}{\mathit{R}}_{\mathit{I}\mathrm{+}\mathit{V}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}\mathit{c}{\mathrm{e}}^{\mathrm{-}\mathrm{(}\mathit{\lambda }\mathrm{-}{\mathit{\lambda }}_{\mathrm{0}}\mathrm{-}\mathit{s}{\mathrm{)}}^{\mathrm{2}}\mathit{/}{\mathit{d}}^{\mathrm{2}}}\mathit{,}& & \\ {\mathit{R}}_{\mathit{I}\mathrm{-}\mathit{V}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}\mathit{c}{\mathrm{e}}^{\mathrm{-}\mathrm{(}\mathit{\lambda }\mathrm{-}{\mathit{\lambda }}_{\mathrm{0}}\mathrm{+}\mathit{s}{\mathrm{)}}^{\mathrm{2}}\mathit{/}{\mathit{d}}^{\mathrm{2}}}\mathit{.}& & \end{array}$The first-order moments M₁(I + V,   β₁) and M₁(I − V,   β₁) provide the centers of gravity of line depths R_I + V(λ) and R_I − V(λ) so that the quantity ${\mathit{A}}_{\mathrm{1}}\mathrm{\equiv }\left[{\mathit{M}}_{\mathrm{1}}\mathrm{(}\mathit{I}\mathrm{+}\mathit{V,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{1}}\mathrm{)}\mathrm{-}{\mathit{M}}_{\mathrm{1}}\mathrm{(}\mathit{I}\mathrm{-}\mathit{V,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{1}}\mathrm{)}\right]\mathit{/}\mathrm{2}$(10)is equal to B_∥. This is the well-known center of gravity method for determining the line-of-sight magnetic field (Rees & Semel 1979). The result in Eq. (10) is independent of c, d and β₁.

3.2. B_⊥ and the function A₂

Let us assume now that B_∥ = 0 and that the azimuth of B _⊥  is parallel to the polarization direction of Q. The line depths of I − Q and I + Q can be written $\begin{array}{ccc}{\mathit{R}}_{\mathit{I}\mathrm{-}\mathit{Q}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}& \mathrm{=}& \mathit{c}{\mathrm{e}}^{\mathrm{-}\mathrm{(}\mathit{\lambda }\mathrm{-}{\mathit{\lambda }}_{\mathrm{0}}\mathrm{-}\mathit{w}{\mathrm{)}}^{\mathrm{2}}\mathit{/}{\mathit{d}}^{\mathrm{2}}}\mathrm{+}\mathit{c}{\mathrm{e}}^{\mathrm{-}\mathrm{(}\mathit{\lambda }\mathrm{-}{\mathit{\lambda }}_{\mathrm{0}}\mathrm{+}\mathit{w}{\mathrm{)}}^{\mathrm{2}}\mathit{/}{\mathit{d}}^{\mathrm{2}}}\\ {\mathit{R}}_{\mathit{I}\mathrm{+}\mathit{Q}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}& \mathrm{=}& \mathrm{2}\mathit{c}{\mathrm{e}}^{\mathrm{-}\mathrm{(}\mathit{\lambda }\mathrm{-}{\mathit{\lambda }}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}\mathit{/}{\mathit{d}}^{\mathrm{2}}}\mathit{,}\end{array}$with $\mathit{w}\mathrm{=}\mathit{k}{\mathit{B}}_{\mathrm{\perp }}\mathit{.}$(13)

3.2.1. Weak transverse components

If the transverse component is weak (w ≪ d), we can write ${\mathit{M}}_{\mathrm{2}}\mathrm{(}\mathit{I}\mathrm{-}\mathit{Q,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{2}}\mathrm{)}\mathrm{\simeq }\frac{\mathrm{(}{\mathit{d}}^{\mathrm{2}}\mathrm{+}{\mathit{w}}^{\mathrm{2}}\mathrm{)}}{\mathit{kd}\sqrt{\mathit{\pi }{\mathit{\beta }}_{\mathrm{2}}}}$(14)and ${\mathit{M}}_{\mathrm{2}}\mathrm{(}\mathit{I}\mathrm{+}\mathit{Q,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{2}}\mathrm{)}\mathrm{\simeq }\frac{\mathit{d}}{\mathit{k}\sqrt{\mathit{\pi }{\mathit{\beta }}_{\mathrm{2}}}}$(15)so that ${\mathit{M}}_{\mathrm{2}}\mathrm{(}\mathit{I}\mathrm{-}\mathit{Q,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{2}}\mathrm{)}\mathrm{-}{\mathit{M}}_{\mathrm{2}}\mathrm{(}\mathit{I}\mathrm{+}\mathit{Q,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{2}}\mathrm{)}\mathrm{\simeq }\frac{{\mathit{w}}^{\mathrm{2}}}{\mathit{kd}\sqrt{\mathit{\pi }{\mathit{\beta }}_{\mathrm{2}}}}$(16)and $\mathit{w}\mathrm{\simeq }{\mathrm{(}\sqrt{\mathit{\pi }{\mathit{\beta }}_{\mathrm{2}}}\mathit{kd}\hspace{0.17em}{\left[{\mathit{M}}_{\mathrm{2}}\mathrm{(}\mathit{I}\mathrm{-}\mathit{Q,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{2}}\mathrm{)}\mathrm{-}{\mathit{M}}_{\mathrm{2}}\mathrm{(}\mathit{I}\mathrm{+}\mathit{Q,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{2}}\mathrm{)}\right]}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(17)From Eq. (15) we see that $\mathit{d}\mathit{/}\mathit{k}\mathrm{=}\sqrt{\mathit{\pi }{\mathit{\beta }}_{\mathrm{2}}}{\mathit{M}}_{\mathrm{2}}\mathrm{(}\mathit{I}\mathrm{+}\mathit{Q,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{2}}\mathrm{)}\mathit{.}$(18)To minimize cross-talk effects between the determination of the longitudinal and transverse field components, we replace M₂(I + Q,   β₂) by the average value W of the widths of I + V and I − V (see Eq. (5)): $\mathit{W}\mathrm{=}\left[{\mathit{M}}_{\mathrm{2}}\mathrm{(}\mathit{I}\mathrm{+}\mathit{V,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{2}}\mathrm{)}\mathrm{+}{\mathit{M}}_{\mathrm{2}}\mathrm{(}\mathit{I}\mathrm{-}\mathit{V,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{2}}\mathrm{)}\right]\mathit{/}\mathrm{2}\mathit{,}$(19)so that (from Eqs. (13), (17) and (18)) ${\mathit{B}}_{\mathrm{\perp }}\mathrm{=}\sqrt{\mathit{\pi }{\mathit{\beta }}_{\mathrm{2}}}{\left(\mathit{W}\left[{\mathit{M}}_{\mathrm{2}}\mathrm{(}\mathit{I}\mathrm{-}\mathit{Q,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{2}}\mathrm{)}\mathrm{-}{\mathit{M}}_{\mathrm{2}}\mathrm{(}\mathit{I}\mathrm{+}\mathit{Q,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{2}}\mathrm{)}\right]\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(20)Combining Eqs. (19) and (20) we see that, for weak transverse magnetic fields, $\begin{array}{ccc}{\mathit{A}}_{\mathrm{2}}& \mathrm{\equiv }& \\ & & {\mathrm{\times }\left({\mathit{M}}_{\mathrm{2}}\mathrm{(}\mathit{I}\mathrm{-}\mathit{Q,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{2}}\mathrm{)}\mathrm{-}{\mathit{M}}_{\mathrm{2}}\mathrm{(}\mathit{I}\mathrm{+}\mathit{Q,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{2}}\mathrm{)}\right)]}^{\mathrm{1}\mathit{/}\mathrm{2}}\end{array}$(21)is a linear function of the transverse component B _⊥ , with a proportionality coefficient of 1/$\sqrt{\mathit{\pi }{\mathit{\beta }}_{\mathrm{2}}}$.

3.2.2. Strong transverse components

In the case of strong transverse fields (w ≫ d) ${\mathit{M}}_{\mathrm{2}}\mathrm{(}\mathit{I}\mathrm{-}\mathit{Q,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{2}}\mathrm{)}\mathrm{\simeq }\mathit{w}\mathit{/}\mathit{k,}$(22)so that the component B _⊥  is provided directly by the moment M₂(I − Q,   β₂). But at the same time M₂(I + Q,   β₂) becomes much smaller than M₂(I − Q,   β₂). We will see later that in realisitic cases (Sect. 5, Fig. 2) the width W of profiles I + V and I − V increases for large B _⊥ , in such a way that A₂ remains close to a linear function of B _⊥  for strong fields as well as for weak fields (see Eq. (21)). It thus becomes possible to use A₂ in the full range of magnetic field strengths. In the following sections, we will express A₂ in Gauss, regardless of the (dimensionless) proportionality coefficient, which is close to (πβ₂) ^− 1/2 for weak fields.

3.3. Magneto-optical effects

We call Φ the angle characterizing magneto-optical effects: $\mathrm{\Phi }\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\arctan \mathrm{\left(}\mathit{U}\mathrm{\right(}{\mathit{\lambda }}_{\mathrm{0}}\mathrm{)}\mathit{/}\mathit{Q}\mathrm{(}{\mathit{\lambda }}_{\mathrm{0}}\mathrm{\left)}\mathrm{\right)}\mathit{,}$(23)with ${\mathit{\lambda }}_{\mathrm{0}}\mathrm{=}\mathit{k}{\mathit{M}}_{\mathrm{1}}\mathrm{\left(}\mathit{P,}\mathrm{1}\mathrm{\right)}$(24)$\mathrm{\Phi }\mathit{U}\mathrm{\left(}{\mathit{\lambda }}_{\mathrm{0}}\mathrm{\right)}\mathrm{\le }\mathrm{0}$(25)and $\mathrm{-}{\mathrm{90}}^{\mathrm{0}}\mathit{<}\mathrm{\Phi }\mathit{<}{\mathrm{90}}^{\mathrm{0}}\mathit{,}$(26)where U(λ₀) and Q(λ₀) are the values of the Stokes parameters U and Q at line center of the intensity profile.

We must note that while magnetic components are related to moments integrated over full profiles, the magneto-optical effects are determined by Stokes parameters observed in one single wavelength. This leads to a higher sensitivity to noise, especially for weak transverse magnetic components, corresponding to low Q and U values. Although other more sophisticated algorithms might be investigated in this respect, we restrict the present calculations to this most simple and least time-consuming method.

4. The RH radiative transfer code

To calibrate the inversion method described here and to test its efficiency, we obtained forward solutions of the polarized radiative transfer equations with the RH radiative transfer code (Uitenbroek 2001, 2003) which is based on the Rybicki & Hummer (1991, 1992) multi-level accelerated lambda iteration scheme. The coupled equations of statistical equilibrium and radiative transfer were solved for a 23 level, 33 lines atomic model of Fe i, including the 630.25 and 630.15 nm lines. These non-LTE iterations were performed in the polarization free approximation to account for the effect of the line profile splitting on the radiative rates (Bruls & Trujillo Bueno 1996). After convergence the four Stokes parameters were calculated with the quasi parabolic DELO method developed by Trujillo Bueno (2003). In the calculation of the absorption and emission coefficients for the four Stokes parameters, the code fully accounts for the actual Zeeman splitting patterns of the atomic levels.

Four different one-dimentional hydrostatic solar atmospheric models were employed in the calculations: FALA, FALC, FALF (Fontenla et al. 1993) and MALTM (Maltby et al. 1986), representing the stratification of the quiet cell interior, averaged quiet sun, solar network and a sunspot umbra, respectively.

5. Functions A₁, A₂, W and Φ for the Fe i 630.2 nm line with the FALC model at disk center

Fig. 1 Functions A1, A2 and Φ (Eqs. (10), (21), (23)), measured in the case of the Fe i 630.2 nm line in the FALC model at disk center. The line styles depend on the values of B ⊥  in panels a) and c), and the values of B∥ in panel b): full lines until 1000 G, dashed lines until 2000 G, dashed-dotted until 3000 G.

We computed values for A₁, W , A₂ and Φ for the Fe i 630.2 nm line in the FALC model at disk center for magnetic fields in the range of  −3000 < B_∥ < 3000 G and 0 < B _⊥  < 3000 G, generally suitable for quiet Sun and solar active regions, with steps of 200 G for both components. The step of wavelength sampling was 2 pm after interpolation of the synthetic profiles. The exponent values in the moment computations (Eqs. (4) and (5)) are $\begin{array}{ccc}{\mathit{\beta }}_{\mathrm{1}}& \mathrm{=}& \mathrm{1}\\ {\mathit{\beta }}_{\mathrm{2}}& \mathrm{=}& \mathrm{2}\mathit{,}\end{array}$according to optimizations that were chosen to minimize noise effects, as detailed in Sect. 8. Figure 1, panels (a) and (b), shows that A₁ and A₂ are nearly linear functions of B_∥ and B _⊥  over all ranges of considered field strengths. Magneto-optical effects are plotted in panel (c) of Fig. 1. In Fig. 2 we note the increase of W, which accounts for the almost linear asymptotic behavior of A₂ for strong transverse fields, and motivated the choice of substituting W for w in the definition of A₂. The increase in W is mainly caused by the increasing depth of the central core, because of the cross-talk of I ± V with the transverse component. Note that the obtained linearity of function A₁ with B_∥ provides clear evidence of the accuracy of the center-of-gravity determination of the line-of-sight component of the magnetic field for a large selection of longitudinal and transverse field component strentghs.

6. Polynomials P_1, FALC, P_2, FALC, and inversion of Fe i 630.2 nm profiles from FALC at disk center

For all synthetic profiles we have computed A₁ and A₂ after azimuth rotation over Φ, as indicated in Sect. 3.3. We first analyze the functions A₁, A₂ and Φ that were computed with the FALC model as functions of B_∥ and B_⊥ in the full range of useful magnetic field components.

A₁ and Φ always have the same sign as B_∥ and A₂ does not depend on the sign of B_∥. To reduce the number of coefficients, the following equations are supposed to be applied to the absolute value of A₁. If A₁ < 0, the final values of B_∥ and Φ must be replaced by  −B_∥ and  −Φ.

Because A₁ is a monotonic, almost linear function of B_∥, we can expand B_∥ as a low-degree polynomial in A₁, and obtain the coefficients with a least-squares fit. The coefficients depend on B_⊥ (see Fig. 1), so that

${\mathit{B}}_{\mathrm{\parallel }}\mathrm{=}\sum _{\mathit{m}}{\mathit{p}}_{\mathrm{1}\mathit{,m}}\mathrm{\left(}{\mathit{B}}_{\mathrm{\perp }}\mathrm{\right)}{\mathit{A}}_{\mathrm{1}}^{\mathit{m}\mathrm{-}\mathrm{1}}\mathit{.}$(29)We use polynomials with m ≤ 7. Because the coefficients p_1,m(B _⊥ ) are known for all values of (B _⊥ ), they can also be expanded as a polynomial in B _⊥ , ${\mathit{p}}_{\mathrm{1}\mathit{,m}}\mathrm{\left(}{\mathit{B}}_{\mathrm{\perp }}\mathrm{\right)}\mathrm{=}\sum _{\mathit{n}}{\mathit{C}}_{\mathrm{1}\mathit{,m,n}}{\mathit{B}}_{\mathrm{\perp }}^{\mathit{n}\mathrm{-}\mathrm{1}}\mathit{,}$(30)again with n ≤ 7 and the coefficients determined by a least-squares fit from the results of the forward modeling in the FALC model. Combining both expansions, we find ${\mathit{B}}_{\mathrm{\parallel }}\mathrm{=}\sum _{\mathit{m}}\sum _{\mathit{n}}{\mathit{C}}_{\mathrm{1}\mathit{,m,n}}{\mathit{A}}_{\mathrm{1}}^{\mathit{m}\mathrm{-}\mathrm{1}}{\mathit{B}}_{\mathrm{\perp }}^{\mathit{n}\mathrm{-}\mathrm{1}}\mathit{.}$(31)In the same way, because A₂ is a monotonic function of B _⊥ , we can determine coefficients p_2,m(B_∥) so that ${\mathit{B}}_{\mathrm{\perp }}\mathrm{=}\sum _{\mathit{m}}{\mathit{p}}_{\mathrm{2}\mathit{,m}}\mathrm{\left(}{\mathit{B}}_{\mathrm{\parallel }}\mathrm{\right)}{\mathit{A}}_{\mathrm{2}}^{\mathit{m}\mathrm{-}\mathrm{1}}\mathit{,}$(32)and we can expand p_2,m(B_∥) as a function of B_∥${\mathit{p}}_{\mathrm{2}\mathit{,m}}\mathrm{\left(}{\mathit{B}}_{\mathrm{\parallel }}\mathrm{\right)}\mathrm{=}\sum _{\mathit{n}}{\mathit{C}}_{\mathrm{2}\mathit{,m,n}}{\mathit{B}}_{\mathrm{\parallel }}^{\mathit{n}\mathrm{-}\mathrm{1}}\mathit{,}$(33)so that ${\mathit{B}}_{\mathrm{\perp }}\mathrm{=}\sum _{\mathit{m}}\sum _{\mathit{n}}{\mathit{C}}_{\mathrm{2}\mathit{,m,n}}{\mathit{A}}_{\mathrm{2}}^{\mathit{m}\mathrm{-}\mathrm{1}}{\mathit{B}}_{\mathrm{\parallel }}^{\mathit{n}\mathrm{-}\mathrm{1}}\mathit{.}$(34)In the following, we denote with P_1,FALC and P_2,FALC the polynomials of coefficients C_1,m,n and C_2,m,n: $\begin{array}{ccc}{\mathit{B}}_{\mathrm{\parallel }}& \mathrm{=}& \\ {\mathit{B}}_{\mathrm{\perp }}& \mathrm{=}& \end{array}$The longitudinal and transverse field components B_∥ and B _⊥  can be derived from the measured quantities A₁ and A₂ by the following fast two-step iteration scheme. In the first step, B_∥ is initialized with A₁ and substituted into Eq. (36). Next the values of the two field components are derived by successive backsubstitutions: $\begin{array}{ccc}\mathit{X}& \mathrm{=}& \\ \mathit{Y}& \mathrm{=}& {\mathit{P}}_{\mathrm{1}\mathit{,}\mathrm{FALC}}\mathrm{\left(}{\mathit{A}}_{\mathrm{1}}\mathit{,X}\mathrm{\right)}\\ {\mathit{B}}_{\mathrm{\perp }}& \mathrm{=}& {\mathit{P}}_{\mathrm{2}\mathit{,}\mathrm{FALC}}\mathrm{\left(}{\mathit{A}}_{\mathrm{2}}\mathit{,Y}\mathrm{\right)}\\ {\mathit{B}}_{\mathrm{\parallel }}& \mathrm{=}& \end{array}$Magneto-optical effects can be deduced from the synthetic function Φ that we can also expand as a polynomial in B_∥ and B _⊥ : $\begin{array}{ccc}\mathrm{\Phi }& \mathrm{=}& \\ \mathrm{\Phi }& \mathrm{=}& \end{array}$Because Φ is not defined for B_⊥ = 0, extrapolated values of P_3,FALC(B_∥,0) are used before computing the coefficients C_3,m,n.

Fig. 2 Width W (Eq. (19)) in the case of the Fe i 630.2 nm line in the FALC model at disk center. Line styles are the same as in Fig. 1, panel b).

Fig. 3 Comparison of longitudinal (panel a)) field strength, transverse field strength (panel b)), and azimuth (panel c)) recovered with the inversion with input values.

To test the inversion code based on moments, we deduce vector magnetic fields from the values of A₁, A₂ calculated with the forward modeling in model FALC and the procedure outlined above in Eqs. (37) to (42). Figure 3 shows the residual errors of B_∥ and B_⊥ resulting from the inversion process. The standard deviations σ_∥ and σ_⊥ are only 1 G and 4 G, respectively, over the full range of vector magnetic fields smaller than 3000 G, when least-squares polynomials are extended to 6th order. For results of the azimuth φ, we restrict the calculation to transverse magnetic fields larger or equal to 200 G, because the azimuth is not defined for zero fields. The standard deviation σ_φ is smaller than 1°.

7. Effects of noise

To test the effects of noise on the results, we added noise to the synthetic profiles at a level of 1.5 × 10^-3 relative to the continuum intensity. This is the level assumed for typical THEMIS/MTR data (Bommier et al. 2007). This noise is simulated by a random number taken between  −3 × 10^-3I_c and 3 × 10^-3I_c. Four different stochastic values are added to quantities I,Q,U and V.

Inversion errors are sensitive not only to the noise level, but also to the wavelength sampling of data. To mimic real observations, we used the wavelength sampling 2 pm, which is also a typical THEMIS/MTR value.

To minimize noise effects in the far line wings for moment calculations, we defined a new continuum value through the intensity at the highest and lowest wavelength in each profile. All data points with higher intensity are set to this value. In a similar way, any negative value of ${\mathit{A}}_{\mathrm{2}}^{\mathrm{2}}$ (Eq. (21)), which may appear because of noise in the case of low transverse magnetic field, is set to zero. Because of the limited number of points in the profiles, and because the results depend on the sets of stochastic numbers used to modify the synthetic data, the standard deviations are obtained by averaging over 100 successive inversions. The results which we present in Fig. 4 have standard deviations of 21 G and 29 G for both magnetic components, and 1.6° for azimuths. As expected, errors in the derived transverse component increase for very weak magnetic fields. The results are similar to the results derived from the UNNOFIT inversion of theoretical profiles (Bommier et al. 2007, their Fig. 4).

Fig. 4 Inversion of synthetic profiles for the FALC model atmosphere at disk center with noise a level of 1.5 × 10-3Ic and a sampling of 2 pm. Line styles are the same as in Fig. 1.

8. Data from several models: interpolation using A_D functions; optimized β₁ and β₂ exponents

8.1. Definition of the function A_D

The widths of the I + V and I − V profiles deduced with Eq. (19) and included in the A₂ functions (Eq. (21)) vary with the different employed model atmospheres. This leads to slight variations of the coefficients in Eqs. (35) and (36), depending on which atmospheric model is used for the polynomial fits, and this leads to variations in the values of the field that are recovered by the inversions. This model dependency needs to be accounted for. So far, we used the choice exponent β₂ = 2 in Eq. (19), which was optimized to decrease the sensitivity of the inversion to noise. Setting the exponent β₂ = 1 appears to yield a quantity that is more sensitive to the model.

To deal with model dependencies more adequately, we define a new quantity A_D similar to a Doppler width: ${\mathit{A}}_{\mathrm{D}}\mathrm{=}\mathrm{\left[}{\mathit{M}}_{\mathrm{2}}\mathrm{\right(}\mathit{I}\mathrm{+}\mathit{V,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{D}}\mathrm{)}\mathrm{+}{\mathit{M}}_{\mathrm{2}}\mathrm{(}\mathit{I}\mathrm{-}\mathit{V,}\hspace{0.17em}{\mathit{\beta }}_{\mathrm{D}}\mathrm{\left)}\mathrm{\right]}\mathit{/}\mathrm{2}\mathit{,}$(43)where we set ${\mathit{\beta }}_{\mathrm{D}}\mathrm{=}\mathrm{1.}$(44)For each model M, A_D,M is a monotonic function of A₂. It can be deduced from A₂ and B_∥ by an equation similar to Eq. (34): $\begin{array}{ccc}{\mathit{A}}_{\mathrm{D}\mathit{,M}}& \mathrm{=}& \\ {\mathit{A}}_{\mathrm{D}\mathit{,M}}& \mathrm{=}& \end{array}$Figure 5 shows the A_D-functions corresponding to the two models FALC and MALTM. They look like the curves of Fig. 2 plotted for FALC with respect to B _⊥ , except that A₂ varies roughly from 0 to 2000 G, while B _⊥  varies from 0 to 3000 G. Accordingly, the relevant model of any solar point can be specified by the comparison between the observed value of A_D and the value derived from Eq. (46), where B_∥ can be replaced by A₁ as a first approximation.

Fig. 5 Widths AD of profiles I ± V computed with two selected models for inversion addressing several models.

8.2. Inversion method

In our experiments we found that our inversion procedure yields sufficiently accurate results when we used a linear combination of polynomials P′ of only two selected models: FALC and MALTM. We therefore propose an inversion procedure using linear combinations of the polynomials P_FALC and P_MALTM determined by fitting the functions A₁, A₂, Φ and A_D for these models so that $\begin{array}{ccc}{\mathit{P}}_{\mathrm{1}\mathit{,\delta }}& \mathrm{=}& \\ {\mathit{P}}_{\mathrm{2}\mathit{,\delta }}& \mathrm{=}& \\ {\mathit{P}}_{\mathrm{3}\mathit{,\delta }}& \mathrm{=}& \end{array}$with the fractional coefficient given by $\mathit{\delta }\mathrm{=}\frac{{\mathit{A}}_{\mathrm{D}}\mathrm{-}{\mathit{P}}_{\mathrm{D}\mathit{,}\mathrm{FALC}}\mathrm{\left(}{\mathit{A}}_{\mathrm{2}}\mathit{,}{\mathit{A}}_{\mathrm{1}}\mathrm{\right)}}{{\mathit{P}}_{\mathrm{D}\mathit{,}\mathrm{MALTM}}\mathrm{\left(}{\mathit{A}}_{\mathrm{2}}\mathit{,}{\mathit{A}}_{\mathrm{1}}\mathrm{\right)}\mathrm{-}{\mathit{P}}_{\mathrm{D}\mathit{,}\mathrm{FALC}}\mathrm{\left(}{\mathit{A}}_{\mathrm{2}}\mathit{,}{\mathit{A}}_{\mathrm{1}}\mathrm{\right)}}$(50)according to Eq. (46), where A₁ approximates B_∥.

We show below that these linear interpolations are accurate enough because the inversion polynomials are relatively unsensitive to model parameters (Sect. 8.6).

With these modifications the inversion scheme now looks as follows:

• 1.

  determine λ₀ (Eq. (4)) and the angle Φ (Eqs. (23), (25) and (26));

• 2.

  determine new Stokes Q parameter after rotation over Φ;

• 3.

  calculate A₁, A₂, A_D and δ (Eqs. (10), (21), (46) and (50));

• 4.

  evaluate polynomials P_1,δ and P_2,δ with known coefficients and δ (Eqs. (47) and (48));

• 5.

  determine B_∥ and B _⊥  (Eqs. (37) to (40)), with P_δ instead of P_FALC;

• 6.

  determine azimuth φ taking into account magneto-optical effects: $\mathit{\phi }\mathrm{=}\mathrm{\Phi }\mathrm{-}{\mathit{P}}_{\mathrm{3}\mathit{,\delta }}\mathrm{\left(}{\mathit{B}}_{\mathrm{\parallel }}\mathit{,}{\mathit{B}}_{\mathrm{\perp }}\mathrm{\right)}\mathit{,}$(51)where P_3,δ is deduced fromp Eq. (49).

Fig. 6 Optimization of exponents β1 and β2 for a low sensitivity to data noise (see Sect. 8.3).

8.3. Optimized exponents β₁ and β₂

The accuracy of inversion from data not only depends on data noise, but also on the values of β₁ and β₂. If the exponents are too small, the largest weight is assigned to the line wings, which are noise-generating, while not very sensitive to magnetic fields. On the other hand, if they are too large, only a few points of the line core are effectively involved in the moments calculations, and the effects of noise increase.

To optimize the exponents β₁ and β₂, we inverted the synthetic profiles of models FALC and MALTM, disturbed by a noise level 1.5 × 10^-3I_c, using different values of β₁ and β₂. We characterize the errors by standard deviations σ_∥ and σ _⊥  corresponding to the mean results from inversions in both models, for field vectors smaller than 3000 G. Figure 6 presents the standard deviations for the parallel (panel a) and perpendicular (panel b) components, respectively, as functions of the parameters β₁ and β₂ for 100 successive realizations in each case. Noise effects are roughly minimized by the choice of integer values β₁ = 1 and β₂ = 2 used in previous sections. In Fig. 6 the results corresponding to β₁ are plotted for β₂ = 2 and the results corresponding to β₂ are plotted for β₁ = 1.

8.4. Results

Table 1 presents the standard deviations between input and recovered values of magnetic field components and azimuths in the cases of synthetic profiles from models FALC, MALTM, FALA and FALF, with and without noise. In all cases the wavelegth step is 2 pm. Averages are extended to all vector magnetic fields less than 3000 G, except for azimuth starting from transverse magnetic fields larger than 200 G.

Mean values  ⟨ δ ⟩  characterizing the widths of line profiles are also presented. As expected, they are close to 0 and 1 for FALC and MALTM respectively. For a given model, they increase generally with noise. The reason is probably that according to the definition of the continuum (Sect. 7), the noise decreases the mean far wing intensities and increases A_D. This increase of line widths through noise may account for the result obtained for the deviation σ_∥ of FALC data with δ = 0.05, which is better than in the case δ = 0 (Fig. 4a).

The standard deviation increases significantly with noise, showing that the errors caused by the inversion process itself can be generally neglected for real observed data.

It must be noted that the wavelength sampling used for the computation of inversion coefficients and for the inverted profiles is exactly the same as in Sect. 5. If additional interpolations are used between both sets of profiles, the σ-values without noise increase by a few Gauss, all the more because δ is also slightly affected through the A_D value. However, the results for noisy profiles are almost unchanged.

8.5. Center-to-limb effects

In Table 1 additional σ-values are given for the inversion of synthetic profiles in the FALC model corresponding to different center-to-limb positions. They are higher than values computed at disk center, but still acceptable. For any μ-value across the solar disk it is always possible to calibrate the exact coefficients of P₁ and P₂ from synthetic profiles at corresponding viewing angles and to obtain very low σ-values.

8.6. Sensitivity to model parameters

To test the sensitivity to the input models used to calibrate the expansion coefficients for the polynomials P, it is interesting to degrade the inversion code and to replace the interpolation between P_FALC and P_MALTM by the mean value (P_FALC + P_MALTM)/2. Corresponding calculations have been performed and results similar to results of Table 1 are presented in Table 2.

For data without noise (lefthand parts of the tables), the higher accuracy of inversions with interpolations is obvious. The only exception is the last line of the tables which shows that the inversion of FALC data from μ = 0.5 is more accurate with the average of P_FALC and P_MALTM from the disk center than with P_FALC from the disk center.

For noisy data (righthand parts) the average standard deviations σ_∥ and σ _⊥  for disk center increase roughly by factors 1.7 and 1.4 only from Tables 1 to 2. This means that for this noise level and this set of models, interpolations that are more sophisticated than the linear one used in Sect. 8.2 do not seem to be required.

9. Comparison with UNNOFIT inversion of THEMIS/MTR data⁵

To test our inversion code with real data, we used THEMIS/MTR observations obtained by Bommier on 2005 September 13, which were already presented in a paper dealing with azimuth ambiguity solution (Bommier et al. 2011). They concern the 630.2 Fe i line and an active region located at L =  −3.8° and b =  −11.7°. Taking into account the slow variation of inversion coefficients around the disk center, as indicated in Tables 1 and 2, we can use the same disk-center polynomials as before to invert the observed profiles.

Before computing the moments for the THEMIS data, we needed to suppress the telluric blend in the red wing of the line, because we deal with ground-based data. We replaced nine points of the observed profile by a least-squares cubic interpolation. The corresponding wavelength interval is roughly 20 pm.

The exponents in the calculations of A₁, A₂ and A_D are β₁ = 1, β₂ = 2 and β_D = 1 respectively. The investigated field of view is 90″  ×  60″.

The UNNOFIT inversion code is based on the Levenberg-Marquardt algorithm applied to the theoretical profile derived from the Unno-Rachkovsky solution of the polarized radiative transfer equation through a Milne-Eddington atmosphere (Landolfi et al. 1984) The best way to use UNNOFIT inversions is to assume inhomogeneous fields with filling factors f different from 1, and to determine the products f ∗ B (Bommier et al. 2007). But our ICM inversion does not include this model so far. Because UNNOFIT codes do exist also with f = 1, we made our test by comparing results obtained in that case.

Fig. 7 THEMIS/MTR data: comparison between ICM and UNNOFIT results assuming a filling factor f = 1.

The ICM computing time is very short, roughly 3 × 10^-4 s per solar point with a Dell R710 computer (2.67 GHz). The results are plotted in Fig. 7. In each interval of 100 G including more than 10 points, the mean ICM values are compared to values deduced from UNNOFIT inversion (Bommier et al. 2007). Error bars of Fig. 7 correspond to standard deviations in each interval. For azimuths, the 180°-ambiguity was taken into account to compare results between  − 90 and  + 90 degrees. Intervals of 10 degrees are used in Fig. 7c, and magnetic fields smaller than 200 G are excluded as in previous sections.

Table 3 shows departures between vector magnetic field amplitudes resulting from both inversion methods in five different ranges of increasing magnetic fields. ϵ-values are mean differences and σ-values standard deviations around them. ϵ-values specify mainly systematic differences and σ-values mainly noise effects.

For weak magnetic fields, the two methods (UNNOFIT and ICM) give results that agree well, taking into account the inaccuracy resulting from the polarimetric noise.

For strong magnetic fields, the mean differences ϵ do not exceed about 50 G for vector amplitudes up to 2500 G. But the σ values are higher. They are probably partly caused by the correction of the telluric blend. Some tests have shown that the results are very sensitive to the wavelength interval used for the cubic interpolation. This might be checked by an inversion of satellite data.

In addition to blend effects, there are two other important points that have not been addressed in our ICM inversion so far:

• 1.

  instrumental profiles were not taken into account. Theyincrease the effective width of line profiles: the computed δ-valuesare between 0 and 4 across the fullfield of view, while the expected values arebetween 0 (quiet regions, FALC)and 1 (spots, MALTM);

• 2.

  limitations caused by spatial resolution and seeing effects are neglected. They also increase the line widths and disturb the magnetic field determinations by mixing solar points of different B-values. A good way to eliminate part of them should be to take into account filling-factor models.

10. Shorter ICM inversion with average model atmosphere

Because of instrumental profiles and limited spatial resolution, the differences between UNNOFIT and ICM inversions of THEMIS data given in Table 3 are larger than the expected standard deviations from noise effects, even in the case of inversion by average polynomials (P_FALC + P_MALTM)/2 (Sect. 8.6, Table 2).

It is interesting to compare the UNNOFIT inversion with the ICM one in the case of average polynomials. The obtained results are quite similar to the previous ones (Table 4). This suggests that before any more sophisticated way to take into account instrumental profiles and spatial resolution limitations is considered, this short ICM inversion with average coefficients may be useful to invert the observed data. The inversion scheme of Sect. 8.2 is reduced: in point 3, A_D and δ were not calculated, and point 4 is suppressed. Only three 7 × 7 arrays of coefficients C_1,m,n, C_2,m,n and C_3,m,n are necessary to perform this fast inversion for a given solar line.

11. Conclusion

Inversion methods based on central moments (ICM) determine vector magnetic fields with very few iterations in realistic model atmospheres.

They use three observable quantities, A₁, A₂ and Φ. A₁ and A₂ are nearly linear functions of magnetic field components in the full range of weak and strong fields. They are independent of mean Doppler velocities, almost independent of other thermodynmical parameters, and optimized to be also almost independent of data noise.

The extremly low number of iterations is possible because when a quantity a is a known function of two quantities b and c, and when a is a monotonic function of b, it is possible to expand b as a function of a and c: $\mathit{a}\mathrm{\left(}\mathit{b,c}\mathrm{\right)}\mathrm{=}\mathrm{\Rightarrow }\mathit{b}\mathrm{\left(}\mathit{a,c}\mathrm{\right)}\mathit{,}$(52)which is the case for A₁(B_∥,B_⊥) and A₂(B_⊥,B_∥).

The ICM inversions can simultaneously address a wide range of models by using additional quantities A_D. In a shorter version, they can be used with an average model atmosphere and three fixed polynomials only.

The ICM methods are very fast (typically less than 30 s for 100 000 solar points) and very suitable for large data sets. It is possible to improve them in several ways. Instrumental profiles can be included. Solar models and an inversion process that include filling factors can also significantly increase the accuracy of the inversions. The magneto-optical effects in weak magnetic fields could be determined in more sophisticated ways that are less sensitive to noise.

We must remember that ICM methods do not require the restriction on the solution that is imposed by the Milne-Eddington inversion, namely that quantities be constant in height. Therefore, ICM methods could be very successful in determining the magnetic field structures in atmospheric conditions where the Milne-Eddington condition is not met. For example, while the ICM method presented in this paper uses the exponents β₁ = 1 and β₂ = 2, the method could be extended by varying these values, which would correspond to probing the height variation of the magnetic field, because higher β-values give more weight to the line center, which forms higher than the line wings. More generally, extended ICM inversions using more parameters may provide promising new opportunities.

Finally, we note that the present method is free of any “first guess” of the solution to initiate the iteration, contrary to most usual inversion codes that are based on the Levenberg-Marquardt algorithm for minimization.

Acknowledgments

We thank the referee for fruitful comments. Han Uitenbroek gratefully acknowledges financial support of the University of Nice Sophia Antipolis during the course of this work. This collaboration has also been supported by the CNRS through a grant of its Programme National Soleil-Terre.

References

All Tables

All Figures

	Fig. 1 Functions A1, A2 and Φ (Eqs. (10), (21), (23)), measured in the case of the Fe i 630.2 nm line in the FALC model at disk center. The line styles depend on the values of B ⊥  in panels a) and c), and the values of B∥ in panel b): full lines until 1000 G, dashed lines until 2000 G, dashed-dotted until 3000 G.
In the text

	Fig. 2 Width W (Eq. (19)) in the case of the Fe i 630.2 nm line in the FALC model at disk center. Line styles are the same as in Fig. 1, panel b).
In the text

	Fig. 3 Comparison of longitudinal (panel a)) field strength, transverse field strength (panel b)), and azimuth (panel c)) recovered with the inversion with input values.
In the text

	Fig. 4 Inversion of synthetic profiles for the FALC model atmosphere at disk center with noise a level of 1.5 × 10-3Ic and a sampling of 2 pm. Line styles are the same as in Fig. 1.
In the text

	Fig. 5 Widths AD of profiles I ± V computed with two selected models for inversion addressing several models.
In the text

	Fig. 6 Optimization of exponents β1 and β2 for a low sensitivity to data noise (see Sect. 8.3).
In the text

	Fig. 7 THEMIS/MTR data: comparison between ICM and UNNOFIT results assuming a filling factor f = 1.
In the text