# Inversion of the radiative transfer equation for polarized light

Jose Carlos del Toro Iniesta
Instituto de Astrofísica de Andalucía (CSIC)
Apdo. de Correos 3004, E-18080, Granada, Spain
email<EMAIL_ADDRESS>
Basilio Ruiz Cobo
Instituto de Astrofísica de Canarias
Vía Láctea, s/n, E-38200, La Laguna (Tenerife), Spain
email<EMAIL_ADDRESS>

###### Abstract

Since the early 1970s, inversion techniques have become the most useful tool
for inferring the magnetic, dynamic, and thermodynamic properties of the solar
atmosphere. Inversions have been proposed in the literature with a sequential
increase in model complexity: astrophysical inferences depend not only on
measurements but also on the physics assumed to prevail both on the formation
of the spectral line Stokes profiles and on their detection with the
instrument. Such an intrinsic model dependence makes it necessary to formulate
specific means that include the physics in a properly quantitative way. The
core of this physics lies in the radiative transfer equation (RTE), where the
properties of the atmosphere are assumed to be known while the unknowns are
the four Stokes profiles. The solution of the (differential) RTE is known as
the direct or forward problem. From an observational point of view, the
problem is rather the opposite: the data are made up of the observed Stokes
profiles and the unknowns are the solar physical quantities. Inverting the RTE
is therefore mandatory. Indeed, the formal solution of this equation can be
considered an integral equation. The solution of such an integral equation is
called the inverse problem. Inversion techniques are automated codes aimed at
solving the inverse problem.

The foundations of inversion techniques are critically revisited with an
emphasis on making explicit the many assumptions underlying each of them. An
incremental complexity procedure is advised for the implementation in
practice. Coarse details of the profiles or coarsely sampled profiles should
be reproduced first with simple model atmospheres (with, for example, a few
physical quantities that are constant with optical depth). If the Stokes
profiles are well sampled and differences between synthetic and observed ones
are greater than the noise, then the inversion should proceed by using more
complex models (that is, models where physical quantities vary with depth or,
eventually, with more than one component). Significant improvements are
expected as well from the use of new inversion techniques that take the
spatial degradation by the instruments into account.

###### Contents

  1. 1 Introduction
  2. 2 Radiative transfer assumptions
    1. 2.1 The non-local thermodynamic equilibrium problem
    2. 2.2 The local thermodynamic equilibrium approximation
    3. 2.3 The Milne–Eddington approximation
    4. 2.4 The weak-field approximation
    5. 2.5 The MISMA hypothesis
  3. 3 Degrees of approximation in the model atmospheres
    1. 3.1 Constant physical quantities
      1. 3.1.1 The Milne–Eddington atmosphere
      2. 3.1.2 The weak-field atmosphere
      3. 3.1.3 Constant vector magnetic field and LOS velocity
    2. 3.2 Physical quantities varying with depth
      1. 3.2.1 Parameterizing the stratifications
      2. 3.2.2 The MISMA atmosphere
      3. 3.2.3 Other special atmospheres
  4. 4 Degrees of approximation in the Stokes profiles
  5. 5 A synthesis approach
    1. 5.1 Constant atmospheres
    2. 5.2 Depth-dependent atmospheres
    3. 5.3 MHD simulations
  6. 6 Response functions
    1. 6.1 Properties of response functions
    2. 6.2 Analytic response functions
  7. 7 Inversion techniques
    1. 7.1 Topology in the space of observables
    2. 7.2 Levenberg–Marquardt based inversions
      1. 7.2.1 Problems in practice
      2. 7.2.2 Automatic selection of nodes
      3. 7.2.3 A non-LTE inversion technique
    3. 7.3 Database-search inversions
    4. 7.4 Other algorithm inversions
      1. 7.4.1 Artificial neural network inversions
      2. 7.4.2 Genetic algorithm inversions
      3. 7.4.3 Bayesian inversions
    5. 7.5 Inversions accounting for spatial degradation
      1. 7.5.1 Spatially-coupled inversions
      2. 7.5.2 Regularized deconvolution inversions
      3. 7.5.3 Sparse inversions
    6. 7.6 Summary of inversion techniques
  8. 8 Discussion on inversion results
    1. 8.1 Increasing complexity in the model atmospheres
    2. 8.2 Inversion retrievals of weak fields
  9. 9 Conclusions
  10. A Appendix. Optimum practical initialization

##  1 Introduction

Unlike other branches of physics, astrophysics cannot apply the third pillar
of the scientific method, experimentation. After observing nature and
conjecturing laws that govern its behavior, astronomers cannot carry out
experiments that confirm or falsify the theory. Experimentation is then
substituted by new observations conducted to check the theoretical
predictions. The intrinsic inability for directly measuring the celestial
objects adds a special difficulty to the astrophysical tasks. We do not have
thermometers, weighing scales, tachometers, magnetometers that can directly
gauge the physical conditions in the object. Rather we have to be content with
indirect evidence or inferences obtained from the only real astrophysical
measurements, namely those related to light. The intensity and polarization
properties for visible light, the associated electric field for radio
frequencies, or the energy or momentum of high energy photons, as functions of
space, wavelength, and time, can be fully quantified with _errors_ that are
directly related to the accuracy of the instruments.111Actual experiments have
also been devised to directly detect neutrinos, other particles, and
gravitational waves, but such projects are way beyond the scope of this paper.
From these real measurements, the observational astronomer must deduce or
infer the physical quantities that characterize the object with
_uncertainties_ that depend on both the experimental errors and the
assumptions that allow him/her to translate light-derived quantities into the
object quantities. Observational astrophysics could hence probably be defined
as the art of _inferring_ the physical quantities of heavenly bodies from real
measurements of the light received from them.

Somehow, these astrophysical tasks can be mathematically seen as a mapping
between two spaces, namely the space of observables and that of the object’s
physical quantities. The success of the astronomer then depends on his/her
ability (_the art_) to characterize not only the mapping but the two spaces.
On the observable side, what really matters is the specific choice of
measurable parameters and how well they are measured; that is, how many light
parameters are obtained (the signal) and which are the measurement errors (the
noise). On the object’s physical condition side, what is substantive is the
selection of quantities to be inferred. Of course, the finer the —affordable—
detail in describing any of the two spaces, the better. The keyword is
_affordable_ because infinite resolution does not exist in the real world: a
compromise is always in order between the number of available observables and
the number of inferred physical quantities. The representation of both spaces
therefore needs approximations that constrain the sub-spaces to be explored
and how they are described: which Stokes parameters, with which wavelength and
time sampling, and with which instrument profile and resolution on the one
hand, and, on the other, which quantities and how they are assumed to vary on
the object with time and space. Concerning the mapping, this should represent
the physics that generates the observables from the given physical conditions
in the object and thus illustrates the dependence of the observables on given
physical quantities. Understanding this physics is crucial if the researcher
is to select observables that are as “orthogonal” as possible; that is, that
depend mostly on one physical quantity and not on the others. Certainly, the
physics mapping needs approximations as well. These approximations depend a
great deal on the observables and on the object’s physical quantities; for
example, the assumptions cannot be the same if you have fully sampled Stokes
profiles or just a few wavelength samples; different hypotheses apply for
physical quantities that do or do not vary with depth in the atmosphere, or
that are expected to present a given range of magnitudes. Therefore, mappings
may include (often over-simplistic) one-dimensional calibration curves between
a given observable parameter and a given physical quantity, or complicated
multidimensional relationships between observables and quantities that require
the definition of a metric or distance in at least one of the two spaces.

Even in the simplest situations, the relationship between observables and
quantities does not have to be linear and may depend on the specific sub-space
of the physical parameters. For example, a calibration curve based on the
weak-field approximation may apply for a given range of magnetic fields but
saturate for stronger ones (see Sections 2.4 and 3.1.2). But, when the problem
can be assumed to be multidimensional, covariances appear because single
observables rarely depend on just a single quantity (see Section 7). For
example, a given spectral line Stokes $V$ profile can seemingly grow or weaken
by the same amount owing to changes in temperature or magnetic field strength
(e.g., Del Toro Iniesta and Ruiz Cobo, 1996). An example can be seen in Figure
1, where two apparently equal $V$ profiles come from two different
atmospheres. With all these ingredients at hand, the astrophysical analysis of
observations is a non-linear, fully involved, topological task where many
decisions have to be made (_the art_) and, hence, cannot be taken for granted.

Figure 1: Left panel: Open circles: Stokes $V$ profile in units of the
continuum intensity of the Fe i line at 630.25 nm synthesized in a model
atmosphere in hydrostatic equilibrium, 2000 K cooler than the Del Toro Iniesta
et al. (1994) model, with a constant longitudinal magnetic field of 800 G, a
gradient in velocity from 2  km s-1at the bottom and 0  km s-1at the top of
the photosphere, and a macroturbulence velocity of 1  km s-1. Solid line:
Stokes $V$ profile of the same line (normalized the same way), synthesized in
a model atmosphere 305 K hotter that the former, 270 G weaker, and with a
higher macroturbulence velocity of 2.06  km s-1. Right panel: $T$ and $B$
stratifications for the two models.

The techniques by which astronomers have obtained information about the
physical conditions in the object have evolved in parallel to technological
advancements; that is, to the available means we have of gathering such
information. The community has gradually enhanced its knowledge from medium-
band measurements including one or several spectral lines to very fine
wavelength sampling of the four Stokes profiles of single or multiple spectral
lines; from old curves of growth for equivalent widths to highly sophisticated
techniques that include the solution of the radiative transfer equation (RTE).
The finer the information, the more complete the physical description.

Following Socas-Navarro (2001), let us consider the simplest case of having a
single observable parameter, the Doppler displacement with respect to the rest
position of the spectral line, $\Delta\lambda$, and a single physical quantity
to derive, the line-of-sight (LOS) velocity, $v_{\rm LOS}$. Imagine that we
measure $\Delta\lambda$ by finding the minimum (or the maximum in the case of
an emission line) of the intensity profile. The biunivocal mapping between the
one-dimensional space of observables —that containing all possible Doppler
displacements— and the one-dimensional space of physical quantities —that of
LOS velocities— is given by the Doppler formula

$v_{\rm LOS}=\frac{\Delta\lambda}{\lambda_{0}}c,$ (1)

where $\lambda_{0}$ stands for the vacuum rest wavelength position of the line
and $c$ for the speed of light. This simple inference relationship requires at
least three implicit physical assumptions for the Doppler displacement to be
properly defined and measured; namely that a) the solar feature is spatially
resolved, b) the line is in pure absorption (or pure emission), and c) $v_{\rm
LOS}$ is constant along the LOS. First, if we have unresolved structures we
cannot ascribe the inferred velocity to any of them. Second, lines with core
reversals, either in absorption or in emission, do not qualify for the
extremum-finding method. And third, as soon as we have an asymmetric profile,
$\Delta\lambda$ can no longer be properly defined for the line but for a given
height through the profile, and then the mapping in Eq. (1) immediately loses
its meaning. While in the case of a constant velocity, we properly infer that
velocity, in the presence of gradients we infer a value corresponding only to
the —in principle unknown— layers where the core of our line has been formed
(typically the highest layers of the atmosphere). _We measure a velocity but
we do not know which one._ Strictly speaking, the same measurement corresponds
to different physical quantities depending on the assumptions. Of course we
could complicate our problem a little and try to determine the stratification
of LOS velocities with height, or simply estimate a gradient, by measuring the
so-called bisector, the geometric position of those points equidistant from
both wings of the profile at a given depth. At that point, our spaces have
increased their dimensions and Eq. (1) is no longer the sole ingredient of our
mapping because we must add some more physical assumptions to interpret the
different displacements of the bisector in terms of velocities at different
heights in the atmosphere. Hence, depending on the assumed physics, the
quantitative results may change. This easy example has been used to illustrate
that even the simplest inference is dependent on physical assumptions. This is
an inherent property of astrophysical measurement and no one can escape from
it: the same observable can mean different things depending on the assumed
underlying physics. Most of the criticisms of the inversion techniques that
are reviewed in this paper often come from this lack of uniqueness of the
results. Many authors claim that the inversion of the RTE is an ill-posed
problem. This being true, one should realize that astrophysics itself is
indeed ill-conditioned, and this is a fact we have to deal with, either
willingly or not.

The physics connecting the object quantities with the observable parameters is
of paramount significance and deserves a little consideration at this point.
Radiative transfer is the discipline encompassing the generation and transport
of electromagnetic radiation through the solar (stellar) atmosphere. Hence,
the mapping between the two spaces will be based upon it and depend on its
degrees of approximation. The specification of the radiation field through a
scattering atmosphere was first formulated as a physical problem by Strutt
(1871a) (Lord Rayleigh); Strutt (1871b) (Lord Rayleigh); Strutt (1881) (Lord
Rayleigh); Strutt (1899) (Lord Rayleigh). In the astrophysical realm, the
problem was posed in the works by Schuster (1905) and Schwarzschild (1906)
without taking polarization into account. After that, although not known to
the astrophysical community, Soleillet (1929) presented a theory of
anisotropic absorption that is nothing but a rigorous formulation of the
radiative transfer equation. Very importantly, he used the formalism proposed
by Stokes (1852) to deal with partially polarized light. It was not, however,
until the works by Chandrasekhar (1946a, b, 1947) that the transfer problem of
polarized light was settled as an astrophysical problem on its own. The Stokes
formalism has regularly been used since then in the astronomical literature.
After Hale’s (1908) discovery of sunspot magnetic fields, the interpretation
of the solar (stellar) spectrum of polarized light became necessary and a full
theory has been developed since the mid 1950s. The first modern formulation of
an equation of radiative transfer for polarized light was presented by Unno
(1956), who also provided a solution in the simplified case of a
Milne–Eddington (ME) atmosphere. Only absorption processes were taken into
account and a complete description had to wait until the works by Rachkovsky
(1962a, b, 1967), who also included dispersion effects (the so-called magneto-
optical effects). These two derivations were phenomenological and somewhat
heuristic. A rigorous derivation of the radiative transfer equation (RTE)
based on quantum electrodynamics was obtained by Landi degl’Innocenti and
Landi degl’Innocenti (1972). Later, four derivations of the RTE from basic
principles of classical physics were published by Jefferies et al. (1989),
Stenflo (1991; see also Stenflo 1994), Landi Degl’Innocenti (1992; see also
Landi Degl’Innocenti and Landolfi 2004), and Del Toro Iniesta (2003b). A
discussion of the RTE and the several assumptions used in various available
inference techniques is deferred to Section 2.

Certainly, any inference has to be based on solutions of the RTE because it
relates the observable Stokes spectrum with the unknowns of the problem;
namely, the physical quantities characterizing the state of the atmosphere
they come from. No matter how simplified such solutions can be, it is natural
to compare the observations with theoretical calculations in prescribed sets
of physical quantities. The comparison of observational and synthetic
parameters results in values for the sought-for quantities that may be refined
in further iterations by changing the theoretical prescriptions. This trial-
and-error method can be practical when the problem is very simple (involving a
few free parameters) but can become unsuitable for practical use if the number
of free parameters is large. Even automated trial-and-error —i.e., Monte
Carlo— methods may fail to converge to a reliable set of physical conditions
in the medium. Some more educated techniques are needed to finally work out
that convergence between observed and synthetic parameters.

Generally speaking, any method in which information about the integrand of an
integral equation is obtained from the resulting value of the integral is
called an inversion method. In our particular case, it is straightforward to
write the synthetic Stokes spectra as an integral involving a kernel that
depends on the physical conditions of the atmosphere (see Eq. 8). In fact, the
emergent formal solution of the RTE is the most basic type of integral
equation, namely a Fredholm equation of the first type, because both
integration limits are fixed. Consequently, we will call inversion codes or
inversion techniques those methods that (almost) automatically succeed in
finding reliable physical quantities from a set of observed Stokes spectra
because we shall understand that they indeed automatically solve that integral
equation. There is a whole variety of flavors depending on the several
hypotheses that can be assumed, but all of them share the characteristic
feature of automatically minimizing a distance in the topological space of
observables. The idea had already been clearly explained in the seminal work
by Harvey et al. (1972): “Solve for $B$ on the bases of best fit of the
observed profiles to the theoretical profiles”. And the free parameters for
such a best fit were found through least squares minimization of the profile
differences. They obtained only an average longitudinal field component
because their Stokes $Q$ and $U$ observations were not fully reliable and
magneto-optical effects were not taken into account, but the fundamental idea
underlying many of the current techniques can already be found in that very
paper, including a simple two-component model to describe the possible
existence of spatially unresolved magnetic fields.

In a thorough study using synthetic Stokes profiles, Auer et al. (1977)
proposed a new inversion method based on Unno’s theory and tested its behavior
in the presence of several realistic circumstances, such as asymmetric
profiles, magnetic field gradients, magneto-optical effects, and unresolved
magnetic features. This technique was later generalized by Landolfi et al.
(1984) to include magneto-optical and damping effects. The numerical check of
the code was fairly successful but neither the original code by Auer et al.
(1977) nor the new one by Landolfi et al. (1984) were applied to observations.
Independently of the latter authors, the preliminary studies by Skumanich and
Lites (1985), Lites and Skumanich (1985), and Skumanich et al. (1985) jelled
in what has been one of the most successful ME inversion codes so far by
Skumanich and Lites (1987), later extended by Lites et al. (1988) to mimic a
chromospheric rise in the source function (see Section 2.3). This code has
been extensively used with observational data, most notably those obtained
with the Advanced Stokes Polarimeter (Elmore et al., 1992).

Based on the thin flux tube approximation, Keller et al. (1990) proposed an
inversion code for extracting physical information not from the Stokes
profiles themselves but from several parameters calculated from $I$ and $V$
observations of a plage and a network. Two years later, Solanki et al. (1992a)
presented a new inversion code whereby from the whole Stokes $I$ and $V$
profiles they selected among a handful of prescribed temperature
stratifications and inferred height-independent magnetic field strength and
inclination, Doppler shift, filling factor (surface fraction in the resolution
element covered by magnetic fields), macro- and microturbulent velocities, and
some atomic parameters of the spectral line. The very same year, Ruiz Cobo and
Del Toro Iniesta (1992) introduced SIR, an acronym for Stokes Inversion based
on Response functions. Like the former codes, SIR ran a non-linear, least-
squares, iterative Levenberg–Marquardt algorithm but with a remarkable step-
forward feature: physical quantities characterizing the atmosphere were
allowed to vary with optical depth. The increase of free parameters can
generate a singularity problem: the variation of some atmospheric parameters
may not produce any change on the synthetic spectra or, in other cases,
different combinations of the perturbation of several parameters may produce
the same change in the spectra. The success of SIR lies in regularizing the
problem through a _tailored_ Singular Value Decomposition method (SVD). This
allows, in principle, to look for any arbitrarily complex atmospheric
stratification. The three components of the magnetic field, the LOS velocity,
the temperature stratification, and the microturbulence may have any height
profile. The code also infers height-independent microturbulent velocity and
filling factor. The possibility exists for also fitting some atomic parameters
(e.g., Allende Prieto et al., 2001). but they are typically fixed in practice.
The code can be applied to any number of spectral lines that are observed
simultaneously. SIR has been successful in a large number of observing cases
and its use is still spreading among the community.

Following SIR’s strategy (that is, using response functions, nodes,
Levenberg–Marquardt, and SVD), an evolution of the Solanki et al. (1992a) code
called SPINOR was presented by Frutiger and Solanki (1998) that also allowed
for height variations of the physical quantities and included the possibility
of multi-ray calculations assuming the thin flux tube approximation. Sánchez
Almeida (1997) proposed an original inversion code under the MISMA (MIcro-
Structured Magnetic Atmosphere) hypothesis (see Sections 2.5 and 3.2.2). In
2000, the codes by Socas-Navarro et al. (2000, NICOLE —NLTE Inversion Code
based on the Lorien Engine—) and by Bellot Rubio et al. (2000; see also 1997)
were presented. The first (based on an earlier code by Socas-Navarro et al.
1998 without taking either polarization or magnetic fields into account)
included non-LTE radiative transfer (see Section 2.1), and the second was
specifically designed for analyzing Stokes $I$ and $V$ profiles in terms of
the thin flux tube approximation by using an analytic shortcut for radiative
transfer proposed by Del Toro Iniesta et al. (1995, see Section 3.2.3). On
their hand, Rees et al. (2000) proposed a Principal Component Analysis (PCA),
which worked by creating a database of synthetic Stokes profiles by means of
an SVD technique. In such a database, given _eigenprofiles_ are obtained that
are later used as a basis for expanding the observed Stokes profiles. Hence,
the description of observations can be made with the help of a few
coefficients, thus speeding up the inversion process. One year later, LILIA
(LTE Inversion based on Lorien Iterative Algorithm), a code with similar
properties as SIR, was presented by Socas-Navarro (2001) and FATIMA (Fast
Analysis Technique for the Inversion of Magnetic Atmospheres), a PCA code, was
introduced by Socas-Navarro et al. (2001). A different technique was proposed
by Carroll et al. (2001, see also  ) that used artificial neural networks
(ANNs) whereby the system was trained with a set of synthetic Stokes profiles.
The structure obtained therefrom finds the solution for the free parameters by
interpolating among the known ones. Although the training can be slow, the
inversion of observational data is very fast. In practice, both the synthetic
training set of ANNs and the synthetic database of PCA have employed ME
profiles to keep the implementation feasible. Otherwise, the number of free
parameters would render the two techniques impracticable. A PCA code to
analyze the Hanle effect in the He i D3 line was developed by López Ariste and
Casini (2003, see also , ).

A substantial modification of the original SIR code, called SIRGAUSS, was
presented by Bellot Rubio (2003) in which the physical scenario included the
coexistence of an inclined flux tube —that is pierced twice by the LOS— within
a background. Such a scenario is used to describe an uncombed field model of
sunspot penumbrae (Solanki and Montavon, 1993). An evolution of this inversion
code, called SIRJUMP, was later used by Louis et al. (2009) that was able to
infer possible discontinuities in the physical quantities along the LOS. A
further code presented by Asensio Ramos (2004) was able to deal with the
Zeeman effect in molecular lines. The very same year, Lagg et al. (2004)
published HeLIx, an ME inversion code that dealt with the Hanle and the Zeeman
effect in the He i line at 1083 nm. Another ME inversion code was presented by
Orozco Suárez and Del Toro Iniesta (2007) with the helpful feature that was
written in IDL, so that it is easily manipulated by relatively inexperienced
users and employed as a routine in high-level programming pipelines. Also in
2007, Bommier et al. took over the Landolfi et al. (1984) method and extended
it to include unresolved magnetic structures. Unfortunately, they fail to
obtain the magnetic field strength and the filling factor separately; only
their product is reliable. Self-consistent levels of confidence in the ME
inversion results were estimated through the code proposed by Asensio Ramos et
al. (2007a) using Bayesian techniques. A rigorous treatment of optical
pumping, atomic level polarization, level crossings and repulsions, Zeeman,
Paschen-Back, and Hanle effects on a magnetized slab was included in HAZEL
(Asensio Ramos et al., 2008), with which analysis of the He i D3 and the
multiplet at 1083 nm can be carried out.

Oriented to its extensive use with the data coming from the Helioseismic and
Magnetic Imager (Graham et al., 2003) aboard the Solar Dynamics Observatory,
Borrero et al. (2011) presented VFISV (Very Fast Inversion of the Stokes
Vector), a new ME code but with several further approximations and simplifying
assumptions to make it significantly faster than other available codes. Mein
et al. (2011) presented an alternative inversion code in which, with a
significant number of simplifying assumptions on top of the ME approximation
(such Stokes $I$ profiles being Gaussians and magneto-optical effects being
almost negligible), some moments of the Stokes profiles are used to retrieve
the vector magnetic field and the LOS velocity. In 2012, a significant step
forward was provided by van Noort, who combined spectral information with the
known spatial degradation effects on two-dimensional maps to obtain a
consistent restoration of the atmosphere across the whole field of view. An
aim similar to van Noort’s is followed by Ruiz Cobo and Asensio Ramos (2013),
who, by means of a regularized method (indeed based on PCA), deconvolve the
spectropolarimetric data that are later inverted with SIR. Based on the
concept of sparsity, Asensio Ramos and de la Cruz Rodríguez (2015) have
proposed a novel technique that allows the inversion of two-dimensional
(potentially three-dimensional) maps at once.

The interested reader can complement this chronological overview with the
reviews by Del Toro Iniesta and Ruiz Cobo (1995, 1996, 1997), Socas-Navarro
(2001), Del Toro Iniesta (2003a), Bellot Rubio (2006), and Asensio Ramos et
al. (2012) and the didactical introductions and discussions by Stenflo (1994),
Del Toro Iniesta (2003b), and Landi Degl’Innocenti and Landolfi (2004). A
critical discussion on the different techniques and the specific
implementations will be developed throught the paper, which is structured as
follows: the basic assumptions of radiative transfer are discussed in Section
2; the following two Sections discuss the approximations used for the model
atmospheres and the Stokes profiles; an analysis of the forward problem,
namely the synthesis of the Stokes spectrum, is presented in Section 5, which
is followed by an analysis of the sensitivities of spectral lines to physical
quantities (Section 6); the basics of inversion techniques are analyzed in
Section 7 and a discussion on inversion results presented in Section 8;
finally, Section 9 summarizes the conclusions. An appendix proposes an optimum
way of initializing the inversion codes through the use of classical
estimates.

##  2 Radiative transfer assumptions

The propagation of electromagnetic energy through a stellar atmosphere —and
its eventual release from it— is a significantly complex, non-linear, three-
dimensional, and time-dependent problem where the properties of the whole
atmosphere are involved. From deep layers up to the stellar surface, the
coupling between the radiation field and the atmospheric matter implies non-
local effects that can connect different parts of the atmosphere. In other
words, the state of matter and radiation at a given depth may depend on that
at the other layers: light emitted at one point can be absorbed or scattered
at another to release part or all of its energy.

The description of the whole system, matter plus radiation field, needs to
resort to the solution of the coupled equations that describe the physical
state of the atomic system and that of the radiation traveling through it.
Therefore, we have to simultaneously solve the so-called statistical
equilibrium equations and the radiative transfer equation. The first
assumption we shall make is that radiative transfer is one dimensional; that
is, that the transfer of radiative energy perpendicular to the line of sight
can be neglected in the matter–radiation coupling. For most solar applications
so far, this assumption has been seen to be valid. Since the purpose of this
paper is not directly related to either of the two systems of equations, let
us simply point out what their main characteristics and ingredients are, and
how the whole problem can be simplified in different situations. We refer the
interested reader to the book by Landi Degl’Innocenti and Landolfi (2004) for
a full and rigorous account of all the details.

Most classical radiative transfer descriptions in the literature do not deal
with polarization. They are typically qualified as radiative transfer studies
for unpolarized light but the name is ill-chosen. Formally speaking, those
analyses are for light traveling through homogeneous and isotropic media (Del
Toro Iniesta, 2003b). As a consequence of that heritage, the community is used
to speak about atomic level populations either calculated through the
Boltzmann and Saha equations (the LTE approximation; see Sect. 2.2) or not
(the non-LTE case; see Section 2.1). These isotropic descriptions of the
transfer problem, however, are not valid when a physical agent such as a
vector magnetic field establishes a preferential direction in the medium,
hence breaking the isotropy. Moreover, the outer layers of a star are a clear
source of symmetry breaking. The exponential density decrease with height
makes the radiation field anisotropic: outward opacity is much smaller than
inward opacity. This should also be the case with collisions between
particles: they are more probable at the bottom than at the top of the
atmosphere. In such a situation, the probability is not zero for the various
degenerate levels of the atom (with respect to energy) to be not evenly
populated and for non-zero coherences or phase relations between them to
exist. The atomic system is then said to be polarized and its state is best
described with the so-called density operator, $\rho$, that provides the
probabilities of the sublevels being populated (hence the populations) along
with the possible _correlations_ or _interferences_ between every pair. In the
standard representation that uses the eigenvectors of the total angular
momentum, $\mbox{\boldmath$J$}^{2}$, and of its third component,
$\mbox{\boldmath$J$}_{z}$, as a basis, the density matrix element

$\rho(\alpha jm,\alpha^{\prime}j^{\prime}m^{\prime})=\langle\alpha
jm|\rho|\alpha^{\prime}j^{\prime}m^{\prime}\rangle$ (2)

represents the coherence or phase interference between the different magnetic
sublevels characterized by their angular momentum quantum numbers. In Eq. (2),
$\alpha$ and $\alpha^{\prime}$ stand for supplementary quantum numbers
relative to those operators that commute with $\mbox{\boldmath$J$}^{2}$ and
$\mbox{\boldmath$J$}_{z}$. Certainly, the diagonal matrix elements
$\rho_{\alpha}(jm,jm)\equiv\rho(\alpha jm,\alpha jm)$ represent the
populations of the magnetic sublevels and the sum

$n_{j}=\sum_{m}\rho_{\alpha}(jm,jm)=\sum_{m=-j}^{j}\langle\alpha
jm|\rho|\alpha jm\rangle$ (3)

accounts for the total population of the level characterized by the $j$
quantum number.

At all depths in the atmosphere, evolution equations for these density matrix
elements have to be formulated that describe their time ($t$) variations due
to the transport of radiation, on the one hand, and to collisions among
particles on the other. All interactions with light —namely, pure absorption
($A$), spontaneous emission ($E$), and stimulated emission ($S$)— have to be
considered. All kinds of collisions —namely, inelastic ($I$), superelastic
($S$), and elastic ($E$) collisions— have to be taken into account. Inelastic
collisions induce transitions between any level $|\alpha jm\rangle$ and an
upper level $|\alpha_{u}j_{u}m_{u}\rangle$ with a consequent loss in kinetic
energy. Superelastic collisions induce transitions to a lower energy level
$|\alpha_{l}j_{l}m_{l}\rangle$ with an increase in the kinetic energy of
collision. Finally, elastic collisions induce transitions between degenerate
levels $|\alpha jm\rangle$ and $|\alpha jm^{\prime}\rangle$; in these, the
colliding particle keeps its energy during the interaction. The statistical
equilibrium equations (4) and (5) that follow for radiative and collisional
interactions, respectively, have slightly different application ranges. The
former are valid for the multi-term atom representation and can even be used
in the Paschen–Back regime, while the latter are only valid for the special
case of the multi-level atom representation (although they can be generalized
to the multi-term representation).222The concepts of multi-level or multi-term
representation of an atomic system basically depend on the assumption or not,
respectively, that coherences can be neglected among magnetic sub-levels that
belong to levels characterized by different quantum numbers $\alpha$ and $j$.
See Landi Degl’Innocenti and Landolfi (2004) for a detailed and rigorous
description. We make them explicit here for illustrative purposes only and
refer the interested reader to the Landi Degl’Innocenti and Landolfi’s (2004)
monograph for details. According to that work, the radiative interaction
equations in the magnetic field reference frame333Where the vector magnetic
field marks the $Z$ direction. can be written as

${\displaystyle\frac{\rm d}{{\rm
d}t}}\rho_{\alpha}(jm,j^{\prime}m^{\prime})=-2\pi{\rm
i}\,\nu_{\alpha}(jm,j^{\prime}m^{\prime})\rho_{\alpha}(jm,j^{\prime}m^{\prime})+$
${\displaystyle\sum_{\alpha_{l}j_{l}m_{l}j^{\prime}_{l}m^{\prime}_{l}}^{\,}}\rho_{\alpha_{l}}(j_{l}m_{l},j^{\prime}_{l}m^{\prime}_{l})\,T_{A}(\alpha
jmj^{\prime}m^{\prime},\alpha_{l}j_{l}m_{l}j^{\prime}_{l}m^{\prime}_{l})+$
${\displaystyle\sum_{\alpha_{u}j_{u}m_{u}j^{\prime}_{u}m^{\prime}_{u}}^{\,}}\rho_{\alpha_{u}}(j_{u}m_{u},j^{\prime}_{u}m^{\prime}_{u})\left[T_{E}(\alpha
jmj^{\prime}m^{\prime},\alpha_{u}j_{u}m_{u}j^{\prime}_{u}m^{\prime}_{u})+T_{S}(\alpha
jmj^{\prime}m^{\prime},\alpha_{u}j_{u}m_{u}j^{\prime}_{u}m^{\prime}_{u})\right]-$
(4)
${\displaystyle\sum_{j^{\prime\prime}m^{\prime\prime}}^{\,}}\left\\{\rho_{\alpha}(jm,j^{\prime\prime}m^{\prime\prime})\left[R_{A}(\alpha
j^{\prime}m^{\prime}j^{\prime\prime}m^{\prime\prime})+R_{E}(\alpha
j^{\prime\prime}m^{\prime\prime}j^{\prime}m^{\prime})+R_{S}(\alpha
j^{\prime\prime}m^{\prime\prime}j^{\prime}m^{\prime})\right]\right.+$
$\rho_{\alpha}(j^{\prime\prime}m^{\prime\prime},j^{\prime}m^{\prime})[R_{A}(\alpha
j^{\prime\prime}m^{\prime\prime}jm)+R_{E}(\alpha
jmj^{\prime\prime}m^{\prime\prime})+R_{S}(\alpha
jmj^{\prime\prime}m^{\prime\prime})]\\}$

where $\nu_{\alpha}(jm,j^{\prime}m^{\prime})$ is the frequency difference
between the two sublevels and the $T$’s and $R$’s are radiative rates of
coherence transfer and relaxation among the sublevels, respectively. Now, the
collisional interactions give

${\displaystyle\frac{\rm d}{{\rm
d}t}}\rho_{\alpha}(jm,jm^{\prime})={\displaystyle\sum_{\alpha_{l}j_{l}m_{l}m^{\prime}_{l}}^{\,}}C_{I}(\alpha
jmm^{\prime},\alpha_{l}j_{l}m_{l}m^{\prime}_{l})\rho_{\alpha_{l}}(j_{l}m_{l},j_{l}m^{\prime}_{l})+$
${\displaystyle\sum_{\alpha_{u}j_{u}m_{u}m^{\prime}_{u}}^{\,}}C_{S}(\alpha
jmm^{\prime},\alpha_{u}j_{u}m_{u}m^{\prime}_{u})\rho_{\alpha_{u}}(j_{u}m_{u},j_{u}m^{\prime}_{u})+$
${\displaystyle\sum_{m^{\prime\prime}m^{\prime\prime\prime}}^{\,}}C_{E}(\alpha
jmm^{\prime},\alpha
jm^{\prime\prime}m^{\prime\prime\prime})\rho_{\alpha}(jm^{\prime\prime},jm^{\prime\prime\prime})-$
(5) ${\displaystyle\sum_{m^{\prime\prime}}^{\,}}\left[\frac{1}{2}X(\alpha
jmm^{\prime}m^{\prime\prime})\rho_{\alpha}(jm,jm^{\prime\prime})+\frac{1}{2}X(\alpha
jm^{\prime}mm^{\prime\prime})^{*}\rho_{\alpha}(jm^{\prime\prime},jm^{\prime})\right.-$
$\left.\frac{1}{2}X_{E}(\alpha
jmm^{\prime}m^{\prime\prime})\rho_{\alpha}(jm,jm^{\prime\prime})+\frac{1}{2}X_{E}(\alpha
jm^{\prime}mm^{\prime\prime})^{*}\rho_{\alpha}(jm^{\prime\prime},jm^{\prime})\right],$

where the $C$’s are collisional transfer rates between levels and the $X$’s
are relaxation rates. The indices refer to the corresponding type of
collisions and the asterisk denotes the complex conjugate.

With the standard notation for the Stokes pseudo-vector
$\mbox{\boldmath$I$}\equiv(I,Q,U,V)^{\scriptscriptstyle{\rm T}}$, where index
T stands for the transpose, the radiative transfer equation can be written as
(e.g., Del Toro Iniesta, 2003b)

$\frac{\rm d\mbox{\boldmath$I$}}{\rm d\tau_{\rm c}}=\mbox{\bf
K}(\mbox{\boldmath$I$}-\mbox{\boldmath$S$}),$ (6)

where $\tau_{\rm c}$ is the optical depth at the continuum wavelength, K
stands for the propagation matrix, and $S$ is the so-called source function
vector. Since the continuum spectrum of radiation can safely be assumed flat
within the wavelength span of a spectral line and non-polarized as far as
currently reachable polarimetric accuracies are concerned, the optical depth,
defined as

$\tau_{\rm c}\equiv\int_{s}^{s_{\rm lim}}\chi_{\rm cont}\,{\rm d}s,$ (7)

is the natural length scale for radiative transfer. Note that the origin of
optical depth ($\tau_{\rm c}=0$) coincides with the outermost boundary of
geometrical distances ($s_{\rm lim}$) and is taken where the observer is
located so that $\tau_{\rm c}$’s are actual depths in the atmosphere. In Eq.
(7), $\chi_{\rm cont}$ is the continuum absorption coefficient (the fraction
of incoming electromagnetic energy withdrawn from the radiation field per unit
of length through continuum formation processes). The propagation matrix deals
with _absorption_ (withdrawal of the same amount of energy from all
polarization states), _pleochroism_ (differential absorption for the various
polarization states), and _dispersion_ (transfer among the various
polarization states). The product of K and $S$ accounts for emission. The RTE
can then be considered _as a conservation equation_ : the energy and
polarization state of light at a given point in the atmosphere can only vary
because of emission, absorption, pleochroism, and dispersion. Equation (6) is
strictly valid only under the assumption that the energy and polarization
state of light are independent of time. To be more specific, we have assumed
that the rate of change of the Stokes parameter profiles is much slower than
the radiative and collisional relaxation time scales involved in the problem.

A formal solution to the general RTE was proposed for the first time by Landi
Degl’Innocenti and Landi Degl’Innocenti (1985), according to whom, the
observed Stokes profiles at the observer’s optical depth ($\tau_{\rm c}=0$)
read

$\mbox{\boldmath$I$}(0)=\int_{0}^{\infty}\mbox{\bf O}(0,\tau_{\rm c})\mbox{\bf
K}(\tau_{\rm c})\mbox{\boldmath$S$}(\tau_{\rm c})\rm d\tau_{\rm c},$ (8)

where O is the so-called evolution operator, and a semi-infinite atmosphere
has been assumed as usual. The solution is called _formal_ because it is not a
_real_ solution as long as the evolution operator (and the propagation matrix
and the source function vector) are not known. Unfortunately, no easy
analytical expression can in general be found for O. Only in some particular
cases, such as that in Sect. 2.3, can a compact form for the evolution
operator and an analytic solution of the RTE be obtained. In all other cases,
numerical evaluations of O and solutions of the transfer equation are
necessary. The emergent Stokes spectrum is obtained through an integral of a
product of three terms all over the whole atmosphere. Claiming that some of
the Stokes parameters are proportional to one of the matrix elements of K is,
at the very least, adventurous. This proportionality can only take place in
very special circumstances (e.g., Sections 2.4 and 3.1.2).

###  2.1 The non-local thermodynamic equilibrium problem

Being a vector differential equation, the RTE should indeed be considered as a
set of four _coupled_ differential equations. These can only be solved
independently in specific media, either isotropic or very simplified ones. But
the situation is far more complicated since both K and $S$ depend on the
material properties described by $\rho_{\alpha}(jm,j^{\prime}m^{\prime})$, as
well as on external fields such as a macroscopic velocity or a magnetic field.
For their part, the radiative and collisional transfer and relaxation rates do
depend on the radiation field. Therefore, Eqs. (4), (5), and (6) describe a
very involved, non-local, non-linear problem, known as the _non-local
thermodynamic equilibrium_ (NLTE) problem and must be consistently solved
altogether. The numerical solution of all those coupled equations requires
iterative procedures that are summarized in Figure 2.

Figure 2: Block diagram of the Stokes profile synthesis under NLTE conditions

By a model atmosphere we understand the set of thermodynamic variables
(usually two, e.g., temperature and pressure, $T$ and $p$), dynamic (the
macroscopic, bulk line-of-sight velocity field, $v_{\rm LOS}$), magnetic (the
vector field $B$, represented by $B$, the strength, $\gamma$, the inclination
with respect to the LOS, and $\varphi$, the azimuth), and possibly some other,
_ad hoc_ variables (such as the micro- and macroturbulence velocities,
$\xi_{\rm mic}$ and $\xi_{mac}$, the filling factor, $f$ —the area fraction of
the resolution pixel that is filled with the unknown atmosphere— and so
forth). All these variables have to be specified as functions of the optical
depth. Numerically, that model can be represented by a vector $x$ of $np+r$
components, $n$ being the number of depth grid points throughout the
atmosphere, $p$ the number of physical quantities varying with depth, and $r$
the number of quantities that are assumed constant throughout the LOS. For
example, one such model atmosphere would look like

$\begin{array}[]{rcl}\mbox{\boldmath$x$}&\equiv&[T(\tau_{1}),T(\tau_{2}),\ldots,T(\tau_{n}),p(\tau_{1}),p(\tau_{2}),\ldots,p(\tau_{n}),B(\tau_{1}),B(\tau_{2}),\ldots,B(\tau_{n}),\\\
&&\gamma(\tau_{1}),\gamma(\tau_{2}),\ldots,\gamma(\tau_{n}),\varphi(\tau_{1}),\varphi(\tau_{2}),\ldots,\varphi(\tau_{n}),\\\
&&v_{\rm LOS}(\tau_{1}),v_{\rm LOS}(\tau_{2}),\ldots,v_{\rm
LOS}(\tau_{n}),\xi_{\rm mic},\xi_{\rm mac},f]^{\scriptscriptstyle{\rm
T}},\end{array}$ (9)

where we have assumed specifically that both micro- and macroturbulence (as
well as the filling factor) are constant with depth. This assumption is based
on the fact that experience teaches that the increase in spatial resolution
reached with new instruments makes less and less necessary the use of such _ad
hoc_ parameters.

Once this model atmosphere is set, the necessary ingredients for the RTE and
the statistical equilibrium equations can be calculated. The solution of the
RTE has to be compared with that coming from it after modification driven by
the new density matrix elements resulting from the solution of the statistical
equations. If the differences are considered small compared with a given
threshold, then a new synthetic set of Stokes parameters has been found. If
not, the equilibrium equations have to be modified in order to iterate the
procedure until convergence is reached. The direct problem of obtaining the
Stokes spectrum of a given line coming out from a given model atmosphere then
turns out to be very complex. It cannot always be computed with the necessary
speed and accuracy. Approximations are, thus, in order.

###  2.2 The local thermodynamic equilibrium approximation

Imagine now that coherences among the Zeeman sublevels can be neglected, and
that all of them are evenly populated. That is, assume that

$\rho(\alpha
jm,\alpha^{\prime}j^{\prime}m^{\prime})=\delta_{\alpha\alpha^{\prime}}\delta_{jj^{\prime}}\delta_{mm^{\prime}}\rho_{\alpha
j},$ (10)

where $\delta$ is Kronecker’s delta. In such conditions,
$n_{j}=(2j+1)\rho_{\alpha j}$. Assume also that $n_{j}$ and the population of
other ionic species can be evaluated through the equations of thermodynamic
equilibrium at the local temperature (the Boltzmann and Saha laws; e.g., Gray,
2005). This assumption will be valid only in the case that the photon mean
free path ($\ell=1/\chi_{\rm cont}$)444For instance, at the bottom of the
photosphere, $\ell\simeq 100$ km. is small compared to the scale of variation
of the physical quantities, i.e., when the atomic populations depend only upon
the values of the local physical quantities. Besides, it can be shown that if
Kirchoff’s law is further assumed, (e.g., Landi Degl’Innocenti and Landolfi,
2004) the source function vector reduces to

$\mbox{\boldmath$S$}=(B_{\nu}(T),0,0,0)^{\scriptscriptstyle{\rm T}},$ (11)

where, $B_{\nu}(T)$ is the Planck function at the local temperature. These are
the conditions of the so-called _local thermodynamic equilibrium_
approximation (LTE) and have automatically decoupled the RTE from the material
equations. Then, if LTE can be supposed for a given spectral line, the
synthesis of its Stokes profiles simplifies significantly because iterative
procedures are no longer needed. This is graphically explained in Figure 3.

Figure 3: Block diagram of the Stokes profile synthesis under LTE conditions.

In some circumstances, it may be useful to relax the fulfillment of the
Boltzman law and, instead, admit that $\rho_{\alpha j}$ deviate from the LTE
values, $\hat{\rho}_{\alpha j}$, so that

$\beta_{j}=\frac{\rho_{\alpha j}}{\hat{\rho}_{\alpha j}}$ (12)

are _departure_ coefficients that measure how far the conditions are from LTE.
Thus, although radiative transfer remains with the LTE scheme sketched in Fig.
3, the second block is affected by Eq. (12) and the $\beta$’s are needed to
calculate the level populations. As we are going to see, this departure-
coefficient approximation can be very useful for formulating NLTE inversion
procedures (see Section 7.2.3).

Figure 4: Stokes $I$, LTE source function for various atmospheric models: the
umbral model E by Maltby et al. (1986, black line), the penumbral model by Del
Toro Iniesta et al. (1994, red line), the plage model by Solanki (1986, blue
line), and the quiet-Sun models by Gingerich et al. (1971, purple line) and
Vernazza et al. (1981, green line).

###  2.3 The Milne–Eddington approximation

An even more simplified approximation is obtained by further assuming that
thermodynamics is sufficiently described with a source function that depends
linearly on the continuum optical depth,

$\mbox{\boldmath$S$}=(S_{0}+S_{1}\tau_{\rm c})\,\mbox{\boldmath$e$}_{0},$ (13)

where $\mbox{\boldmath$e$}_{0}\equiv(1,0,0,0)^{\scriptscriptstyle{\rm T}}$,
and that the other physical quantities ($B$, $v_{\rm LOS}$, etc.) in the model
are constant throughout the atmosphere, hence defining a constant K. Figure 4
shows the LTE source function (the first component of the vector in Eq. 11) at
525 nm for several realistic model atmospheres, namely, the umbral model E by
Maltby et al. (1986, black line), the penumbral model by Del Toro Iniesta et
al. (1994, red line), the plage model by Solanki (1986, blue line), and the
quiet-Sun models by Gingerich et al. (1971, yellow line) and Vernazza et al.
(1981, green line). The hypothesis of linearity does not seem very accurate
for all the models. Nevertheless, in spite of its seemingly unrealistic
nature, when we are dealing with a weak spectral line, the optical depth
interval at which the line is sensitive to the atmospheric quantities is
usually small enough to consider that a linear source function is not a bad
approximation. There is wide experience in showing how useful the ME
approximation is for inferring average values of the magnetic field vector and
the LOS velocity, starting with the paper by Skumanich and Lites (1987, for a
check with other approaches see , ). The key point is that the RTE has an
analytic solution (Stokes $I$ at $\tau_{\rm c}=0$) under these assumptions
(e.g., Del Toro Iniesta, 2003b):

$\mbox{\boldmath$I$}(0)=(S_{0}+\mbox{\bf
K}^{-1}S_{1})\,\mbox{\boldmath$e$}_{0}.$ (14)

The analytic character of the solution helps in grasping many of the relevant
features in line formation; it cannot reproduce Stokes line asymmetries,555By
Stokes line asymmetries or Stokes profile asymmetries we mean deviations from
the even (Stokes $I$, $Q$, and $U$) or odd (Stokes $V$) functional shape about
the central wavelength of the line. This is commented on in several places in
this review, e.g., Secs. 2.5, 3.1.1 and 3.2, and discussed in Section 5.
though (Auer and Heasley, 1978). Using this useful feature, Landi
Degl’Innocenti and Landi Degl’Innocenti (1985) had the clever idea of
tailoring the functional shape of the source function so that it might be used
to synthesize chromospheric line profiles while preserving an analytic
solution because of the constancy with depth of the propagation matrix. Atomic
polarization is neglected in this modeling. The so-called “field-free
approximation” is assumed. The latter grants substitution of the scalar
components of the source function for those corresponding to the same atom in
the absence of a magnetic field (Rees, 1969). Later on, Lites et al. (1988)
elaborated Landi Degl’Innocenti and Landi Degl’Innocenti’s idea and proposed a
new source function that was incorporated into their inversion code to
interpret the observed profiles of the Mg i b lines at 517.27 and 518.36 nm.
Specifically, they wrote the RTE in terms of the line center optical depth,
$\tau_{0}$, which remains the same as in Eq. (6) but substituting K by
$\mbox{\bf K}^{\prime}\equiv r_{0}\mbox{\bf K}$, where $r_{0}$ is the
continuum-to-line absorption coefficient ratio and with a new source function
$\mbox{\boldmath$S$}^{\prime}$ that follows from two distinct continuum and
line source functions given by

$\begin{array}[]{c}\mbox{\boldmath$S$}_{\rm cont}=\mbox{\boldmath$S$},\\\
{\displaystyle\mbox{\boldmath$S$}_{\rm
lin}=\mbox{\boldmath$S$}-\sum_{i=1}^{2}A_{i}{\rm
e}^{-\varepsilon_{i}\tau_{0}}},\end{array}$ (15)

where $S$ is defined in Eq. (13). The exponential shape of the last two terms
in $\mbox{\boldmath$S$}_{\rm lin}$ tries to mimic the consequences in the
source function of the actual chromospheric rise of temperature. The $A$’s and
$\varepsilon$’s are free parameters that can be tuned to fit the observed
profiles. With this formulation, the analytic solution of the transfer
equation (at $\tau_{0}=0$) turns out to be

$\mbox{\boldmath$I$}(0)=\left[S_{0}+\mbox{\bf
K}^{\prime-1}S_{1}-\sum_{i=1}^{2}A_{i}(\mbox{\bf
K}^{\prime}+\varepsilon_{i}{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}})(\mbox{\bf
K}^{\prime}-r_{0}{\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm 1\mskip-5.0mul}})\right]\mbox{\boldmath$e$}_{0},$ (16)

where ${\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}$ stands for the identity $4\times 4$
matrix.666Note that this is not a non-LTE inversion technique but a
phenomenological approach that can help in fitting the profiles of
chromospheric lines that are indeed formed under conditions far from local
thermodynamic equilibrium.

Further exploiting the analytic character of the Milne–Eddington solution,
slight modifications in the assumptions were also suggested by Landolfi and
Landi Degl’Innocenti (1996) to deal with small velocity gradients and even
with discontinuities along the LOS. In summary, we can say that approximations
to the RTE predicated on keeping the K matrix constant or almost constant are
useful and still a field for exploitation in observational work.

###  2.4 The weak-field approximation

A further simplification of radiative transfer is sometimes used. When the
magnetic field can be assumed constant with depth and weak enough, the
resulting Stokes $V$ profile of many lines turns out to be proportional to the
longitudinal component of the field, regardless of the remaining physical
quantities (see Section 3.1.2). Under this assumption (and for not extremely
weak fields since linear polarization is zero to first order approximation),
the ratio between Stokes $U$ and $Q$ is proportional to the tangent of twice
the field azimuth. The weakness of the field is guaranteed provided that
(e.g., Landi Degl’Innocenti and Landolfi, 2004)

$g_{\rm eff}\frac{\Delta\lambda_{\rm B}}{\Delta\lambda_{\rm D}}\ll 1,$ (17)

where $g_{\rm eff}$ is the effective Landé factor of the line,
$\Delta\lambda_{\rm B}$ is the Zeeman splitting, and $\Delta\lambda_{\rm D}$
is the Doppler width of the line. The effective Landé factor is given by

$g_{\rm
eff}=\frac{1}{2}(g_{u}+g_{l})+\frac{1}{4}(g_{u}-g_{l})[j_{u}(j_{u}+1)-j_{l}(j_{l}+1)],$
(18)

where $g_{u}$ and $g_{l}$ are the Landé factors of the upper and lower level
of the transition, respectively. In $LS$ coupling, those factors are functions
of the quantum numbers:

$g=\frac{3}{2}+\frac{s(s+1)-l(l+1)}{2j(j+1)}.$ (19)

The Zeeman splitting is given by

$\Delta\lambda_{\rm B}=\frac{\lambda_{0}^{2}e_{0}B}{4\pi mc^{2}},$ (20)

where $\lambda_{0}$ is the central, rest wavelength of the line, $e_{0}$ and
$m$ are the charge and mass of the electron, $B$ is the magnetic field
strength, and $c$ stands for the speed of light. For its part, the Doppler
width is given by

$\Delta\lambda_{\rm D}=\frac{\lambda_{0}}{c}\sqrt{\frac{2kT}{m_{a}}+\xi_{\rm
mic}^{2}},$ (21)

where $T$ is the temperature, $k$ is the Boltzmann constant, and $m_{a}$ is
the mass of the atom.

From a formal point of view, Eq. (17) is a good conditioning inequality.
However, in practical terms, one should establish what is meant by _much less
than_ 1\. This is addressed in Section 3.1.2 but we can be sure that the wider
the line, the more the weak-field approximation applies. Hence, broad
chromospheric lines are good candidates for using it. One of the first
attempts at measuring a magnetic field with a chromospheric line, known to the
authors of this review, was carried out as early as 1990 by Martínez Pillet et
al. who (photographically) observed Stokes $I$ and $V$ profiles of the Ca ii H
line and interpreted them in terms of the weak-field approximation. This
approach remains useful as interest in the chromosphere increases (e.g., de la
Cruz Rodríguez et al., 2013).

###  2.5 The MISMA hypothesis

Driven by the ubiquitous appearance of Stokes profile asymmetries in
observations, Landi Degl’Innocenti (1994) suggested considering the
atmospheric physical quantities, instead of deterministic stratifications, to
have stochastic distributions about mean values _with possible correlation
effects_ among them. Assuming that the source function nevertheless varies
linearly with depth through the whole atmosphere and that the propagation
matrix stays constant at the spatial scale of each of the realizations of such
a common stochastic distribution, he found an analytic solution for the
transfer equation. Certainly inspired by the Landi Degl’Innocenti’s proposal,
Sánchez Almeida et al. (1996) put forward a new approach. Realizing that the
wavelength symmetries in the propagation matrix elements do indeed avoid such
Stokes profile asymmetries in the absence of LOS velocity gradients in the
regular formulation of the transfer problem (Landi Degl’Innocenti, 1992), they
proposed that the solar atmosphere may be pervaded by MIcro-Structured
Magnetic Atmospheres (MISMAs). The hypothesis implies a highly inhomogeneous
atmosphere at scales much smaller than the photon mean free path whereby the
integration of Eq. (6) turns out to be very difficult. An alternative
formulation is thus in order by locally averaging the propagation matrix and
the emission vector. The resulting equation reads

$\frac{{\rm d}\mbox{\boldmath$I$}}{{\rm d}s}=-\left<\mbox{\bf
K}^{\prime}\right>(\mbox{\boldmath$I$}-\mbox{\boldmath$S$}^{\prime}).$ (22)

It formally looks very much like the regular RTE but is formulated in terms of
geometrical distances, $s$; $\mbox{\bf K}^{\prime}=\chi_{{\rm cont}}\mbox{\bf
K}$;

$\mbox{\boldmath$S$}^{\prime}\equiv\left<\mbox{\bf
K}^{\prime}\right>^{-1}\left<\mbox{\bf K}^{\prime}\mbox{\boldmath$S$}\right>;$
(23)

and the averages are taken over a distance $\Delta s$ that may vary along the
optical path. The distance $\Delta s$ is supposed to be still smaller than
$\ell$ for Stokes $I$ to be assumed constant within its range. In addition,
the averages are considered to vary smoothly along the line of sight.

With all these assumptions, Eq. (22) is formally the same as Eq. (6). All the
mathematical tools developed to solve the latter can be used to find a
solution to the former. This is so despite the (numerically) inconvenient
formulation in terms of geometrical distances: it requires either non-equally-
spaced grid points or an increase in computation time. The good news is that,
since correlations may exist among the physical parameters of the
microstructures, the symmetry properties of matrix $\left<\mbox{\bf
K}^{\prime}\right>$ are automatically destroyed. Hence, asymmetric Stokes
profiles can appear naturally.

##  3 Degrees of approximation in the model atmospheres

Provided that physical atmospheric quantities are bounded functions of the
optical depth, we can safely expect that they are either continuous or have
some jump (Heaviside-like) discontinuities throughout the line formation
region. Therefore, except for the discontinuity points, a Taylor expansion
approximation seems simple and sensible. The good feature of Taylor expansions
is that you can keep them at a given order of approximation that can be
subsequently increased if needed. The sequential approach is of great help in
following the principle of Occam’s razor — _lex parsimoniae_ — which, in our
opinion, should prevail in the interpretational work. The question arises as
to whether an order of approximation is useful or whether it should be
increased to give account of the observations. The answer must be found in the
degree of accuracy with which we are trying to reproduce the observables.
Hence, it has to do with the balance between the signal and the noise: if the
next order of approximation only introduces variations that are below, say,
three times the rms noise, $\sigma$, then its use is discouraged. If, on the
contrary, the difference between the observed and synthetic profiles is
greater than $3\sigma$, its use may be advisable.777By adopting $\sigma$ as a
measure of noise we are assuming that the noise statistics is Gaussian and
this seems a common and sensible assumption as well. Requiring signals to be
larger than $3\sigma$, therefore, implies more than 99.7 % certainty in the
detection. We refer the reader to Del Toro Iniesta and Martínez Pillet (2012)
for a discussion on polarimetric accuracy and signal-to-noise ratio. For
Bayesian selection among model atmospheres, see Asensio Ramos et al. (2012).
Let us postpone the discussion to the following sections and present here the
various atmospheres we are considering. We start with the zeroth order
approximation and assume that physical quantities are constant with depth to
continue with gradients, higher order variations, and jumps or
discontinuities.

Figure 5: Examples of ME Stokes profiles of the Fe i line at 617.3 nm as
observed with an instrument whose Gaussian spectral PSF has a FWHM of 6 pm.
Two model atmospheres are used that differ only in the magnetic field
strength: $B=1200$ G for the black lines and 200 G for the red ones.

###  3.1 Constant physical quantities

Let us distinguish among three possibilities, namely, the Milne–Eddington
approximation, the weak-field approximation, and an atmosphere where $B$ and
$v_{{\rm LOS}}$ are constant but where thermodynamics is properly accounted
for with a realistic stratification of temperature.888Indeed, two variables
are needed for specifying the thermodynamical state of the medium. However,
most of the spectral lines used in typical observations present a very low
dependence, if any, on pressure. Therefore, we shall assume that pressure is
stratified according to hydrostatic equilibrium throughout the paper.

####  3.1.1 The Milne–Eddington atmosphere

As commented on in Sect. 2.3, a Milne–Eddington atmosphere provides an
analytic solution to the RTE. With nine parameters, the Stokes profiles of a
spectral line can be synthesized. The model parameters are the three
components of the magnetic field, $B$, $\gamma$, and $\varphi$, the LOS
velocity, $v_{\rm LOS}$, and the so-called thermodynamic parameters: the line-
to-continuum absorption coefficient, $\eta_{0}$ ($=1/r_{0}$), the Doppler
width of the line, $\Delta\lambda_{\rm D}$, the damping parameter, $a$, and
the two coefficients for the source function, $S_{0}$ and $S_{1}$. The actual
values of $\eta_{0}$, $\Delta\lambda_{\rm D}$, and $a$ may vary significantly
throughout the atmosphere. Therefore, assigning one single value for each may
be, say, risky. Experience, however, indicates that this is possible.
Reasonable fits to actual data can be obtained with this approximation and we
can even understand the relationship between the single-valued parameters and
their actual stratification (Westendorp Plaza et al., 1998). Only Stokes
profiles with definite symmetry properties can be formed in an ME atmosphere.
Stokes $I$, $Q$, and $U$ are even functions of wavelength while Stokes $V$ is
odd. This is a consequence of the absence of velocity gradients (Auer and
Heasley, 1978) and will be discussed later in Section 5. Figure 5 shows two
examples of ME profiles corresponding to the Fe i line at 617.3 nm as observed
with an instrument whose (Gaussian) spectral profile (point spread function,
PSF) has a full width at half maximum (FWHM) of 6 pm. The thermodynamic model
parameters are $\eta_{0}=5.06$, $\Delta\lambda_{\rm D}=2.6$ pm, $a=0.22$,
$S_{0}=0.1$, and $S_{1}=0.9$; they come from a fit to the FTS spectrum (Kurucz
et al., 1984; Brault and Neckel, 1987). The magnetic inclination and azimuth
are both equal to 30∘; $B=1200$ G for the black lines and 200 G for the red
ones.

####  3.1.2 The weak-field atmosphere

As stated in Sect. 2.4, when $B$ is constant with depth and very weak, then
the Stokes $V$ profile turns out to be proportional to the longitudinal
component of the magnetic field independently of the remaining quantities. It
can be shown (e.g., Landi Degl’Innocenti and Landolfi, 2004) that

$V(\lambda)\simeq-g_{\rm eff}\,\Delta\lambda_{\rm B}\cos\gamma\,\frac{\partial
I_{\rm nm}}{\partial\lambda},$ (24)

where $I_{\rm nm}$ is the non-magnetic Stokes $I$ profile, corresponding to
the line in the absence of a magnetic field. Equation (24) has been key for
many magnetic inferences. In fact, written as $V=CB_{\parallel}$, it is known
as the magnetographic equation since it provides a calibration of the
magnetographic signal. When magnetographs used only one or two wavelength
samples of the circular polarization, the magnetographic equation was indeed
the only means of obtaining estimates of the component of the magnetic field
along the line of sight. Nowadays, with modern magnetographs providing more
samples in all four Stokes parameters, that equation is still useful for
morphological, qualitative estimates but cannot be trusted everywhere and
under all circumstances. The modern way to evaluate $C$ indeed implies some
radiative transfer calculations in given model atmospheres (e.g., Martínez
Pillet et al., 2011), and these calculations readily show that the
approximation saturates at low magnetic field strengths. In the left panel of
Fig. 6, we plot the maximum of the Stokes $V$ profile as a function of the
field strength (the field is along the LOS, $\gamma=0^{\circ}$) with an
instrumental profile FWHM of 6 pm (asterisks) and of 8.8 pm (diamonds). In
solid lines, the linear (red) and quadratic (blue) fits are also shown. Only
strengths up to 600 G are plotted because the relationship is evidently
nonlinear above that threshold. For weaker fields, it is apparent that the
instrumental broadening of the profiles helps linearity to hold as differences
between the linear and quadratic fits are smaller for the broader PSF. Those
differences are for most of the points above $3\cdot 10^{-3}I_{\rm c}$; that
is, more than $3\sigma$, with $\sigma$ being the noise level of the
polarization continuum signal of typical observations. Such differences are
clearly detectable by current means. Hence, the approximation loses validity
for yet weak fields. Deviations from linearity are even clearer if one sees
the green lines in the figure, which correspond to linear fits including only
data points for which $B$ is less than 200 G. In our example, the weak field
approximation for the Stokes $V$ peaks breaks down at fields stronger than 300
G with a FWHM of 6 pm and stronger than 400 G with a FWHM of 8.8 pm.
Certainly, if the instrument has a narrower spectral PSF or if the noise is
smaller, the approximation fails earlier. The approximation clearly worked
better for older instruments.

Figure 6: Maximum of the Stokes $V$ profile as a function of the magnetic
field strength for a longitudinal field (left panel). Maximum of the Stokes
$Q$ profile as a function of the square magnetic field strength (right panel).
Asterisks correspond to an instrumental profile FWHM of 6 pm and diamonds to a
FWHM of 8.8 pm. Red lines represent linear fits to the points; blue lines
display quadratic fits; green lines correspond to fits for fields weaker than
200 G.

Further arguments can be supplied for the user to be cautious about weak field
assumptions with typical, visible photospheric lines. The first one is that
Eq. (24) is hardly applicable, as shown in Fig. 7, not only because Stokes $V$
does not follow it but because Stokes $I$ deviates from $I_{{\rm nm}}$ even
sooner (and, up to first order, $I=I_{{\rm nm}}$ must hold for Eq. (24) to be
valid; e.g., Landi Degl’Innocenti and Landolfi, 2004). In the left column of
the figure, the differences between the left-hand and the right-hand members
of the equation are plotted. Colors correspond to 600 G (black), 500 G (red),
400 G (blue), 300 G (green), 200 G (purple), and 100 G (dark green). The
dashed, horizontal purple lines mark the $3\sigma$ level. The upper rows are
for a FWHM of 6 pm and the bottom row is for a FWHM of 8.8 pm. The plots in
the left column are of course consistent with the results from Figure 6. Those
in the right column are illustrative of how Stokes $I$ varies with the
magnetic field strength. Differences between the various profiles can easily
be discerned above the $3\sigma$ level. When the profiles themselves are
affected by noise, unlike in these plots, detecting the differences may be
more difficult but the message is clear: contrary to the common belief, the
Stokes $V$ profile is not the only tool for estimating the longitudinal
component of weak magnetic fields; Stokes $I$ helps a lot and should not be
forgotten.

The second argument concerns the diagnostic capability for typical lines to
disentangle $B$ from $\gamma$ in the weak-field regime. Most statements about
the only accurate retrieval to be the longitudinal magnetic field component
are based on Eq. (24), as if it were the only available tool from radiative
transfer. Stokes profiles other than $V$ are often obliterated. It is easy to
understand (e.g., Landi Degl’Innocenti and Landolfi, 2004), however, that the
mere deviations between $I$ and $I_{\rm nm}$ we have seen in Fig. 7 should
imply the appearance of linear polarization signals (provided that the
inclination is different from zero): such Stokes $I$ deviations from $I_{\rm
nm}$ are second order terms in an expansion of all four Stokes profiles.999The
expansion is in terms of powers of a dimensionless parameter that scales the
vector magnetic field, and that is only valid when $\Delta\lambda_{\rm
B}\rightarrow 0$. At second order, Stokes $Q$ and $U$ are no longer zero (or
below the noise) either and start to provide additional information. It can
also be proven (e.g., Landi Degl’Innocenti and Landolfi, 2004) that $Q\propto
B^{2}\sin^{2}\gamma$, as shown in the right panel of Fig. 6, where the maximum
of Stokes $Q$ is plotted against $B^{2}$ for a field that is inclined
$45^{\circ}$ with respect to the vertical.101010Stokes $Q$ is assumed to be
defined here in the reference frame where Stokes $U$ is zero (constant
magnetic azimuth). Here, deviations between linear and quadratic fits are
smaller than for the $V$ case (note that the $Y$ scale is an order of
magnitude smaller) but the interesting point is that, above $B=200$ G, linear
polarization signals begin to be larger than $3\sigma$ and, hence, detectable.

Figure 7: Differences between the Stokes $V$ profile and its weak-field
approximation (left column) and differences between the Stokes $I$ profile and
that for a zero field strength. Colors indicate values of the longitudinal
component of the field. The dashed horizontal lines mark the $3\sigma$ level
of typical, modern observations. The upper row is for a FWHM of 6 pm and the
bottom one for a FWHM of 8.8 pm. Colors correspond to 600 G (black), 500 G
(red), 400 G (blue), 300 G (green), 200 G (purple), and 100 G (dark green).

A third argument we want to bring to the reader’s attention is related to the
common belief that weak fields are hardly distinguished from strong fields
(say above 1 kG) with a filling factor significantly smaller than 1. We will
return to this issue in Secs. 6.2 and 8.2, as the problem has already been
discussed in the literature (e.g., Del Toro Iniesta et al., 2010). Let us only
mention here that the loss of linearity of Stokes $V$ above, say, 400 G and,
most importantly, the behavior of Stokes $I$ are reasons enough for the two
types of atmospheres to be distinguished by observational means.

If Eq. (24) were universally accepted, then it would indicate that the RTE is
almost useless since the emergent profile is proportional to one of the matrix
elements of K. Elementary mathematics readily explain that this is not
possible except for, perhaps, a small value range of fields. In summary, we
must acknowledge that Stokes $V$ is not proportional to the longitudinal
component of the magnetic field.

####  3.1.3 Constant vector magnetic field and LOS velocity

There is still a third option to deal with constant $B$ and $v_{\rm LOS}$.
Imagine that the atmosphere is a regular one as far as thermodynamics is
concerned but where the magnetic and dynamic quantities do not vary with
depth. Since the propagation matrix is no longer constant, no analytic
solution of the RTE is available.111111The clause is not very rigorous but is
true in practical terms. Indeed, one can conceive other K stratifications that
still allow an analytic solution of the RTE (Landi Degl’Innocenti and Landi
Degl’Innocenti, 1985). One is then led to use numerical techniques to
synthesize the spectrum. The atmosphere, however, is greatly simplified since
the number of parameters is reduced. This can be very helpful for quicker
analyses of the data or as a makeshift for more elaborate subsequent
approaches that include variations of $B$ and $v_{\rm LOS}$ with the optical
depth. This is the approximation used for the first version of the SPINOR code
(Solanki et al., 1992a) or as an option in the SIR code (Ruiz Cobo and Del
Toro Iniesta, 1992).

###  3.2 Physical quantities varying with depth

The community has gathered a great deal of evidence about variations of $B$
and $v_{\rm LOS}$ along the optical path everywhere over the solar disk. In
addition, physical laws such as those of magnetic flux and mass conservations
demand that these quantities vary with optical depth in a number of
structures. The approximations in the former subsections cannot then be
considered but as first-step approaches or simplified descriptions of reality.
In any case, we can safely assume that stratifications of the physical
quantities are bounded functions of $\tau_{\rm c}$ (or whichever variable
parameterizing the optical path), as we admitted in the beginning of this
Section.

A historical landmark for the full acknowledgement of LOS velocity gradients
from an observational point of view was established by the discovery by Mickey
and Orrall (1974) and Illing et al. (1974a, b, 1975) of a broadband circular
polarization in sunspots. The true explanation was already suggested in the
last of those papers, although schematically founded on the assumption of two
slabs with _different velocity_ and magnetic field strengths. The broadband
observations were soon related to spectral line net circular polarization (the
integral of the Stokes $V$ profile over the wavelength span of the line):
Grigorjev and Katz (1975) computed all four Stokes profiles in the presence of
an LOS velocity gradient and certainly obtained asymmetric profiles; later on,
Auer and Heasley (1978) demonstrated that a necessary and sufficient condition
for such a net circular polarization had to be found in velocity gradients
along the line of sight, although they were neglecting magneto–optical
effects. Rigorous derivations (including dispersion effects) have later been
obtained and can be found, for example, in the elegant work by Landi
Degl’Innocenti and Landi Degl’Innocenti (1981). The symmetry properties of the
propagation matrix elements predict no net circular polarization (or Stokes
$V$ area asymmetry) in the absence of an LOS velocity gradient. Other
mechanisms such as insufficient spatial resolution that implies mixtures of
individual atmospheres within a pixel, may produce asymmetries in the peaks
(the so-called amplitude asymmetries) but the integral of $V$ will remain
zero. Therefore, any net circular polarization is unambiguous observational
evidence for the presence of velocity gradients. And Stokes $V$ area
asymmetries are observed practically everywhere. Unfortunately, no such
unambiguous evidence exists for the presence of magnetic field gradients,
although we know on physical grounds there are plenty of them, such as those
through magnetic canopies where a magnetic layer is overlaying a non-magnetic
one.

####  3.2.1 Parameterizing the stratifications

Among the numerical codes relevant to this review (see Sect. 7) there are some
that acknowledge variations of $B$ and $v_{\rm LOS}$. We deal here with what
might be called “normal” or “regular” stratifications, such as those employed
by Ruiz Cobo and Del Toro Iniesta (1992), Frutiger and Solanki (1998), Socas-
Navarro et al. (2000), and Socas-Navarro (2001), and leave some others,
devoted to specific solar features, to the following paragraphs.

Since the number of depth grid points used for the numerical integration of
the RTE can be high, it may be advisable to reduce the degrees of freedom of
the variations with depth of the physical quantities. As commented on above, a
reasonable approach would be to follow higher order polynomials in a stepwise
form. From constant values to linear, parabolic, third-order polynomial
dependences, and so on. Then, if we assume, for instance, that $v_{\rm LOS}$
is linear with $\tau_{\rm c}$, we only need to specify the velocity at two
grid points (_nodes_ in SIR’s terminology) and three if it is parabolic, hence
reducing the number of free parameters of the model. We do not need to specify
$T$, $B$, and $v_{\rm LOS}$ at every single point we use for solving the RTE
but only at a few of them. We shall see in Sect. 7 that one can go even
further with this kind of approach and consider more involved optical depth
dependences if necessary.

####  3.2.2 The MISMA atmosphere

As we explained in Sect. 2.5, the MISMA hypothesis guarantees the appearance
of Stokes profile asymmetries but at the expense of introducing a significant
number of extra free parameters. In fact, even in the simplest MISMA
atmosphere (Sánchez Almeida and Landi Degl’Innocenti, 1996), where all the
micro-structures are described by ME atmospheres, one has in principle as many
as ten free parameters per needed component (also known as micro-structure).
To the nine regular ME parameters, the volume occupation fraction for each
micro-structure must be added. In more complicated MISMAs, the number of
parameters is even higher (Sánchez Almeida, 1997). Moreover, in spite of the
very detailed physical description where equilibrium equations are required
for slender flux tubes, the inclination and azimuth of the magnetic field are
kept constant throughout the whole atmosphere, which does not seem very
realistic. (Canopies are found almost everywhere owing to the fanning out of
magnetic field lines with height.) Last, but not least, when the structuring
of the atmosphere is established at sizes comparable to $\ell$, the average
propagation matrix does not result in an RTE as in Eq. (22), which is no
longer valid. Modern observations with continuously increasing spatial
resolution do indeed show this kind of structuring both in quiet and active
regions and sunspots. For example, single magnetic flux tubes of approximately
150 km size have been fully resolved by Lagg et al. (2010); their evolution
followed for half an hour by Requerey et al. (2014); and the internal
structure of network magnetic structures revealed (Martínez González et al.,
2012) with Sunrise/IMaX observations (Martínez Pillet et al., 2011; Barthol et
al., 2011). In our opinion, the MISMA hypothesis, being a clever idea for
producing asymmetries, is advisable as a “when-all-else-fails” atmosphere but
there are yet conventional radiative transfer treatments that provide
reasonable interpretation of the observations.

####  3.2.3 Other special atmospheres

This subsection is devoted to three special cases where the physical scenario
envisaged to explain the observations requires a specific configuration that
is not intended to be universally valid. Those specific configurations,
however, help in interpreting the Stokes profiles emerging from given solar
features.

##### Interlaced atmospheres

Imagine that you can assume that your line of sight is piercing a number $n$
of alternate boundaries
$\left\\{s_{i}\right\\}_{i=1,\ldots,n}(s_{1}<s_{2}<\ldots<s_{n})$ between two
distinct atmospheres, as when observing from a side two identical thin flux
tubes that are close but not stuck to each other. In such a scenario, the
structuring of the atmosphere is comparable in size with $\ell$ and,
therefore, the MISMA hypothesis does not hold. If you happen to know the
solution, $\mbox{\boldmath$I$}_{\pm 1}$, of the RTE in each of the two
atmospheres, labeled $\pm 1$, Del Toro Iniesta et al. (1995) found out that
the formal solution to the problem is

$\mbox{\boldmath$I$}(s)=\mbox{\boldmath$I$}_{+1}(s)+\sum_{i=1}^{n}(-1)^{n-i}\left[\prod_{j=i}^{n}\mbox{\bf
O}_{(-1)^{n-j}}(s_{j+1},s_{j})\right]\Delta\mbox{\boldmath$I$}(s_{i}),$ (25)

for any $s\in[s_{n},s_{\rm lim}]$, where the $+1$ atmosphere is assumed to be
the outermost one,
$\Delta\mbox{\boldmath$I$}(s_{i})\equiv\mbox{\boldmath$I$}_{-1}(s_{i})-\mbox{\boldmath$I$}_{+1}(s_{i})$,
and $\mbox{\bf O}_{\pm 1}$ are the evolution operators for both atmospheres
(e.g. Landi Degl’Innocenti and Landi Degl’Innocenti, 1985). Equation (25) is
at the root of the flux-tube inversion code by Bellot Rubio et al. (1997, see
,  as well). A different treatment of discontinuities along the line of sight
was proposed by Borrero et al. (2003) where the density of depth grid points
is increased in the discontinuity neighborhood.

##### Atmospheres with Gaussian profiles

The existence of net circular polarization in the penumbrae of sunspots was
also the driver for Bellot Rubio (2003) to propose an implementation of the
uncombed model by Solanki and Montavon (1993). The scenario is based on two
components; namely, a magnetic component and a penumbral magnetic flux tube,
the latter occupying a fractional area of the resolution element. The model
parameters of the penumbral tube are built by Gaussian modifications (in
depth) of those in the background. All the Gaussians have the same width and
are located at the same depth, but their amplitudes depend on the specific
model parameter, of course. With these premises, the SIR code was modified
into the so-called SIRGAUS code, which has been used, among others by Jurčák
et al. (2007), Jurčák and Bellot Rubio (2008), Ishikawa et al. (2010), and
Quintero Noda et al. (2014).

##### Atmospheres with jump discontinuities

Discontinuities can be treated numerically by decreasing the depth grid step
or by using Eq. (25). A specific implementation of such discontinuities was
first used by Louis et al. (2009) for an analysis of sunspot light bridges.
Like SIRGAUS, it is based on a modification of the SIR code to take this
particular scenario into account. In it, two magnetic atmospheres coexist in
the resolution pixel: a background atmosphere whereby $B$ and $v_{\rm LOS}$
are constant with depth, and another magnetic atmosphere where those
quantities have a Heaviside-like discontinuity. This code (called SIRJUMP) has
also been used in practice, e.g., by Martínez González et al. (2012) and Sainz
Dalda et al. (2012).

##  4 Degrees of approximation in the Stokes profiles

Since the ultimate goal of inversions is the _bona fide_ reproduction of
observed profiles, an analysis of the properties of Stokes spectra as
functions of the wavelength is in order. Such an analysis should be aimed at
finding the most conspicuous characteristics of the profiles in order for
these characteristics to be the best reproduced among all the features. In
other words, if, for instance, a given Stokes $I(\lambda)$ profile shows only
small deviations from a Gaussian, we should aim to obtain the Gaussian that
best simulates the profile and identify the model parameters responsible for
this bulk behavior. In some cases we may be satisfied just with this “coarse”,
or not very detailed, description and leave small deviations or nuances to
further, in-depth analysis that might even be carried out separately. As we
are going to see, this approximation of incremental complexity for the
profiles is well in line with the successive approximations we have described
for the model atmospheres in Section 3.

The Stokes $Q$, $U$, and $V$ profiles and Stokes $I$ in line depression; that
is,

$I_{\rm d}\equiv 1-\frac{I}{I_{\rm c}},$ (26)

as functions of $x\equiv\lambda-\lambda_{0}$ (where $\lambda_{0}$ is the
central wavelength of the line), can be decomposed as sums of even and odd
functions of $x$, as any other function defined over ${\mathchoice{\rm
I\mskip-4.0muR}{\rm I\mskip-4.0muR}{\rm I\mskip-4.5muR}{\rm
IRmskip-5muR}}$.121212Strictly speaking, we should take a mean LOS velocity
wavelength shift into account as well, if it is large enough. Specifically, if
we call $S(x)$ any one of the profiles, then

$S(x)=S_{+}(x)+S_{-}(x),$ (27)

where

$S_{+}(x)\equiv\frac{S(x)+S(-x)}{2}\,\,\,{\rm
and}\,\,\,S_{-}(x)\equiv\frac{S(x)-S(-x)}{2}.$ (28)

By construction, $S_{+}$ is even and $S_{-}$ is odd.131313The property is also
valid for Stokes $I$. We have, however, chosen $I_{\rm d}$ for reasons that
will become clear a little later in this Section.

Figure 8: Differences between the Stokes profiles of the Fe i line at 617.3 nm
as synthesized in two model atmospheres that differ in the LOS velocity. See
text for details.

This parity property is very interesting because, as we have seen in former
Sections, the Stokes profiles of any line formed in the absence of velocity
gradients have definite symmetry (parity) properties. Since asymmetries in
regular profiles are relatively small, that is, the profiles usually display a
predominant parity character (even for Stokes $I$, $Q$, and $U$, and odd for
Stokes $V$), a sum of even and odd profiles may give account of the observed
spectra as if the opposite parity component was indeed a _perturbation_
related to velocity gradients. This can explain the success of ME inversion
codes for fitting many observations (cf. Westendorp Plaza et al., 1998; Orozco
Suárez et al., 2010). The ME atmosphere accounts for the main bulk of the
observed Stokes profiles. In Fig. 8 we plot the differences among the Stokes
profiles of the Fe i line at 617.3 nm as synthesized in two model atmospheres.
Both have the HSRA (Gingerich et al., 1971) stratification of temperature with
$B=1500$ G and $\gamma=\varphi=30$∘. One of the models has a constant $v_{\rm
LOS}=1.87$ km s-1and the other a small gradient from $v_{\rm LOS}=2$ km s-1at
the bottom of the atmosphere through 1.75  km s-1at the top. Both have
$\xi_{\rm mac}=1$ km s-1and have been convolved with the IMaX PSF. Note that
these differential profiles display almost the opposite parity character to
their corresponding Stokes profiles.

Figure 9: First six eigenprofiles for Stokes $I$, $Q$, $U$, and $V$. They are
obtained from observations of the Fe i line pair at 630.1 and 603.2 nm.
Adapted from Rees et al. (2000).

A description with only $S_{+}$ or $S_{-}$ can thus provide a first approach
analysis to a large number of observations. Indeed, $S_{+}$ should be good for
$I$, $Q$, and $U$, and $S_{-}$ for $V$ as the differences are smaller than our
“nominal” noise of $10^{-3}\,I_{\rm c}$. This, of course, cannot always be the
case. Very peculiar Stokes profiles are often observed as our polarization
accuracy increases. For instance, Sigwarth et al. (1999) first reported the
observation of one-lobed $V$ profiles that were later studied in detail by
Grossmann-Doerth et al. (2000) and Sigwarth (2001). Most of these profiles are
found in the internetwork (e.g., Sainz Dalda et al., 2012).

A different description of the Stokes profiles as functions of wavelength was
proposed by Rees et al. (2000) who suggested that they can be described as
sums of given _principal components_ or _eigenprofiles_. If those
eigenprofiles are contained in a database and are properly selected, they can
increasingly give account of the profile shapes just by increasing the number
of principal components in the expansion. An example of such eigenprofiles is
given in Fig. 9. By adding these components properly weighted, the
corresponding Stokes profiles are synthesized. This is the basis for all the
PCA inversion techniques presented so far and the concept is fairly simple.

A similar approach to that of PCA was proposed by Del Toro Iniesta and López
Ariste (2003), based on the fact that Stokes $I_{\rm d}$, $Q$, $U$, and $V$
belong to ${\mathchoice{\rm I\mskip-4.0muL}{\rm I\mskip-4.0muL}{\rm
I\mskip-4.5muL}{\rm IRmskip-5muL}}^{2}$, the space of square integrable
functions over ${\mathchoice{\rm I\mskip-4.0muR}{\rm I\mskip-4.0muR}{\rm
I\mskip-4.5muR}{\rm IRmskip-5muR}}$. Since ${\mathchoice{\rm
I\mskip-4.0muL}{\rm I\mskip-4.0muL}{\rm I\mskip-4.5muL}{\rm
IRmskip-5muL}}^{2}$ is a Hilbert space with a well defined scalar product, an
exact, infinite expansion of the profiles is possible in terms of any of the
several bases of the space. Among those basis systems, Del Toro Iniesta and
López Ariste selected the family of Hermite functions, $h_{n}(x)$, because of
the similarity between the shapes of the first few elements of the family and
the Stokes profiles (see Figure 10). Somehow, the Hermite functions (see the
aforementioned paper for a definition) provide a suitable basis for
approximating the observed profiles with finite expansions of a few terms.
Apart from their possible use in inversion codes that has not been
investigated so far, the expansion of Stokes profiles in terms of Hermite
functions has been used by Toussaint et al. (2012) for compressing observed
Stokes profiles and by Harker (2012) for an automatic solar active region
detection. The first authors, after expansion of the Stokes profiles, only
keep the coefficients compressed with a conventional algorithm. This way they
reduce the storage space by a factor 20 while keeping most of the information
virtually noise free. The latter author discriminates the different active
regions after looking at the complexity of the emerging Stokes profiles as
described by their Hermite-function expansion coefficients.

Figure 10: First six Hermite functions. The abscissa has to be understood as a
normalized wavelength (by the Doppler width of the line, for instance).

##  5 A synthesis approach

As described in the introduction, an approximate knowledge as to how the
Stokes profiles react to the various model parameters is advisable as it helps
to select the adequate observables as “orthogonal” as possible. The ideal way
to explore the diagnostic capabilities of Stokes profiles is by means of
response functions (see Sect. 6). Tackling the problem head-on, that is,
synthesizing the profiles in different model atmospheres, may help grasp basic
ideas on the Stokes profile behavior, though. The idea is to study how the
Stokes profiles vary when the model parameters are modified. This is the aim
of this section.

###  5.1 Constant atmospheres

As we have been doing in the two previous sections, let us start with the
easiest case of atmospheres that do not vary with optical depth and,
specifically, with ME atmospheres, since their analytic solution of the RTE
enables a quick numerical overview of the space of model parameters. Figure 11
shows Stokes $I$, $Q$, and $V$ for the Fe i line at 617.3 nm with the same
thermodynamic parameters used in Figs. 5 and 7. Since the linear polarization
$L^{2}\equiv Q^{2}+U^{2}$ is rotationally invariant and $\varphi$ is constant
throughout the ME atmosphere, we are assuming to have selected the preferred
reference frame where $U$ is identically zero, so that $L=Q$. In practical
terms we have selected $\varphi=0^{\circ}$ for all the profiles. The magnetic
field strength is $B=300$, 500, 900, and 1200 G for the four rows from top to
bottom. The magnetic inclination is encoded in color: $\gamma=0$ (dark green),
15 (purple), 30 (pink), 45 (green), 60 (blue), 75 (red), and $90^{\circ}$
(black). If we again assume a typical noise $\sigma=10^{-3}I_{\rm c}$, then
the small differences in the core of Stokes $I$ for $B=300$ may not be
detected and a neat distinction between $\gamma=0$ and $15^{\circ}$ or
$\gamma=90$ and $75^{\circ}$ may hardly be reachable with Stokes $Q$ and $V$
at a $3\sigma$ level. Most certainly, however, there should not be any problem
to distinguish between 0, 30, 60, and $90^{\circ}$. Of course, the dependences
of $I$ and $V$ on $\gamma$ are significant enough when the field is stronger.
One should not be restricted to longitudinal components even for this small
field strength, and the situation improves further for lower noise levels.
This is a well-known issue: in polarimetric observations we are photon
starved, indeed as much as night-time astronomers may be for detecting very
faint objects.

Figure 11: Stokes $I$, $Q$, and $V$ as functions of the magnetic field
inclination: dark green is for $\gamma=0$, purple for 15, pink for 30, green
for 45, blue for 60, red for 75, and black for $90^{\circ}$. The magnetic
field strength is different for the various rows: from top to bottom, $B=300$,
500, 900, and 1200 G. The magnetic azimuth is identically zero.

Figure 12: Stokes $I$, $Q$, and $V$ as functions of the magnetic field
strength: black is for $B=1200$, red for 900, blue for 500, green for 300, and
purple for 200G. The magnetic field inclination is different for the various
rows: from top to bottom, $\gamma=15$, 30, 60, and $75^{\circ}$. The magnetic
azimuth is identically zero.

The alternative way to gauge the sensitivity of the Stokes profiles to $B$ and
$\gamma$ by synthesizing the profiles is shown in Figure 12. The rows
correspond to different inclinations: $\gamma=15$, 30, 60, and $75^{\circ}$
from top to bottom. The magnetic field strength is this time encoded in
colors: $B=1200$ (black), 900 (red), 500 (blue), 300 (green), and 200 G
(yellow). This complementary view shows that the dependence on $B$ is in fact
stronger than on $\gamma$. Therefore, properly sampled Stokes profiles with
enough polarimetric accuracy should be able to provide the required
information to infer the magnetic field strength and inclination separately
for most of the strength spectrum. Weaker fields will have bigger
uncertainties for sure, but they should not imply a theoretical inability. As
a matter of fact, the weaker the fields we want to explore, the smaller the
noise we need in our observations, but this is somehow obvious.

###  5.2 Depth-dependent atmospheres

In the solar atmospheres, physical quantities do vary with depth.
Acknowledging such variations almost always implies resorting to numerical
solutions of the transfer equation. That was first done by Beckers (1969a, b).
Numerical results by Staude (1969, 1970) and by Wittmann (1971, see Fig. 13)
soon appeared and numerical codes were described (e.g. Wittmann, 1974, Landi
Degl’Innocenti, 1976). Those first numerical codes capable of synthesizing the
Stokes profiles were based on the fourth-order Runge–Kutta algorithm that is
very accurate at the price of being very computationally expensive. A
generalization of the method by Feautrier (1964) to polarized light was
proposed by Auer et al. (1977) and later modified by Rees and Murphy (1987) in
order to take magneto-optical effects into account. A fast solution of the
RTE, after being reformulated as an integral Volterra equation (first
suggested by Staude, 1969), was formulated by Rees et al. (1989) with their
so-called DELO method.

Figure 13: Comparison between the observed and computed Stokes $I$ profile of
the Fe i line at 630.25 nm. Adapted from Wittmann (1971).

An improvement in accuracy and computational speed was obtained by Bellot
Rubio et al. (1998) with an Hermitian method based on developing the Stokes
vector as a fourth-order polynomial with depth. In Fig. 14 we show an example
of a synthesis of the same line as in Fig. 13, where clear asymmetries in
wavelength in the Stokes profiles can be seen. According to Sections 3.2 and
4, such asymmetries have been produced by the variation with depth of physical
quantities.

Learning how the various profile features of the four Stokes parameters depend
on the many model parameters is certainly difficult and cannot be summarized
in this paper. Experience, however, can train a researcher to be able to
deduce —many times after a quick glance (the _art_)— a specific stronger
$v_{\rm LOS}$ or $B$ here or there in the atmosphere. The situation is
therefore much more complicated than for the ME case and one should rely upon
inversions.

Figure 14: Stokes profiles of Fe i line at 630.25 nm as synthesized with an
Hermitian method. Adapted from Bellot Rubio et al. (1998).

###  5.3 MHD simulations

The advent of magnetohydrodynamic (MHD) simulations such as those by Vögler
(2003), Vögler et al. (2005), and Rempel (2012) has opened a new window on the
exploration of the solar photosphere. They have enabled calculations that may
help to envision what is expected from observations and to interpret them. The
simulations also provide predictions that can be confronted with them. The
enrichment has been remarkable because the realistic atmospheres resulting
from the simulations have physical quantities varying along the optical path
without any _a priori_ assumptions and may be closer to the actual Sun than
other simplified atmospheres. MHD simulations have been used to test the
reliability of inversion techniques (Orozco Suárez et al., 2010; de la Cruz
Rodríguez et al., 2012). In the first of these works, a confirmation of the
predictions by Sánchez Almeida et al. (1996) was found: if your inference
technique assumes magnetic fields and velocities constant with depth and you
use it on data coming from an atmosphere where these quantities are depth
dependent, the result is just the average of the actual stratification
weighted with the generalized response function to perturbations of that
quantity. In the second of these papers, the non-LTE inversion code called
NICOLE is tested. Here we report on preliminary results by Harker et al.
(2016) to illustrate the role of simulations as a tool to determine the
optimum wavelength sample of Stokes profiles. This is a particularly
interesting topic that is very relevant to the observational (and hence
interpretational) work. Are the available samples enough for capturing all the
information encoded in the Stokes profiles? What is the optimum sampling one
should use with a new instrument under development depending on the goals such
an instrument aims to fulfill?

Figure 15 shows the Stokes $I$ (in depression), $Q$, $U$, and $V$ profiles of
the Fe i line at 630.25 nm across a slit over a sunspot simulation by Rempel
(2012) (left column panels) and their corresponding power spectra (right
column panels). The simulation contains the transition from the quiet Sun (at
both sides of the $X$ dimension) through the penumbra and the umbra of a
sunspot. The power spectra of Stokes $V$, $Q$, and $U$ are wider than that of
Stokes $I$ as a natural consequence of their shapes. Therefore, a cut-off
frequency is better found in the polarization profiles. In this example,
Shanon’s critical sampling interval is around 1.25 pm/pixel with the
remarkable fact that no convolution with an instrumental PSF has been applied
to the profiles. This value is coarser than that provided by several ground-
based spectrographs with resolutions about $R\simeq 10^{6}$ that would be
considered too fine for the required diagnostics. As soon as the profiles are
observed by an instrument with a finite width PSF, the Nyquist frequency will
shrink to smaller values. Hence, the critical sampling will be coarser. These
kinds of calculations can therefore help in designing new instruments and, as
far as this paper is concerned, in deciding the adequate spectral sampling for
any synthesis or inversion code: if the sampling is too fine, we are wasting
computational time; if the sampling is too coarse, we are neglecting available
information.

Figure 15: Stokes profiles ($I$ is in line depression) across a sunspot
simulation by Rempel (2012) (left column) and power spectra (right column).
Adapted from Harker et al. (2016).

##  6 Response functions

The needle in an old-fashioned ammeter is useful as far as it moves when a
current circulates through the electric circuit. If, for instance, the needle
does not move when a given very weak or very strong current passes through,
then the ammeter is useless or _insensitive_ to those intensities. Keeping
this analogy, our observables, the Stokes profiles are useful for inferring
given physical quantities as long as they change when those physical
quantities vary. Of course, to be detectable, the change should be larger than
the noise. Therefore, the correct question to ask of a given
spectropolarimetric proxy whether it senses, for example, $T$, $B$, or $v_{\rm
LOS}$ is how much it modifies when $T$, $B$, or $v_{\rm LOS}$ change. The
preceding section has been an attempt in that direction: we have been changing
the various atmospheric quantities and checking the modifications in the
Stokes profiles. We have proceeded in the _direct_ way, that is, through the
solution of the _differential_ RTE. This direct approach is not very useful in
practice, though, as we already recognized in Section 5.2. Which quantity is
to be modified first, at which optical depths, and by how much? If the problem
is the simplest we talked about in the introduction (that of measuring
velocities from the line core wavelength), then there may be some room for the
direct approach. If not, the diagnostic capabilities of the Stokes profiles
have to be further explored with the final goal of proceeding the _inverse_
way, that is, of solving the _integral_ equation known as the formal solution
of the RTE (Eq. 8).

Here, we have a difficult problem where the observables depend nonlinearly on
the unknowns. The nonlinear character is clear: the observables —on the left-
hand side of Eq. (8)— are equal to an integral of the product of three terms,
each depending strongly non-linearly on the physical quantities that
characterize the model atmosphere (e.g., Del Toro Iniesta, 2003b). Changes in
the Stokes spectrum are then very difficult to predict when modifications in
the physical parameters occur. As in many other branches of physics, the
diagnostic tools come out trough a linearization analysis. We can assume, for
instance, that in a very special regime, when perturbations are small enough,
changes occur linearly. These are the basics of linearization that, in the
realm of solar physics were introduced for non-polarized light by Mein (1971,
see also , , and , ) through the so-called weighting functions, although the
name of _response functions_ (RFs) did not appear in the literature until the
work by Beckers and Milkey (1975). Since polarization was not taken into
account, those analyses were only strictly valid for isotropic media or, as
far as we are concerned, for non-magnetic atmospheres. Response functions were
introduced within polarized radiative transfer by the brothers Landi
Degl’Innocenti and Landi Degl’Innocenti (1977) but RFs were paid little
attention to until the works by Landi Degl’Innocenti and Landolfi (1982,
1983), Grossmann-Doerth et al. (1988), Sánchez Almeida (1992), and Ruiz Cobo
and Del Toro Iniesta (1992, 1994). The latter found that the perturbations
$\delta x_{i}$ exerted to the $p+r$ atmospheric quantities ($p$ of them
varying with height and $r$ constant) induce modifications
$\delta\mbox{\boldmath$I$}(0)$ to the observed Stokes profile given by

$\delta\mbox{\boldmath$I$}(0)=\sum_{i=1}^{p+r}\int_{0}^{\infty}\mbox{\boldmath$R$}_{i}(\tau_{\rm
c})\,\delta x_{i}(\tau_{\rm c})\,{\rm d}\tau_{\rm c},$ (29)

where

$\mbox{\boldmath$R$}_{i}(\tau_{\rm c})\equiv\mbox{\bf O}(0,\tau_{\rm
c})\,\left[\mbox{\bf K}(\tau_{\rm
c})\frac{\partial\mbox{\boldmath$S$}}{\partial x_{i}}-\frac{\partial\mbox{\bf
K}}{\partial x_{i}}(\mbox{\boldmath$I$}-\mbox{\boldmath$S$})\right].$ (30)

Therefore, the modification of $\mbox{\boldmath$I$}(0)$ is given by a sum of
terms, each related to one of the atmospheric quantities characterizing the
medium. The terms are integrals over the whole atmosphere of the model
atmospheric quantities weighted by the RFs. The physical meaning in Eq. (29)
is straightforward: imagine that we change only a given quantity ($T$, $B$,
$v_{\rm LOS}$, or any other) with a magnitude unity (i.e., 1 K, 1 G, 1  km
s-1, etc.) in the narrow surroundings of a given continuum optical depth
$\tau_{0}$; then, the subsequent modification in the emergent Stokes spectrum
is just the value of the corresponding RF at that optical depth:

$\delta\mbox{\boldmath$I$}(0)=\mbox{\boldmath$R$}_{i}(\tau_{0}).$ (31)

Then, since the Stokes profiles are usually recorded normalized to some
reference value (e.g., the average, —unpolarized— continuum intensity of the
quiet Sun), units for RFs are inverse units of the corresponding quantity.
That is, the response function to perturbations of temperature is measured in
K-1; the response to perturbations in the magnetic field strength is measured
in G-1; and so on. Thus, a response function can be defined as the
modification that the Stokes spectrum experiences when the medium undergoes a
unit perturbation of a given physical quantity at a given very narrow region
in optical depth. Equation (30) tells us that these modifications build upon
the variations of the propagation matrix and the source function vector with
respect to the physical quantities and their evolution through the atmosphere
as driven by the evolution operator. The two variations have an opposite sign.
This means that they are somehow competing as one could expect. While $S$
represents the sources of photons, K represents the sinks. Indeed we know that
we do not only speak about photon removal but also about pleochroism and
dispersion but, certainly, the propagation matrix role is somehow similar to
that of a withdrawal. The counterbalancing between $S$ and K is very important
to understanding radiative transfer because some analyses forget it and only
account for the effects of K (absorption in the non-polarized case). This was
clearly pointed out and explained by Ruiz Cobo and Del Toro Iniesta (1994).

Equation (29) suggests that RFs play the role of partial derivatives of the
observed Stokes profiles with respect to the atmospheric quantities once they
have been discretized. This role is even more clear when we go down to the
real world of a quadrature formula for that equation. Model atmospheres are
usually described numerically by a grid of points that are spaced in
logarithmic optical depth. Let $\Delta(\log\tau_{\rm c})$ be that spacing. If
we call $x_{i,j}\equiv x_{i}(\tau_{j})$ and
$\mbox{\boldmath$R$}_{i,j}\equiv\mbox{\boldmath$R$}_{i}(\tau_{j})$, then Eq.
(29) can be written as

$\delta\mbox{\boldmath$I$}(0)=\sum_{i=1}^{p}\sum_{j=1}^{n}a_{j}\mbox{\boldmath$R$}_{i,j}\,\delta
x_{i,j}+\sum_{k=1}^{r}\mbox{\boldmath$R$}^{\prime}_{k}\,\delta x_{k},$ (32)

where $a_{j}=\Delta(\log\tau_{\rm c})\ln 10\,c_{j}\tau_{j}$, with $c_{j}$
being the quadrature coefficients. Therefore, if we include $a_{j}$ in the
RFs, as one usually does in graphical representations, then Eq. (32) shows the
Stokes spectrum modifications as _linear_ expansions of the new variables
$x_{i,j}$ and $x_{k}$. The first term on the right-hand side corresponds to
those physical quantities that vary with depth; the second stands for those
that are assumed to be constant.141414For the specific meaning of the RF to
perturbations of a constant quantity, $\mbox{\boldmath$R$}^{\prime}$, see
Sect. 6.1 below. In summary, we can say that RFs are indeed partial
derivatives of $\mbox{\boldmath$I$}(0)$ with respect to the (numerical)
atmospheric parameters and, thus, they directly provide the sensitivities of
the Stokes spectrum to perturbations of the physical conditions in the medium.
Examples of these RFs are plotted in Figs. 16, 17, 18, and 19. They have been
evaluated for Stokes $I$, $Q$, $U$, and $V$, respectively, of the Fe i line at
630.25 nm to perturbations of the temperature (top row panels), of the
magnetic field strength (middle row panels), and of the LOS velocity (bottom
row panels). The RF values are multiplied by $10^{6}\,$K-1, $10^{6}\,$G-1, and
$10^{4}\,(\\!\mbox{$\>$km$\,$s${}^{-1}$})^{-1}$. The two columns correspond to
two different model atmospheres. That in the left-hand columns has the
temperature stratification of the HSRA model (Gingerich et al., 1971), a
constant $B=2000$ G, $\gamma=30\mbox{${}^{\circ}$}$, and
$\varphi=60\mbox{${}^{\circ}$}$; the plasma is at rest in this model. That in
the right-hand columns has a 500 K cooler temperature, and a magnetic field
500 G weaker, 20∘less inclined, and an azimuth of 10∘; $v_{\rm
LOS}=1.58+0.3\log\tau_{\rm c}$ ( km s-1).

Equation (32) hints at a way for calculating RFs through what could be called
_the brute force method_. This method is a four-step procedure: 1) synthesis
of the Stokes spectrum in a given model atmosphere; 2) perturbation of just
one of the (numerical) atmospheric parameters by a small amount and synthesis
of the spectrum in the new model atmosphere; 3) calculation of the ratio
between the difference of the two spectra and the perturbation; 4) repetition
of steps 2) and 3) for each optical depth, for each wavelength sample, and for
the remaining Stokes parameters. This is a formidable calculation as soon as
the number of free parameters is large. Fortunately, Eq. (30) provides a
shortcut since the evolution operator, the propagations matrix, and the source
function vector have to be calculated anyway in every synthesis of the
spectrum. With only the added calculations of the derivatives, one can easily
calculate RFs at the same time as the RTE is solved. This property is
extremely useful for inversion codes, as we shall see in Section 7.

Equation (32) also offers an explicit explanation of the astrophysical ill
conditioning we commented on in Section 1: the same modification of
$\mbox{\boldmath$I$}(0)$ may be produced by perturbations of different
quantities or by perturbations of a single physical quantity but at several
optical depths. That is, the effects of temperature can be similar to those of
the magnetic field strength or the effects of perturbing $B$ at $\log\tau_{\rm
c}=-0.5$ can be the same as those of perturbing $B$ at $\log\tau_{\rm c}=-3$.
Therefore, we cannot say that the changes $\delta\mbox{\boldmath$I$}(0)$ are
produced by perturbations of this physical parameter or that other without
considering all of them at the same time. Cross-talk among some parameters may
appear and, then, the retrieval of those parameters will be less reliable
(see, e.g., Section 6.2).

###  6.1 Properties of response functions

A glance at Figs. 16, 17, 18, and 19 readily tells us that some RFs are bigger
than the others. This means that our line is more sensitive to some physical
quantities than to others. However, the fact that RFs are measured in inverse
units of those for their corresponding parameters makes it difficult to
compare their relative sensitivity. In this regard, _relative_ RFs shed some
light. If we consider relative perturbations $\delta x_{i,j}/x_{i,j}$, then we
can define relative RFs
$\tilde{\mbox{\boldmath$R$}}_{i,j}\equiv\mbox{\boldmath$R$}_{i,j}x_{i,j}$.
Hence, $\tilde{\mbox{\boldmath$R$}}_{i,j}$ speak about the response of the
Stokes spectrum to relative (i.e., dimensionless) perturbations. Experience
shows that relative RFs to $T$ perturbations are the biggest at all depths and
wavelengths, clearly indicating that temperature is the most important
quantity in line formation. (Indeed, temperature is _the_ physical quantity
that governs the thermodynamical state of the material medium because we
assume that hydrostatic equilibrium prevails throughout our model atmospheres.
After this assumption, pressure, the necessary second thermodynamic variable
gets automatically prescribed.) Response functions to temperature
perturbations start being different from zero at the deepest layers when
compared to the remaining quantities. This is because the second term in the
right-hand side of Eq. (30) goes to zero as the continuum optical depth tends
to infinity. This was explained by Ruiz Cobo and Del Toro Iniesta (1994). This
physical fact implies that spectral lines tend to be insensitive at these low
layers to the other physical quantities.

The Stokes profile wavelength symmetries are preserved in RFs: in the absence
of velocity gradients, RFs of Stokes $I$, $Q$, and $U$ to any perturbation are
even functions of wavelength and RFs of Stokes $V$ are odd. This means that,
in fact, velocity gradients increase the diagnostic capabilities of spectral
lines. In their absence, half of the profile is useless since the information
they provide is exactly the same as the other half.151515Indeed, such
redundant information should help in decreasing the uncertainties in the
retrievals by reducing the noise by a factor two.

For given purposes, we can conceive constant perturbations with depth, in
spite of the quantity being depth dependent. Owing to their nature, some
physical quantities may be assumed constant with depth (e.g., macro- and
microturbulent velocity, or any of the ME free parameters). In such cases,
constant perturbations are in order. If so, then, Eq. (29) tells us that the
resulting modification in the Stokes spectrum is given by the product of such
a constant perturbation times the integral of the corresponding RF over the
whole atmosphere. Hence, we can say that the RF to a constant perturbation,
$\mbox{\boldmath$R$}^{\prime}$ in Eq. (32), is directly the integral of the
regular response function or, in numerical terms,

$\mbox{\boldmath$R$}^{\prime}_{k}\equiv\sum_{j=1}^{n}a_{j}\mbox{\boldmath$R$}_{k,j}.$
(33)

As shown by Ruiz Cobo and Del Toro Iniesta (1994), RFs play the role of a PSF
in the general theory of linear systems. Under this general theory, our system
—the Stokes spectrum— experiences an input (the perturbation) and provides an
output, $\delta\mbox{\boldmath$I$}(0)$. If the input is a Dirac delta, then
the output is the corresponding value of the response function. If the input
is harmonic throughout the atmosphere, then the response is the Fourier
transform of the RF.

Response functions are model dependent. This property is extremely important
in our understanding of spectral line sensitivities and, thus, in the
inversion of the RTE. Instead of being a drawback, such a model dependence
helps in disentangling the effects produced by the different quantities in
distinct model atmospheres. Even in fixed model atmospheres, it is very
difficult to discard one quantity or the other at once, however, and most of
them have to be retrieved at the same time. Once this is carried out, one can
theoretically understand the meaning of measurements (Sánchez Almeida et al.,
1996, Del Toro Iniesta, 2003b).

Figure 16: RFs of Stokes $I$ to perturbations of $T$ (top panels), $B$ (middle
panels), and $v_{\rm LOS}$ (bottom panels). Units are $10^{-6}\,$K-1,
$10^{-6}\,$G-1, and $10^{-4}\,(\\!\mbox{$\>$km$\,$s${}^{-1}$})^{-1}$.The two
columns correspond to two different model atmospheres. That in the left-hand
column has the temperature stratification of the HSRA model (Gingerich et al.,
1971), a constant $B=2000$ G, $\gamma=30\mbox{${}^{\circ}$}$, and
$\varphi=60\mbox{${}^{\circ}$}$; the plasma is at rest in this model. That in
the right-hand column has a 500 K cooler temperature, and a magnetic field 500
G weaker, 20∘less inclined, and an azimuth of 10∘; $v_{\rm
LOS}=1.58+0.3\log\tau_{\rm c}$ ( km s-1).

Figure 17: Same as Fig. 16 for Stokes $Q$.

Figure 18: Same as Fig. 16 for Stokes $U$.

Figure 19: Same as Fig. 16 for Stokes $V$.

Last, but not least, the linear nature of RFs allows us to generalize them to
any linear combination of Stokes profile wavelength samples. This property
helps us understand what can be extracted from different proxies that are
traditional in solar polarimetry but, most importantly, helps the astronomer
in taking influence of the instrument into account. Since most instruments act
as linear systems on light, the detected spectrum is a convolution of the
actual spectrum with the instrument spectral PSF. Convolution is linear and,
thus, one can easily conclude (Ruiz Cobo and Del Toro Iniesta, 1994) that the
RFs of the convolved spectrum are nothing but those of the original spectrum
convolved as well with the PSF.161616Convolutions with the spatial PSF of the
instrument are also taken into account in modern inversions that take spatial
degradation into account. See Section 7.5.

###  6.2 Analytic response functions

When all the terms on the right-hand side of Eq. (30) can be calculated
analytically (see Sect. 2.3), RFs are necessarily analytic and then we can use
them to gain some physical insight into the diagnostic capabilities of the
Stokes spectrum about the physical quantities that characterize the medium.
This is the case of the ME approximation, where all the quantities are
constant with depth. There, index $j$ in Eq. (32) drops and (after inclusion
of the coefficient into the RF) we can properly write:

$\mbox{\boldmath$R$}_{i}(\lambda)=\frac{\partial\mbox{\boldmath$I$}(\lambda)}{\partial
x_{i}}.$ (34)

That is, that RFs are strict partial derivatives of the Stokes spectrum with
respect to the free parameters of the problem (Orozco Suárez and Del Toro
Iniesta, 2007). In Eq. (34) we have removed the $\tau_{\rm c}=0$ indicator in
the emergent Stokes spectrum and, rather, we have made explicit the dependence
of RFs on wavelength. We have just seen in the former subsection that these
RFs are indeed integrals over $\tau_{\rm c}$ and, hence, the dependence on it
disappears. This analytic character is particularly important in the
controversial discussion about the possibility of disentangling $B$, from
$\alpha$ (the filling factor) in the case that our magnetic features are not
fully spatially resolved. In Fig. 20 we reproduce Fig. 3 from Del Toro Iniesta
et al. (2010). It shows RFs to (constant) perturbations in those two
quantities plus $v_{{\rm LOS}}$ in a ME atmosphere. As one can clearly see,
both Stokes $V$ RFs to $\alpha$ and to $B$ perturbations are almost
proportional among themselves and to the Stokes $V$ profile itself. This can
easily be traced back to the expected behavior from the weak field
approximation, as expressed in Equation (24). From this proportionality we
should conclude that it is indeed very difficult to discern the values of $B$
and $\alpha$ separately. If it were not for Stokes $I$, the widespread belief
that only $\alpha B$ or $\alpha B\cos\gamma$ can be retrieved would be true.
However, the Stokes $I$ RFs to $\alpha$ and $B$ are neatly different from each
other and this necessarily implies that we have the means of inferring the two
quantities independently. Stokes $Q$ and $U$ also help in disentangling the
field strength and the filling factor for similar reasons as soon as they are
above the noise level.

Figure 20: Stokes I (left panel) and V (right panel) RFs to $v_{\rm LOS}$
(solid, black lines), to $B$ (dotted, blue lines), and to $\alpha$ (dashed,
red lines) in the weak-field case. Perturbations of 10  m s-1for $v_{\rm
LOS}$, of 10 G for $B$, and of 0.1 for $\alpha$ have been assumed. Adapted
from Del Toro Iniesta et al. (2010).

Among other features of RFs, the latter authors showed how the thermodynamical
parameters of the ME atmosphere can have cross-talk among themselves: their
RFs are fairly similar in shape, so that their effects can be misinterpreted
by the inversion codes. Notably, the RFs to perturbations of $B$, $\gamma$,
$\varphi$, and $v_{\rm LOS}$ are markedly different from one another and with
respect to those of $\eta_{0}$, $\Delta\lambda_{\rm D}$, and $a$. This
explains the good result of ME inversion codes in accurately retrieving the
magnetic and dynamic parameters while the thermodynamic parameters are
sometimes wrong. Our conclusion is also consistent with, and explains, the
findings by Lagg et al. (2004) who decided to leave the damping parameter
fixed with minor changes in the fitted magnetic and dynamic parameters while
$\Delta\lambda_{\rm D}$ and $\eta_{0}$ were significantly affected. Linearity
can also be useful for spatially coupled inversion techniques (see Sect. 7.5.1
below).

##  7 Inversion techniques

Once we have discussed all the ingredients and assumptions, we can face the
main problem in astrophysics, namely that of making theory and observations
compatible. In other words, we can face the inversion problem by deriving the
unknown physical quantities through comparison between observed and synthetic
Stokes profiles.171717The inversion problem could be thought of as the
“observational” part of the compatibility game. The “theoretical” part
involves the choice of hypotheses included in the physical scenario which,
according to Del Toro Iniesta et al. (2010), defines the assumed model
atmosphere. One has to decide whether the radiative transfer is LTE or NLTE,
whether the physical quantities are dependent on the optical depth, whether
macro- or microturbulence are needed, etc. All these assumptions settle the
theoretical framework of the problem. Figure 21 describes how the problem gets
complicated as compared to the mere forward problem in Figure 2.

Figure 21: Block diagram of the inversion problem under NLTE conditions.

One can clearly see how a new overarching loop is present that indicates the
needs for changing the model atmosphere if the synthetic Stokes spectra do not
properly fit the observed ones. The problem turns out to be formidable and
requires new, specific assumptions that make it tractable. In particular, some
of the quantities have to be calculated from the strict NLTE conditions (see
Section 7.2.3).

Even the simpler LTE problem gets complicated, and indeed becomes iterative,
regardless of acknowledging the stratification in the physical quantities of
the atmosphere (see Figure 22). The needs for modifying such a model
atmosphere according to the deviations between observed and synthetic profiles
makes a loop necessary after calculating both the synthetic spectra and their
derivatives with respect to the free parameters. Fortunately, we know how to
calculate these derivatives through RFs at the same time as we synthesize the
Stokes spectrum with little extra computational effort.

Figure 22: Block diagram of the inversion problem under LTE conditions.

Looking for convergence means measuring the distance between observed and
synthetic profiles in the space of observables. Any inversion procedure must
have a threshold below which the user can consider that convergence has been
reached because the fit cannot be further improved within the current
assumptions.

###  7.1 Topology in the space of observables

The topological problem depicted in the Introduction for the inversion problem
(or any astrophysical inference) needs to be substantiated in minimizing a
distance: a metric in the space of the observables. Since Stokes $I_{\rm d}$,
$Q$, $U$, and $V$ belong to ${\mathchoice{\rm I\mskip-4.0muL}{\rm
I\mskip-4.0muL}{\rm I\mskip-4.5muL}{\rm IRmskip-5muL}}^{2}$, the quadratic
norm of the difference turns out to be the natural distance between any two
profiles. We want to approximate two sets of profiles, so that all four Stokes
parameters should be taken into account. Therefore, when in practice we deal
with discrete samples, the sought-for distance can be written as

$\chi^{2}(\mbox{\boldmath$x$})\equiv\frac{1}{N_{\rm
f}}\sum_{s=0}^{3}\sum_{i=1}^{q}\left[I_{s}^{\rm obs}(\lambda_{i})-I_{s}^{\rm
syn}(\lambda_{i};\mbox{\boldmath$x$})\right]^{2}\,w_{s,i}^{2},$ (35)

where index $s$ runs for the four Stokes parameters, we assume $q$ wavelength
samples, and $N_{\rm f}$ stands for the number of degrees of freedom, that is,
the difference between the number of observables ($4q$) and that of the free
parameters (the number of elements in $x$, $np+r$; see Section 2.1 and
Equation 9).

$\chi^{2}(\mbox{\boldmath$x$})$ is a merit function of the atmospheric
quantities that measures the distance between the observed and the synthetic
profiles and has to be minimized in order to achieve a good fit. Having a
normalized merit function to the degrees of freedom is useful to warn the user
not to use an unreasonably large number of free parameters as compared with
the number of observables. In such a case, $\chi^{2}$ would turn out to be
always too big. The _weights_ $w_{s,i}$ can be used to favor some data more
than the others. For instance, one can set them to the inverse of the
measurement errors. For many applications they are simply kept at unity.

We can look at $\chi^{2}(\mbox{\boldmath$x$})$ as a scalar field in an
($np+r$)-dimensional space. Since the number of dimensions may be too large,
the minimization problem may turn out to be intractable. Before going to
specific techniques that make it affordable, let us consider the paths through
which we can look for the minimum of the merit function. That is, we have to
find the derivatives of $\chi^{2}$ with respect to the atmospheric free
parameters. Ruiz Cobo and Del Toro Iniesta (1992, see , , and ,  as well)
showed that such derivatives are directly given by the RFs:

$\frac{\partial\chi^{2}}{\partial x_{m}}=\frac{2}{N_{\rm
f}}\sum_{s=0}^{3}\sum_{i=1}^{q}\left[I_{s}^{\rm obs}(\lambda_{i})-I_{s}^{\rm
syn}(\lambda_{i};\mbox{\boldmath$x$})\right]\,w_{s,i}^{2}\,R_{m,s}(\lambda_{i}),$
(36)

where, for the sake of a more compact notation, index $m$ runs from 1 to
$np+r$ (including constant and variable physical quantities), the quadrature
coefficients in Eq. (32) are assumed to be included in the RFs when needed,
and no distinction is made between $R$’s and $\mbox{\boldmath$R$}^{\prime}$’s.
The same authors also demonstrated that the second derivatives can be
approximated by

$\frac{\partial^{2}\chi^{2}}{\partial x_{m}\partial
x_{k}}\simeq\frac{2}{N_{\rm
f}}\sum_{s=0}^{3}\sum_{i=1}^{q}w_{s,i}^{2}\,\left[R_{m,s}(\lambda_{i})\,R_{k,s}(\lambda_{i})\right].$
(37)

Regardless of the way we approach its minimum in the hyperspace of parameters,
that $\chi^{2}$ is the natural metric is reinforced by the fact that other
metrics have been tried (e.g. Lagg et al., 2004) that have finally converged
to almost the same formulation as in Eq. (35) (Lagg et al., 2007).