Issue 
A&A
Volume 629, September 2019



Article Number  A24  
Number of page(s)  16  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201935692  
Published online  02 September 2019 
FIRTEZdz
A forward and inverse solver of the polarized radiative transfer equation under Zeeman regime in geometrical scale
^{1}
Leibniz Institut für Sonnenphysik, Schöneckstr. 6, 79104 Freiburg, Germany
email: apy@leibnizkis.de
^{2}
Instituto de Astrofísica de Canarias, C. Vía Láctea s/n, 38205 La Laguna, Tenerife, Spain
^{3}
Departamento de Astrofísica de Canarias, Universidad de La Laguna, Avda Astrofísico Sánchez s/n, 38206 La Laguna, Tenerife, Spain
Received:
13
April
2019
Accepted:
9
July
2019
We present a numerical code that solves the forward and inverse problem of the polarized radiative transfer equation in geometrical scale under the Zeeman regime. The code is fully parallelized, making it able to easily handle large observational and simulated datasets. We checked the reliability of the forward and inverse modules through different examples. In particular, we show that even when properly inferring various physical parameters (temperature, magnetic field components, and lineofsight velocity) in optical depth, their reliability in heightscale depends on the accuracy with which the gaspressure or density are known. The code is made publicly available as a tool to solve the radiative transfer equation and perform the inverse solution treating each pixel independently. An important feature of this code, that will be exploited in the future, is that working in geometricalscale allows for the direct calculation of spatial derivatives, which are usually required in order to estimate the gas pressure and/or density via the momentum equation in a threedimensional volume, in particular the threedimensional Lorenz force.
Key words: methods: data analysis / methods: numerical / techniques: polarimetric / polarization / radiative transfer
© ESO 2019
1. Introduction
In observational sciences, meaningful physical properties of an object are inferred from available observations. In the context of the stellar atmospheres, and in particular for this work, the solar photosphere, most of this inference is performed by means of different approximations to the solution of the radiative transfer equation (hereinafter referred to as RTE) applied to the spectra and polarized spectra emitted by the object.
The RTE characterizes the interaction of the matter with the light that traverses through it. This interaction leaves its “fingerprints” in many different ways, for instance in the shape of the continuum of the spectra coming from the object, or in the presence or absence and the shape of spectral lines. Properly modeled, these features allow approximate characterization of the studied object. In recent decades, the advent of polarimetric observations has allowed the generalization of the RTE to polarized light (Unno 1956; Rachkovsky 1967; Landi Degl’Innocenti & Landi Degl’Innocenti 1972). This step was of great significance, as the polarization of the light encodes information on the magnetic field configuration as well as the geometry of the system, which allows a deeper insight into the system being studied.
In solar physics, there are different approaches to solving the polarized RTE depending on the magnetic field regime adopted: atomic, Hanle, Zeeman, and/or Paschen (see e.g., Landi Degl’Innocenti & Landolfi 2004). In this work, we focus on the solution of the polarimetric RTE under the Zeeman magnetic regime, which is usually the dominant one at this atmosphere layer.
There are already several methods to address the solution of the RTE, each with its respective accuracy (see, e.g., del Toro Iniesta & Ruiz Cobo 2016). There are several computational codes that allow one to infer the optical depth (τ) dependence of the physical properties in the solar atmosphere; for example, SIR (Ruiz Cobo & del Toro Iniesta 1992), SPINOR (Frutiger et al. 2000), NICOLE (SocasNavarro et al. 2015), SNAPI (Milić & van Noort 2018), and STIC (de la Cruz Rodríguez et al. 2019). These are of particular relevance here as they allow the furthest insight into the properties of the solar atmosphere in a threedimensional volume. In these cases, the conversion between optical depth τ and geometrical height z is then performed by considering the relation dτ = −ρκ_{c}dz, where ρ is the density and κ_{c} the continuum opacity, with the latter being a nonlinear function of the temperature T and gas pressure P_{g}. Since the proper calculation of the gas pressure (and density) involves the whole threedimensional volume, for which there is no account in the aforementioned codes, the gas pressure and density inferred by them are unrealistic, and hence also the conversion between τ and z. Deviations from hydrostatic equilibrium involve, among others, timevariation of the velocity, advection and viscous terms and, the most relevant at this atmospheric layer, the Lorentz force. The latter is responsible for the Wilson depression, a depression of the height at which photospheric spectral lines are formed with respect to the quietSun because of the presence of magnetic field. Because the geometrical scale z is then inaccurate, so is the calculation of the spatial derivatives, such as those involved in the determination of ∇ ⋅ B = 0 and the electric currents . Consequently, it is not possible to establish whether the inferred magnetic field verifies Maxwell’s equations.
Up to now there have only been a handful of attempts to address these problems and limitations. One of them was presented by Puschmann et al. (2010), who used the SIR code to infer the various physical parameters in a τscale before transforming these parameters into a selfconsistent geometricalheight scale z in a threedimensional volume. To that end, these latter authors used a more detailed gas pressure (P_{g}) and density (ρ) than the one given by the SIR code as they were calculated by minimizing the static equation of momentum in ideal magnetohydrodynamics (MHD) including the Lorentz force, together with the ∇ ⋅ B = 0 condition. In doing so, the new gas pressure and density are not necessarily compatible with the observed polarization spectra from which the physical parameters were inferred. Another possibility to gain access to the geometrical scale z was given by Riethmüller et al. (2017), where threedimensional MHD simulations were used to reproduce the polarization signals of an observed area on the solar surface. Since an MHD code was employed to obtain the physical parameters, the gas pressure and density verify the MHD equation of motion in the whole of the domain employed. Moreover, physical parameters are already prescribed on the geometrical scale, and thus no τ → z conversion is needed. However, this approach requires a huge computational effort as several MHD simulations have to be run. Also, under this approach, the result of the inversion might be accurate as long as the actual physical model giving rise to the observed spectra is attained in the simulation.
Here, we put forward an alternative way to tackle the aforementioned problems. If one solves the RTE in geometrical scale, then spatial derivatives are handled directly. This strategy allows the straightforward inclusion of the Lorentz force terms in the equation of motion as well as the null divergence condition in the magnetic field to be fulfilled: ∇ ⋅ B = 0. Thus it would be possible to solve the RTE selfconsistently with physical constraints in three dimensions. Here we present a first step in this approach, and we implement a onedimensional polarized RTE solver for the Zeeman regime in the geometrical scale z. The code is parallelized and can be used to calculate synthetic polarized spectra from MHD simulations and invert observed spectropolarimetric data in order to infer the various physical parameters (temperature, three magnetic components of the magnetic field, and lineofsight (LOS) velocity) of the solar atmosphere. The inference of these parameters is accurate in terms of τ and depends on the ability of the user to infer the gas pressure and density in geometrical scale.
2. Numerical model
2.1. Forward modeling
This numerical code faces the problem of solving the RTE in a geometrical scale under the assumption of local thermodynamic equilibrium (LTE), complete redistribution approximation, and a plane parallel atmosphere. This problem is a coupled system of equations that can be expressed as:
where s is the geometrical distance along the ray path, 𝕀(λ, s) is a wavelengthdependent pseudovector of the four Stokes parameters characterizing the polarization state of the light, is the socalled propagation matrix, a wavelengthdependent 4 × 4 matrix that characterizes the interaction between light and the medium it travels through, and 𝕊(λ, s) is a wavelengthdependent pseudovector for the source function. In the following, and for the sake of clarity, we avoid the explicit dependence in λ.
The derivation of the solution of the polarized RTE in the Zeeman regime has been covered extensively in the literature (see e.g., Landi Degl’Innocenti & Landolfi 2004, and references therein) and several codes have been implemented to solve it in τscale. Here we do not revisit those derivations. Instead, we outline the general building blocks of the formalism employed for the resolution of this problem in geometrical scale.
This equation has a formal solution in terms of a formal operator 𝕆 (Landi Degl’Innocenti & Landi Degl’Innocenti 1985):
where s_{0} is the lower boundary, is the Stokesevolution operator between the lower boundary and the integration limit, and 𝕀(s_{0}) is the illumination from the lower boundary. This expression evaluates the leaving polarimetric properties of the light given the lower boundary and the physical properties of the medium. In this work, we consider that the medium through which light propagates is made of N slabs, each of which has constant physical properties along the whole slab. Under such an assumption, 𝕆(s, s′) has an analytic expression (see Eq. (10), Landi Degl’Innocenti & Landi Degl’Innocenti 1985). In this case, the light leaving one of these slabs (index k) is given, following Eq. (2), by
where 11 is the 4 × 4 identity matrix, 𝕀(s_{k}) stands for the Stokes vector leaving slab k, 𝕊(s_{k}) is the source function in slab k, which is constant throughout the slab, and 𝕀(s_{k − 1}) is the incident Stokes vector from below the slab. In the general case, we consider a set of N such slabs with a proper boundary condition. In this case, it can be proven that the light leaving this medium is given by the concatenation of the various slabs:
Equation (4) expresses the light leaving the medium in terms of the physical properties of each slab (by means of the Stokesevolution operator and the emissivity in each slab) and on the lower boundary condition 𝕀(s_{0}). This equation is equivalent to that presented by del Toro Iniesta et al. (1995, Eq. (16)).
For the boundary condition 𝕀(s_{0}), we note that in a deepenough layer the temperature, the gas pressure, and the density are high enough to cause the opacity to be sufficiently large so as to ensure that the Stokesevolution operator becomes negligible. Therefore, from Eq. (3), if the Stokesevolution operator is small enough, from that deepenough layer, the leaving beam of such a deep layer depends only in its own physical properties through the source function 𝕀(s_{0}) = 𝕊(s_{0}). As we are assuming LTE, 𝕊(s_{k}) = (B(T(s_{k})), 0, 0, 0)^{†}, where B(T(s_{k})) is the Planck function for the temperature at s_{k}.
With this formalism one can solve the forward problem for the RTE. To do so, one needs to specify some additional features:

An equation of state (EOS): This is set to be given by an ideal gas equation with a variable molecular weight to account for free electrons.

Continuum opacities: We follow a reduced scheme of the socalled NICOLE opacity package (see Sect. 4.1.3 in SocasNavarro et al. 2015), excluding the contributions from Rayleigh scattering on molecular hydrogen and the contributions of boundfree and freefree terms due to magnesium.

Line opacity: This is calculated following Landi Degl’Innocenti (1976), excluding the nonLTE departure coefficients as here we restrict ourselves to the LTE case.

Partition functions: These are implemented following Wittmann (1974) who allowed analytical expressions for both the partition functions and their derivatives with respect to the temperature. This set of functions includes 92 species and 3 ionization states. We note that for the hydrogen atom, the first ionization stage involves the negative hydrogen.

Line broadening due to collisions: Two different treatments are available: (1) The ABO theory (Anstee & O’Mara 1995; Barklem et al. 1998), which requires additional parametric information to estimate the collisional broadening, namely the cross section (σ) and the velocity parameter (α). These can be supplied or they can be internally calculated (provided the energy limit for each energy level involved in the transition is known) from the tabulated tables (Anstee & O’Mara 1995; Barklem & O’Mara 1997; Barklem et al. 1998). (2) In the absence of the aforementioned data, the classic Unsöld approximation (Unsold 1955) is used. In this case, the line broadening is calculated without any additional atomic parameters. However, this formalism is known to underestimate the collisional broadening giving rise to narrower spectral lines (Barklem 1998).

Magnetic field effect: At this point the code takes into account only the Zeeman effect for the Stokes parameters, following the formalism of Landi Degl’Innocenti (1976) but changing the sign of Stokes V so that the code possesses the same sign criteria for the polarity of magnetic field as SIR, SPINOR, or NICOLE.

Spectral line blends are treated selfconsistently.
With these ingredients one is able to perform the forward modeling for any solar or stellar atmosphere prescribed in N slabs in geometrical scale. Furthermore, given this parametric model in geometrical scale, it is possible to face the inverse solver: namely, to infer the various physical parameters as a function of the geometrical scale given a set of spectropolarimetric observations. To the best of our knowledge, this is the first code capable of doing so. It is worth emphasizing that this RTE solver, as prescribed here, has analytical expressions for all the model parameters. This is extremely advantageous because it allows one to obtain analytical derivatives of the spectropolarimetric spectra (i.e., Stokes pseudovector) with respect to the various physical properties of the medium. These are the socalled response functions, and are of critical importance as they are usually required in the resolution of the inverse problem.
2.2. Response functions
2.2.1. General considerations
Response functions give the variation of the leaving Stokes pseudovector, 𝕀(s_{N}), due to the variation of any physical parameter x at a given height s_{k} (referred to as x_{k} = x(s_{k})). The various x_{k} refer to the temperature (T), the gas pressure (P_{g}), the density (ρ), the lineofsight velocity (v_{LOS}), and the three cartesian components of the magnetic field (B_{x}, B_{y}, and B_{z}) at the klayer of the atmospheric model. In our case, taking derivatives in Eq. (4):
We now take into account that the only parameters that depend on x_{k} are and 𝕊(s_{k}), and so it can be shown that
where δ_{k} stands for . As mentioned before, under the assumption of LTE, the source function at s_{k}, 𝕊(s_{k}), corresponds to the blackbody emission and depends only on the temperature at that height T_{k} = T(s_{k}). Therefore, derivatives of 𝕊(s_{k}) with respect to any other physical parameter other than the temperature vanish. The derivatives of the Stokesevolution operator with respect to the physical parameters can be determined by applying the chain rule to the expressions provided by Landi Degl’Innocenti & Landi Degl’Innocenti (1985, Eq. (10)). It can be shown that these derivatives depend only on the derivatives of the various elements of the propagation matrix : η_{I, Q, U, V}(s) and ρ_{Q, U, V}(s). For the sake of clarity and completeness we present them here in terms of the cartesian components of the magnetic field:
where κ_{c} and κ_{l} are the continuum and line opacity, respectively, B is the magnetic field modulus, and the various ϕ_{b, p, r} and ψ_{b, p, r} are the absorption and dispersion profiles, respectively, whose formulation can be found elsewhere (see e.g., Chapt. 5 of Landi Degl’Innocenti & Landolfi 2004).
From these equations, it can be shown that the derivatives of the various expressions in Eq. (7) with respect to v_{LOS}, B_{x}, B_{y}, and B_{z} are straightforward by means of the application of the chain rule and taking into account that κ_{c} and κ_{l} do not depend on these physical parameters. The derivatives involving thermodynamic properties, however, require additional consideration.
2.2.2. Thermodynamic derivatives of the propagation matrix elements
Here we focus on the calculation of the partial derivatives of the Stokes pseudovector with respect to the thermodynamic parameters: temperature, T, gas pressure, P_{g} and density ρ. These derivatives have to take into account that there is an additional constraint (EOS) that relates them to one another. Therefore, in contrast to the derivatives with respect to the lineofsight velocity, v_{LOS}, and the cartesian components of the magnetic field vector, B_{x}, B_{y}, and B_{z}, where all remaining parameters are assumed constant, for the derivatives with respect to T, P_{g}, and ρ, there are always two nonconstant parameters. For instance, if we have the derivative with respect to T, and we assume that ρ is constant; subsequently, P_{g} must change so that the EOS is fulfilled. For this reason, derivatives with respect to T, P_{g}, or ρ explicitly state which of them is assumed constant, that is, the Stokes derivative with respect to T at constant ρ is expressed as: . Therefore, here we are interested in obtaining
The selection of the thermodynamical parameter to be kept constant during the inversion is an arbitrary decision as they give equivalent results. It is important to bear in mind however that the computation of is computationally cheaper than that of , as the former does not involve changes in temperature, which is by far the most frequent parameter in the equations to be differentiated. For the derivatives with respect to T, there is no difference in using either or .
For the partial derivative , we need the calculation of δ_{k}𝕊(s_{k}) in Eq. (6). In addition, the partial derivatives with respect to the thermodynamic parameters in Eq. (8), require the calculation of the derivatives of the Stokes evolution operator: in Eq. (6). At the same time, these derivatives of the folloing Stokes evolution operator rely upon the calculation of δ_{k}{η_{I},η_{Q},η_{U},η_{V},ρ_{Q},η_{U},η_{V}}. In this particular case, these involve the determination of the derivatives of the continuum κ_{c} and line κ_{l} absorption coefficients, which in turn depend on the number of electrons per cubic centimeter (i.e., electron density) n_{e}. Therefore, after applying the chain rule over the Stokes evolution operator, what we need to evaluate is
These derivatives require additional algebra, yet they can be analytically calculated, as is shown in Appendices A–C.
2.3. Inverse solver
Up to here, we have faced the problem of obtaining the polarized spectra provided a set of heightdependent physical stratifications 𝕀(λ_{i})^{syn} = f(λ_{i}, X), where λ_{i} are the independent variables (the number of wavelength times the number of Stokes profiles observed, n_{ν}) and X the set of free parameters defining the forward model. It is sometimes interesting to proceed the other way around, that is, one records the polarized spectra from an object and aims at inferring some meaningful physical properties of the plasma from where the light comes from. This is called the inverse problem, where one looks for the model parameters that best fit, under some approximations, the observations 𝕀(λ_{i})^{obs} by defining a merit function:
where n_{f} stands for the number of free parameters and σ_{i} for the uncertainty and/or weighting of each observed data point. There are several strategies and algorithms when facing such an optimization problem, both local and globally. In this work we have chosen to implement a LM (Levenberg–Marquardt, Levenberg 1944 and Marquardt 1963) algorithm to perform this step. This kind of algorithm starts from a provided guess of the set of model parameters (X) and looks for the optimization by iteratively modifying these parameters. To do so, the perturbations of the free parameters are evaluated following
where δ is the perturbation of the set of the model free parameters, H ≈ J(λ)^{T} J(λ) is an approximation to the socalled Hessian matrix, J(λ) is the Jacobian matrix, and Λ is the dampening factor. Here it is worth emphasizing that the calculation of the perturbation requires the Jacobian matrix, which in our code is calculated analytically at the same time as the forward problem is solved by means of the response functions described in Sect. 2.2. It is important to note that the Jacobian makes use of dimensionless response functions, which are calculated by multiplying in Eq. (6) by a typical value of the physical parameter. We return to this particular point in Sect. 3.1. This point speeds up the calculations since otherwise they have to be numerically evaluated solving the forward model additional times for each iteration of the LM algorithm. Each time a perturbation, δ, to the model parameters is evaluated, the objective function of these new parameters is checked as a ruler to approve or reject the proposed perturbation.
An important point to notice here is that the first term on the righthandside of Eq. (11) may involve (as it does for instance in the examples presented in the following section) the inversion of large and quasisingular matrices. In order to tackle this point we follow the strategy introduced by Ruiz Cobo & del Toro Iniesta (1992) and use the Singular Value Decomposition (SVD) method^{1}. This way, it is possible to solve an illposed problem by reducing its dimensions. This reduction is controlled by the number of eigenvalues that one keeps when calculating the inverse dampenedHessian matrix. In this fashion, it is possible to find the optimum solution, in terms of leastsquares, for that number of eigenvalues. The choice of the threshold value is, for Stokes inversion codes, somewhere between 10^{−2} and 10^{−6} times the largest eigenvalue. In all examples in this paper we have used 10^{−4}.
The calculation on the uncertainties in the inference of the various X′ is performed assuming that each free parameter is independent of the others. In such a case, the covariance matrix is diagonal and given by the inverse of the diagonal elements of the Hessian. Therefore, the uncertainties are given by:
where . Since the uncertainties calculated this way neglect possible correlation between the various free parameters, the result is a lower limit for the uncertainties of each inferred parameter.
In principle, one can invert the whole set of physical parameters defining the atmosphere at once. However, this strategy can lead to inversion results that might not actually fit the profiles or at least, not as well as one would expect. This is because the minimization algorithm might get stuck in a local minimum which is not as deep as others might be. This is usually the case if one is not fairly close to the actual minimum and there are several free parameters (for instance 256 or 1280 as below). Instead, as introduced by Ruiz Cobo & del Toro Iniesta (1992), it is better to start first with a reduced problem where the number of free parameters is much lower than in the original problem and is then progressively increased as one approaches the minimum for that reduced problem. For this reason, we perform the fitting process in a sequential manner: we progressively increase the number of equivalent slabs perturbed, δ, from a constant in height perturbation for each parameter, until all the slabs are perturbed individually. Since response functions are given for all heights, we employ the concept of equivalent response functions (see e.g., Appendix A of Ruiz Cobo & del Toro Iniesta 1992) to linearly interpolate from the inversion grid into the model height grid at each iteration step.
2.4. Parallelization
Modern solar observing facilities such as Hinode (Kosugi et al. 2007), the Solar Dynamics Observatory (Scherrer et al. 2012), the Swedish Solar Telescope (Scharmer et al. 2003), the Goode Solar Telescope (Cao et al. 2010), and GREGOR (Schmidt et al. 2012) provide large amounts of data that currently can barely be processed at the needed speed (Borrero et al. 2011). In the future this issue will be aggravated with the arrival of new observatories such as the Daniel K. Inouye Solar Telescope (DKIST; Rimmele et al. 2010) and the European Solar Telescope (EST; Collados et al. 2010) whose data throughput will be several orders of magnitude larger than that of the current facilities. For this reason, it is currently compulsory to build new forward or inverse solvers for the RTE that are parallelized. Moreover, it is nowadays common practice for theoreticians and modelers, in particular developers of threedimensional magnetohydrodynamic codes such as MURAM (Rempel et al. 2009), CO5BOLD (Freytag et al. 2012), and MANCHA (GonzálezMorales et al. 2018), and even observers, to employ forward solvers for the polarized RTE to produce simulated observations from snapshots of the MHD simulations (Danilovic et al. 2008, 2010a,b; Borrero et al. 2010, 2014; Lites et al. 2017) in order to compare with real observations and to make predictions, so as to study under what circumstances a particular feature seen in the simulations could be detected observationally. The large number of available snapshots and the ever increasing size of the simulation boxes call also for fast (i.e., parallel) codes to solve the RTE.
With this in mind the forward and inverse code presented in this work has been conceived, from its initial stages, to being fully parallelized. The implementation is straightforward as each pixel on the observed solar surface (x, y) is treated independently and requires no information from neighboring pixels. Hence, we have chosen a masterslave parallelization scheme where the master process reads, distributes the work, and writes the outputs. Slave processes are used to solve the RTE and/or solve the inverse problem for each pixel independently.
3. Results
3.1. Forward modeling
Given the implementation detailed in the previous section, one can solve the RTE for polarized light under LTE and in the Zeeman regime for any atmospheric stratification. In Fig. 1 (top panel) we show the continuum intensity of a snapshot from the Rempel (2012) sunspot simulation. The size of the simulation box is 4096 × 512 × 256 with a grid size of Δx = 10, Δy = 10, and Δz = 8 km. We additionally (remaining panels in Fig. 1) show three different simulated Stokes profiles, 𝕀^{syn}(λ), together with their physical parameter stratification chosen from this snapshot. The synthetic profiles here, and all the forthcoming tests in this paper, are synthesised using the same spectral lines and wavelength grid. We used the two widely used Fe I spectral lines at 6301.5 and 6302.5 Å. The wavelength grid runs from −700 mÅ to +1800 mÅ with respect to the 6301.5 Å spectral line with a step of 5 mÅ. This setup is chosen so that the simultaneous inversion of the temperature, LOS velocity, and the three cartesian components of the magnetic field can be performed using the original height grid of the MHD simulation (256 heights). In a more realistic case, one can limit the number of free parameters setting the maximum number of heights at which perturbations of the model are calculated. Subsequently, all the intermediate heights are linearly interpolated from the previous ones.
One can also retrieve the response functions for a supplied atmosphere. As stated in Sect. 2.2, these functions encode information about where Stokes profiles are sensitive, to what they are sensitive, and to what degree. To the best of our knowledge, this is the first code capable of calculating response functions analytically at geometrical scale height. As an example, Fig. 2 shows the Stokes profiles (top row) in a penumbral pixel along with the corresponding dimensionless response functions to temperature T (second row), gas pressure P_{g} (third row), and density ρ (fourth row), three components of the magnetic field B_{x}, B_{y}, and B_{z} (fifth to seventh rows, respectively), and v_{LOS} (eighth row). The normalization factors employed to make the response functions dimensionless are: K, gr cm^{3}, dyn cm^{−2}, G (likewise for the other components of the magnetic field), and cm s^{−1}. It is important to emphasize that the response functions to gas pressure and density are not negligible compared with the response functions to other physical parameters. The high sensitivity to these two parameters could be surprising, as when solving the RTE in optical depth scale, the Stokes profiles are only slightly sensitive to them. The difference resides in the fact that solving the RTE in geometrical scale (as is the case here), the optical depth scale is affected by the combination of temperature, gas pressure, and density. Consequently, the thermodynamic response functions allow for a better insight into how (i.e., what combination of thermodynamic parameters) makes the optical depth scale behave as it does. For instance, far away from the line core of the two spectral lines in Fig. 2, the response functions of Stokes I to P_{g} and ρ are not zero and behave more like the response functions to the temperature and unlike the response function to B_{x}, B_{x}, B_{x}, and so on. Consequently, the gas pressure and density have an important effect on the continuum intensity by changing the geometrical height z at which the continuum level τ_{5} = 1 is formed.
Fig. 1.
Examples of three simulated Stokes profiles as RTE is solved for a snapshot of the Rempel simulation of a sunspot (Rempel 2012). Top panel: monochromatic continuum intensity map in which three pixels are highlighted in black, red, and lightblue squares as representatives of low magnetized, penumbral, and umbral areas, respectively; their spectral profiles are presented following this same color code in panels a–d for Stokes I, Q, U, and V, respectively. The physical parameters (T(z), P_{g}(z), ρ(z), v_{LOS}(z), B_{x}(z), B_{y}(z), and B_{z}(z)) that lead to these profiles are shown in panels e–k, respectively. 
Fig. 2.
Analytical dimensionless response functions for, from the second row to the last one, T(z), P_{g}(z), ρ(z), B_{x}(z), B_{y}(z), B_{z}(z), and v_{LOS}(z) for the penumbral pixel highlighted in Fig. 1 (red square) and whose profiles are shown on the first row. Columns, from left to right: Stokes I, Q, U, and V. 
3.2. Inverse solver
The inverse solver aims at inferring the vertical stratification of the various physical parameters of the atmosphere X(z), provided a set of observed Stokes profiles, 𝕀(λ)^{obs}. To test the performance of our inverse solver we take as observed Stokes profiles the simulated ones 𝕀^{syn}(λ) of the spectral lines described in Sect. 3.1, in particular the penumbral pixel (red case in Fig. 1). The number of data points is 2000 (500 wavelength points for each of the four Stokes parameters). In principle, the inference can be done for any combination of: temperature, T(z), gas pressure, P_{g}(z), density, ρ(z), the three cartesian components of the magnetic field, B_{x, y, z}(z) and the LOS velocity, v_{LOS}(z). The only restriction is that it is not possible to simultaneously invert all three thermodynamic parameters: T(z), P_{g}(z) and ρ(z). This occurs because provided two of them, the third one can be determined through the EOS. Therefore, inverting the three of them simultaneously can potentially lead to thermodynamic physical parameters that are not consistent with the EOS.
In this paper we focus mostly on the influence that the thermodynamic parameters have on an inversion code that works in the geometrical scale z instead of the usual opticaldepth τ. Because of this, most of the forthcoming discussion avoids the performance of the code in the determination of the three components of the magnetic field and LOS velocity.
As a first test we only invert one of the three thermodynamic physical parameters at a time. The other two thermodynamic parameters as well as the three components of the magnetic field and v_{LOS} are not inverted at this point. In addition, and this is the case for all the forthcoming tests but the last, we set the three components of the magnetic field and v_{LOS} to a constant in height value, meaning they do not encode any information about the height. The inverted parameters are allowed to change at each of the 256 vertical heights (same as the MHD simulation employed to produce 𝕀^{syn}(λ)) and therefore in this test we have 256 free parameters.
In Fig. 3 we depict the real (black solid lines), starting guess (thin color lines), and inferred (solid color lines) stratifications for the temperature T(z) (purple), gas pressure P_{g}(z) (red), and density ρ(z) (blue) for the three different inversions. As can be seen, the recovered parameters are in good agreement in the region where the profiles are sensitive to: z ∈ [1000, 1300] km as obtained from the heights at which z(lgτ_{5} = 0)≈1000 and z(lgτ_{5} = −3)≈1300 km (see horizontal color bars in this plot). Shadowed areas represent the lower limit of the uncertainties on the inference of each physical parameter. These error bars are tiny because of the absence of noise in the profiles together with the large number of data points considered.
Fig. 3.
Results of the inversion of three thermodynamic parameters. Three different inversions are presented, in each of which a single physical parameter is inverted: the temperature, T(z), (solidpurple line), the gas pressure, P_{g}(z), (dashedred), and the density, ρ(z), (dashdottedblue). Left panels: simulated spectral lines inverted (𝕀^{syn}(λ), in black) coming from the same penumbral pixel highlighted in Fig. 1 (red square) together with the results of the three inversions. Right panels: thermodynamic parameters (synthetic X^{syn}(z) and inverted X′(z)), from top to bottom, T(z), P_{g}(z), and ρ(z), (the initialization for each inversion is shown in thin lines for each thermodynamical parameter). Lineofsight velocity and magnetic field components are taken to be constant with height. For the sake of clarity, in each panel we only show the physical parameter inverted, as the other ones are taken as known. The colored shadows represent the uncertainties as derived by the code (see Sect. 2.3). Right panels include a horizontal line for each physical parameter (following the same color code) with ticks for, from left to right, the heights at which lgτ_{5} = 0, −1, −2, −3. 
An important drawback when working in the geometrical scale z, compared to optical scale τ, is that there are two intrinsic degeneracies involving the thermodynamic properties. (1) Any shift in height of the whole atmosphere leads to exactly the same Stokes spectra. (2) There are different combinations of the thermodynamic parameters that can lead to the same absorption coefficient. For instance, one could increase the density (increasing the number of atoms) and at the same time decrease the temperature (decreasing the population of an ionization stage) so that the absorption coefficient remains the same. To better illustrate this second point (as the first one is straightforward), Fig. 4 shows the inversion of a simulated Stokes profile of the same spectral line configuration employed until now, 𝕀^{syn}(λ), which is obtained through a given atmospheric model (black lines). The strategy is to invert the gas pressure P_{g}(z) while forcing the temperature stratification T(z) to be different (see blue and red color lines) from the one used for the synthesis. This way, if we can fit, up to a given threshold, the Stokes profiles for the different fixed temperature stratifications, it means that there exists more than one combination of thermodynamic parameters as a function of the geometrical scale z that yield the same Stokes profiles. As shown, this is exactly the case, as for the different temperatures the gas pressure and densities have changed (the former through the inversion and the latter through the EOS) so that the Stokes profiles ( and ) are approximately equally well fit. It is worth noticing that, even for the extreme temperatures provided, the fit of the Stokes profiles is quite good, in particular for the first case (red). The second one (blue) is slightly worse because, as the supplied temperature gradient is steeper than the reference one, the spectral line formation heights are comprised in a smaller number of slabs and so the solution of the RTE is not as accurate as for the simulated case, giving rise to small discrepancies in the stokes profiles, mostly close to the core of the line in any Stokes parameter. The number of free parameters in this test is the same as in the previous one: 256.
One possibility to circumvent this degeneracy is to assume hydrostatic equilibrium. Here, the momentum equation (see, e.g., Priest 1987), simplifies to dP_{g}/dz = − ρg, where g refers to the surface gravity of the Sun. Under such an assumption and together with the equation of state, the gas pressure and density, P_{g}(z) and ρ(z), can be obtained provided the temperature T(z) plus a boundary condition for the gas pressure. This can be done numerically by means of any onedimensional solver for firstorder differential equations (e.g., Runge–Kutta). By including this additional constraint, one imposes a tight restriction on the thermodynamic stratifications so that only T(z) and a boundary condition for the gas pressure remain as free parameters during the inversion. An equivalent problem can be formulated using the density ρ(z) as free parameter instead of the temperature. This assumption considers the most simplified version of the MHD equation of motion and is often made for RTE solver codes working in optical depth τ such as SIR (Ruiz Cobo & del Toro Iniesta 1992), SPINOR (Frutiger et al. 2000), and NICOLE (SocasNavarro et al. 2015).
Fig. 4.
Two different inversions of the gas pressure, P_{g}(z), for the penumbral pixel highlighted in Fig. 1 by a red square. Left column, from top to bottom: Stokes I, Q, U, and V for the simulated case (𝕀^{syn}(λ) solidblack lines) and for the two inversions ( solidred lines and dashedblue lines) with fixed and “different from the simulated case” temperatures: T_{1}(z) and T_{2}(z), respectively. Right column: thermodynamic parameters (synthetic X^{syn}(z) and inverted X′(z)) of the various models, following the same color code. Upper panel: horizontal line for each model atmosphere (following the same color code) with ticks for, from left to right, the heights at which lgτ_{5} = 0, −1, −2, −3. 
In Fig. 5 we test the uniqueness of the solution in the thermodynamic parameters using the additional restriction of hydrostatic equilibrium described above. To do so we supplied a temperature stratification T(z) as well as a value for the gas pressure in the uppermost atmospheric layer P_{g}(z_{max})=25.4 dyn cm^{−2}, with z_{max} = 2048 km, and calculate the gas pressure and density stratifications, P_{g}(z) and ρ(z) using hydrostatic equilibrium. The resulting atmosphere, which we use to obtain simulated Stokes profiles, 𝕀^{syn}(λ), are denoted by black lines in Fig. 5. Once this atmosphere and its resulting Stokes profiles are set, we perform different inversions for the temperature stratifications, T(z), under the assumption of hydrostatic equilibrium, while fixing different values for the gas pressure at the top boundary, dyn/cm^{2} and dyn/cm^{2}. Since only one physical parameter is inverted for all available heights, the number of free parameters is still 256. If 𝕀^{syn}(λ) is successfully fitted, this readily implies that the temperature T(z) is able to compensate for the different values of the boundary condition in the gas pressure. As is seen in the middle row of Fig. 5, this is exactly the case, as there are several thermodynamic stratifications that fit the spectra equally well, up to an acceptable error.
Fig. 5.
Inversion of Stokes profiles in hydrostatic equilibrium for different values of the gas pressure top boundary condition, P_{g}(z_{max} = 2048 km). Left column: Stokes parameters, from top to bottom, I, Q, U, and V, for the simulated case, 𝕀^{syn}(λ) (solidblack lines), the inversion with dyn cm^{−2} (, solid red lines) and the inversion with dyn cm^{−2} (, dash blue lines). Middle column: thermodynamic parameters used for the simulated data (X^{syn}(z)) together with the ones obtained for each inversion ( and , following the same color code). Right column: same thermodynamic parameters (, , and ), but setting a common height origin for the three cases h(τ_{5} = 1)=1000 km. In the uppercenter and right panels we include horizontal lines with ticks for, from left to right, the heights at which lgτ_{5} = 0, −1, −2, −3 for each atmosphere (following the same color code). 
We note that, in spite of the seemingly very different stratifications of the thermodynamic parameters (middle panels in Fig. 5), they agree well if shifted to share a common height origin. In Fig. 5 (rightmost panels), we perform a shift of the zscale such that the geometrical height at which the τ_{5} = 1 forms is the same in all three (black –– and blue and red ––) curves. The shifts applied in this case are Δz^{1} = 83 and Δz^{2} = −228 km, for the cases of and , respectively, as the top boundary conditions in the gas pressure. This result implies that the assumption of hydrostatic equilibrium breaks the aforementioned degeneracy regarding the thermodynamic parameters (point 2 in Sect. 3.2), and allows us to infer the correct temperature stratification T(z) (as long as P_{g} and ρ are obtained in hydrostatic equilibrium as well), save for a shift in the absolute reference for the origin of the zscale. This was expected since, as we explained above (point 1), any shift in height of the atmosphere model as a whole leads to exactly the same Stokes parameters. Setting a common origin would require taking into account additional terms in the momentum equation, in addition to the hydrostatic one, such as the magnetic pressure and tension terms (see Sect. 2.7 of Priest 1987). This cannot be done in a onedimensional approach like the one used in this paper, where each pixel (x,y) on the solar surface is treated independently, and therefore setting a common origin lies out of the scope of the present investigation. Yet, some general scenarios are discussed later.
Finally, we perform a similar inversion as the one we have just carried out but this time applied to the Stokes profiles coming from the same snapshot shown in Fig. 1. Here the physical parameters used to produce the simulated Stokes pseudovector, 𝕀^{syn}(λ), do not rely on the assumption of hydrostatic equilibrium, but rather they verify the full set of MHD equations. Therefore, the inversion, assuming that hydrostatic equilibrium applies, will provide insight into how accurate the inferences of the various heightdependent parameters are. Contrary to all previous tests, we now include in the inversion the three components of the magnetic field, B_{x, y, z}, as well as the LOS velocity v_{LOS} and we do not force them to be constant in height. The inversion is performed for the whole height scale z grid and thus there are 1280 free parameters.
We consider the same pixel as before (red square in Fig. 1). The physical parameters () together with the Stokes simulated spectra they give to (𝕀^{syn}(λ)) are shown in Fig. 6 (black lines). The results of the inversion for the temperature T(z), the LOS velocity v_{LOS}, and the three cartesian components of the magnetic field are shown in red (). We recall here that we employ two different approaches to deal with the large number of free parameters involved in the inverse problem (1280). On the one hand, the inversion is performed in a sequential manner. Thus, in a first step or cycle, each physical quantity inverted (T, v_{LOS}, and B_{x, y, z}) was modified by a constant value across all zslabs. In the following steps or cycles, the perturbation for each parameter was linearly interpolated from the calculated ones at 8, 64, and finally all 256 different heights (equally spaced in z). On the other hand, the regularization introduced during the inversion of the dampened Hessian, (highly) reduces the dimension of the problem. For instance, in the last step of the inversion, there are 256 free parameters per physical quantity, but as is evident from middle column in Fig. 6, perturbations in slabs below 900 km are ignored by the SVD method as the eigenvalues of the dampened Hessian associated to these heights are below the 10^{−4} threshold. Something similar applies to higher layers where the density is so small that Stokes profiles are highly insensitive to the physical parameters on these slabs. In addition to these, the regularization also implies a smoothing of the solution in the sense that it reduces the effective height resolution. The smaller the threshold the larger the smoothing. For the inversion shown in Fig. 6, of the 1280 possible free parameters only about 350 are really considered in the inversion after applying a SVD threshold of 10^{−4}. The result, as can be seen from the agreement between the simulated (black lines on the leftmost column in Fig. 6) and fitted (red lines in the same panels) Stokes profiles, is good.
Fig. 6.
Inversion of a penumbral pixel (highlighted in Fig. 1 with a red square). In black, the simulated Stokes parameters, 𝕀^{syn}(λ), (left column) and atmosphere model parameters inverted in geometrical scale (X^{syn}(z), middle column) and in optical depth scale (X^{syn}(lgτ_{5}), rightmost column). The results of the inversion in hydrostatic equilibrium for the third and fourth cycle (see text), ) and ), are in solid red and dashed blue, respectively. The physical parameters inverted are T(z), v_{LOS}(z), B_{x}(z), B_{y}(z), and B_{z}(z). Since the height origin is not absolute (see Sect. 3.2), each atmosphere model is forced to have h(τ_{5} = 1) = 1000 km. The atmosphere models are only shown between lg(τ_{5}) = 1.5 and lg(τ_{5}) = − 5 as the Stokes profiles, given the spectral line setup, are sensitive inside these optical depths, though the whole atmosphere is inverted. The shadowed area shows the lower limit errors for each physical parameter. 
If we focus on the physical parameters we first see that inferred parameters sometimes show large variations between slabs, mostly in LOS velocity and the magnetic field components. The reason is that since we are using the original height grid from the MHD simulation, the effect of a single slab over the overall Stokes profile appearance might be tiny, and so large excursions are possible with this inversion setup. Also, we see that the assumption of hydrostatic equilibrium leads to a general misfit of the z scale (middle row in Fig. 6), that is mostly seen as a different slope between the simulated and inferred temperature and as a compression of the inferred (red) LOS velocity and magnetic field components with respect to the simulated ones (black). We notice that, as mentioned before, there is no absolute heightscaleorigin and therefore we set the height at which τ_{5} = 1 as 1000 km. In contrast, the rightmost column of Fig. 6 shows that, in optical depth scale and in between lgτ_{5} = 0 and lgτ_{5} = −3, the physical parameters are nicely inferred. It is interesting to point out that the excellent recovery of the optical depth stratification is the result of the combined action of the two points mentioned above. First, if there is no regularization applied while calculating the inverse of the dampened Hessian, there is no convergence at all as the inverted matrix is (almost certainly for this large number of free parameters) quasisingular. Second, without the sequential approach to the solution by progressive increase of the number of free parameters, the inversion results are strongly determined by the initial guess model. In this sense, the first inversion steps or cycles with a reduced number of free parameters can be seen as a way to get a good enough guess for the final one, as they do not capture all the complexity of the profiles but they are better suited to retrieve the general atmospheric behavior. For instance, in Fig. 6, we have included in red and blue the final atmospheric model for the third and fourth cycles of the inversion, respectively (i.e., the ones with 64 and 256 perturbations per physical parameter). The results after the third cycle are already very close to the real one (black lines) but there are still some features that are slightly improved in the fourth step or “cycle” of the inversion, where the perturbation of the model atmosphere is calculated over the whole geometrical scale grid. The sequential increase of the number of free parameters must be carried out with caution when applying the inversion to real observations. This is illustrated in Fig. 7, where the differences between the fitted and simulated profiles after the third and fourth inversion cycles are displayed (following the same color code as in Fig. 6). These differences are already close to ∼10^{−3}, 10^{−4}I_{c} (typical noise level in modern spectropolarimeters) after the third cycle (solid red lines) and therefore a small number of free parameters might be kept in real applications.
Fig. 7.
Difference between the Stokes profiles of the inversion results and the inverted ones for the example considered in Fig. 6 (following the same color code). 
One may wonder as to why the model atmosphere in optical depth is so good while in geometrical height it is not. The reason is that in order to fit the Stokes profiles (𝕀^{syn}(λ)), the temperature, magnetic field components and LOS velocity in optical depth must be the same (up to a given accuracy) to that used for the synthesis. Yet, since hydrostatic equilibrium involves a very strong restriction, as it ignores all the other forces present in the MHD equation of motion, the gas pressure and density obtained in this fashion are unrealistic and so is the continuum opacity that relates geometrical and optical depth scales. Thus, in geometrical scale, the various physical parameters inverted, are inaccurately inferred. Here, it is clear that the ability to determine the various physical parameters in a threedimensional volume from a set of polarimetric observations relies on the accuracy with which the gas pressure is known or estimated.
4. Conclusions and discussion
We present a numerical code that solves the RTE for polarized light (in the Zeeman regime) in geometrical scale z in a onedimensional setup, which is made publicly available^{2}. It includes, for the first time, the analytical calculation of the Response functions in height scale z, which allows a better understanding as to how the optical depth, and as a consequence the spectral lines, are formed in z. It is also possible to solve the inverse problem, inferring all relevant physical parameters from a given observation. The code is parallelized so that large datasets can be easily tackled both in the forward and the inverse solver. Here, we have tested the performance of the inversion using a huge number of free parameters, yet this number can be easily reduced if desired.
We see that solving the RTE in geometrical scale z yields, when facing the inverse problem, to two intrinsic degeneracies: (1) the absence, when analyzing each pixel independently, of an absolute height scale origin, and (2) solving the RTE in geometrical scale leads to a degeneracy concerning the thermodynamical parameters: T(z), P_{g}(z), and ρ(z). The latter appears because there are different combinations of these physical parameters that yield the same absorption coefficient. In other words, there is no guarantee that, provided an observation, different attempts at the simultaneous inversion of two thermodynamical parameters will lead to the same geometrical scale z stratification.
One possibility to avoid this multiplicity on the solution of the inverse problem is to include an additional constraint such as hydrostatic equilibrium. Under this assumption, provided a boundary condition for the gas pressure, P_{g}(z_{max}), and either T(z) or ρ(z), the remaining stratification of the thermodynamic parameters can all be obtained without degeneracy. The inclusion of this constraint in the solution of the inverse problem guarantees the uniqueness of inferred stratifications as a function of z, except for an absolute heightscale origin. The absence of such an origin cannot be addressed in a onedimensional approach, as it requires that additional terms be considered (e.g., magnetic pressure and tension), in addition to the hydrostatic one, in the momentum equation. This is evident as the Lorentz force involves the calculation of spatial derivatives that can only be evaluated in a threedimensional approach.
Finally, we tested the effects of the assumption of hydrostatic equilibrium on the inference of heightdependent physical parameters X(z). We see that since the gas pressure obtained under this assumption is not very accurate, the inferred zdependences of the physical parameters deviate from the original ones. However, the inference of the temperature, the three components of the magnetic field vector, and LOS velocity is accurate in the optical depth scale τ.
These two open problems, namely the lack of an absolute height origin when treating each pixel individually (i.e., onedimensional approach), and the unreliability of the inferred gas pressure and density (in both the z and τ scales), can be addressed by including more terms, in addition to the hydrostatic one, in the momentum equation. This would allow a more realistic gas pressure, P_{g}(z), and density ρ(z) to be retrieved. Once this is achieved, as we demonstrate here, all other physical parameters (including the magnetic field) can be reliably determined as a function of the geometrical height z, which in turn would allow a more realistic determination of electric currents, , and to infer a proper solenoidal magnetic field ∇ ⋅ B = 0. Since the proper determination of the currents and magnetic fields relies on the accuracy of the P_{g}(z) and, at the same time, the latter depends on the currents and magnetic fields, this process will require an iterative approach. This would represent a major leap in the study of the magnetic topology in the solar atmosphere.
In the past, some attempts have been made with this aim (Solanki et al. 1993; Martinez Pillet & Vazquez 1993; Mathew et al. 2004), but it remains unclear how to proceed from observations to threedimensional physical parameters. An ambitious possibility, presented by Riethmüller et al. (2017), is to directly use stateoftheart MHD codes to iteratively retrieve a simulation run that gets closer and closer to an observation. This way, since one is already solving the whole threedimensional physical volume consistently, the final atmosphere has taken into account all the information available. However, this approach is computationally expensive, as it requires running MHD simulations and there is no guarantee that the simulation reproduces all the observed features. A different approach was presented by Puschmann et al. (2010) where they inverted an area of the solar surface under a onedimensional approach and then included additional constraints by means of the hydrodynamic and magnetic terms of the MHD equation together with the ∇ ⋅ B = 0 condition. However, in this approach the new and final threedimensional atmosphere is not fully compatible with the observations. In addition, since it relies on evaluations of vertical displacements in each pixel, for large fieldsofviews the number of free parameters becomes extremely large. A possibility to circumvent the latter was presented by Löptien et al. (2018), where they used a modification of the method of Puschmann et al. (2010), only accounting for the ∇ ⋅ B = 0 condition and using a Fourier decomposition so that the free parameters can be taken to a workable size. In this work, we present a numerical code that allows this problem to be faced in a different way, yet follows the idea presented by Puschmann et al. (2010), namely including the hydrodynamic and magnetic terms of the MHD equation and the ∇ ⋅ B = 0 condition so that the threedimensional atmosphere is horizontally consistent. Here however there are two advantages: (1) the horizontal derivatives required that the additional terms included in the equation of motion can be directly handled, and (2) the results of solving these additional equations, which lead to a different gas pressure, can be used to feed again the inversion code to retrieve a threedimensional atmosphere that would be fully consistent with the observations.
Acknowledgments
We thank the referee for many interesting suggestions that improved the paper. This work was supported by the German Deutsche Forschungsgemeinschaft, DFG project number 321818926. JMB acknowledges financial support from the Spanish Ministry of Economy and Competitiveness (MINECO) under the 2015 Severo Ochoa Program MINECO SEV20150548.
References
 Anstee, S. D., & O’Mara, B. J. 1995, MNRAS, 276, 859 [NASA ADS] [CrossRef] [Google Scholar]
 Barklem, P. S. 1998, PhD Thesis, The University of Queensland, Australia [Google Scholar]
 Barklem, P. S., & O’Mara, B. J. 1997, MNRAS, 290, 102 [Google Scholar]
 Barklem, P. S., Anstee, S. D., & O’Mara, B. J. 1998, PASA, 15, 336 [NASA ADS] [CrossRef] [Google Scholar]
 Borrero, J. M., Rempel, M., & Solanki, S. K. 2010, Astron. Nachr., 331, 567 [NASA ADS] [CrossRef] [Google Scholar]
 Borrero, J. M., Tomczyk, S., Kubo, M., et al. 2011, Sol. Phys., 273, 267 [Google Scholar]
 Borrero, J. M., Lites, B. W., Lagg, A., Rezaei, R., & Rempel, M. 2014, A&A, 572, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cao, W., Gorceix, N., Coulter, R., Coulter, A., & Goode, P. R. 2010, in Groundbased and Airborne Telescopes III, Proc. SPIE, 7733, 773330 [CrossRef] [Google Scholar]
 Collados, M., Bettonvil, F., Cavaller, L., et al. 2010, in Groundbased and Airborne Telescopes III, Proc. SPIE, 7733, 77330H [CrossRef] [Google Scholar]
 Danilovic, S., Gandorfer, A., Lagg, A., et al. 2008, A&A, 484, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Danilovic, S., Schüssler, M., & Solanki, S. K. 2010a, A&A, 509, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Danilovic, S., Schüssler, M., & Solanki, S. K. 2010b, A&A, 513, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de la Cruz Rodríguez, J., Leenaarts, J., Danilovic, S., & Uitenbroek, H. 2019, A&A, 623, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 del Toro Iniesta, J. C., & Ruiz Cobo, B. 2016, Liv. Rev. Sol. Phys., 13, 4 [Google Scholar]
 del Toro Iniesta, J. C., Ruiz Cobo, B., Bellot Rubio, L. R., & Collados, M. 1995, A&A, 294, 855 [NASA ADS] [Google Scholar]
 Freytag, B., Steffen, M., Ludwig, H.G., et al. 2012, J. Comput. Phys., 231, 919 [Google Scholar]
 Frutiger, C., Solanki, S. K., Fligge, M., & Bruls, J. H. M. J. 2000, A&A, 358, 1109 [NASA ADS] [Google Scholar]
 GonzálezMorales, P. A., Khomenko, E., Downes, T. P., & de Vicente, A. 2018, A&A, 615, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kosugi, T., Matsuzaki, K., Sakao, T., et al. 2007, Sol. Phys., 243, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Landi Degl’Innocenti, E. 1976, A&AS, 25, 379 [NASA ADS] [Google Scholar]
 Landi Degl’Innocenti, E., & Landi Degl’Innocenti, M. 1972, Sol. Phys., 27, 319 [NASA ADS] [CrossRef] [Google Scholar]
 Landi Degl’Innocenti, E., & Landi Degl’Innocenti, M. 1985, Sol. Phys., 97, 239 [NASA ADS] [CrossRef] [Google Scholar]
 Landi Degl’Innocenti, E., & Landolfi, M. 2004, in Polarization in Spectral Lines, Astrophys. Space Sci. Lib., 307 [Google Scholar]
 Levenberg, K. 1944, Quart. Appl. Math., 2, 164 [Google Scholar]
 Lites, B. W., Rempel, M., Borrero, J. M., & Danilovic, S. 2017, ApJ, 835, 14 [Google Scholar]
 Löptien, B., Lagg, A., van Noort, M., & Solanki, S. K. 2018, A&A, 619, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marquardt, D. W. 1963, J. Soc. Ind. Appl. Math., 11, 431 [Google Scholar]
 Martinez Pillet, V., & Vazquez, M. 1993, A&A, 270, 494 [NASA ADS] [CrossRef] [Google Scholar]
 Mathew, S. K., Solanki, S. K., Lagg, A., et al. 2004, A&A, 422, 693 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Milić, I., & van Noort, M. 2018, A&A, 617, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Priest, E. R. 1987, Solar Magnetohydrodynamics (Dordrecht: D. Reidel) [Google Scholar]
 Puschmann, K. G., Ruiz Cobo, B., & Martínez Pillet, V. 2010, ApJ, 720, 1417 [NASA ADS] [CrossRef] [Google Scholar]
 Rachkovsky, D. N. 1967, Izvestiya Ordena Trudovogo Krasnogo Znameni Krymskoj Astrofizicheskoj Observatorii, 37, 56 [NASA ADS] [Google Scholar]
 Rempel, M. 2012, ApJ, 750, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Rempel, M., Schüssler, M., Cameron, R. H., & Knölker, M. 2009, Science, 325, 171 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Riethmüller, T. L., Solanki, S. K., Barthol, P., et al. 2017, ApJS, 229, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Rimmele, T. R., Wagner, J., Keil, S., et al. 2010, in Groundbased and Airborne Telescopes III, Proc. SPIE, 7733, 77330G [CrossRef] [Google Scholar]
 Ruiz Cobo, B., & del Toro Iniesta, J. C. 1992, ApJ, 398, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Scharmer, G. B., Bjelksjo, K., Korhonen, T. K., Lindberg, B., & Petterson, B. 2003, in Innovative Telescopes and Instrumentation for Solar Astrophysics, eds. S. L. Keil, & S. V. Avakyan, Proc. SPIE, 4853, 341 [NASA ADS] [CrossRef] [Google Scholar]
 Scherrer, P. H., Schou, J., Bush, R. I., et al. 2012, Sol. Phys., 275, 207 [Google Scholar]
 Schmidt, W., von der Lühe, O., Volkmer, R., et al. 2012, Astron. Nachr., 333, 796 [NASA ADS] [CrossRef] [Google Scholar]
 SocasNavarro, H., de la Cruz Rodríguez, J., Asensio Ramos, A., Trujillo Bueno, J., & Ruiz Cobo, B. 2015, A&A, 577, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Solanki, S. K., Walther, U., & Livingston, W. 1993, A&A, 277, 639 [NASA ADS] [Google Scholar]
 Unno, W. 1956, PASJ, 8, 108 [NASA ADS] [Google Scholar]
 Unsold, A. 1955, Physik der Sternatmospharen, MIT besonderer Berucksichtigung der Sonne (Berlin: Springer) [Google Scholar]
 Wittmann, A. 1974, Sol. Phys., 35, 11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Derivative of the number of electron density with respect to temperature at constant gas density
We focus on the derivative of the density of electrons with respect to the temperature at constant gas pressure. The number of electron density is given by the sum over the whole set of atomic elements (index i) and for the ionization levels (index j):
where n_{i, j} gives the number of atoms per cm^{3} of the species i at the ionization stage j. In the current implementation of the code, N_{ion} = 3 and C_{i, j} encodes the additional electrons freed (absorbed) by each species at each ionization stage. For hydrogen, j = 1 refers to the negative hydrogen ion so C_{1, 1} = −1, j = 2 is the neutral hydrogen (C_{1, 2} = 0), and j = 3 is the ionized hydrogen (C_{1, 3} = 1). For the rest of elements, j = 1 is the neutral atom (C_{i ≠ 1, 1} = 0) and j = 2 and j = 3 refer to the once and twice ionized atoms (C_{i ≠ 1, 2} = C_{i ≠ 1, 3} = 1). Therefore, the derivative remains as:
Now n_{i, j} is given by the product of the density of i atoms multiplied by the fraction of these atoms in the ionization state j:
where n_{1} is the density of hydrogen atoms and A_{i} the abundance of species i, which is defined in the usual astronomical scale: A_{i} = lg(n_{i}/n_{H})+12. We now refer to f_{i, j} as the fraction of atoms of the species i at the ionization state j with respect to the total number of the same species i:
where we assume that there is no element that is ionized more than three times. Since we are assuming LTE, the various ionization states are given by the Saha equation, so it can be seen that the population of each ionization state is:
where
Here, u_{i, j} refers to the partition function of the species i at the ionization state j, , ϵ_{i, j} is the ionization potential, and h, m_{e}, and k_{B} are the Planck, the electron mass, and the Boltzmann constants, respectively.
Now, one can substitute Eqs. (A.2) and (A.3) into Eq. (A.1) and get
where we can substitute n_{1}, provided that our equation of state is the ideal gas equation, so that its derivative is
where s′=∑_{i}10^{Ai − 12}m_{i}. We focus now on the last derivative appearing in Eq. (A.6) and we define so that
From Eq. (A.4), it is evident that both derivatives on the righthand side of the Eq. (A.8) depend on the derivative we want to calculate (). Let us write and as a sum of a part that depends on and another that does not depend on it:
where the definitions for 𝒟_{i}, ℱ_{i}, 𝒢_{i, j}, and ℋ_{i, j} are detailed below. If we introduce Eq. (A.9) into Eq. (A.8) and then this latter together with Eq. (A.7) into (A.6). These, we obtain
We can now group terms with and get the analytical expression for it:
This allows the calculation of the derivative of the number of electron density with respect to the temperature at constant density. We just need to specify the expressions for: 𝒟_{i}, ℱ_{i}, 𝒢_{i, j}, and ℋ_{i, j}. If we write
then
where the various terms for different ionization stages are gathered in Table A.1, and
Actual mathematical expressions for 𝒜_{i, j}, ℬ_{i, j}, 𝒞_{i, j}, 𝒢_{i, j}, and ℋ_{i, j} terms for the various ionization stages.
Appendix B: Derivative of the number of electron density with respect to gas pressure at constant temperature
For this derivation, we can reuse the first equations appearing in Appendix A up to Eq. (A.6) substituting with and so
The derivative can be calculated making use of the equation of state:
and for the second derivative of the righthand side of Eq. (B.1), we can follow a similar strategy as in Appendix A:
where we now express the derivatives of the righthand side as function of a term that depends on plus another one that does not:
Subsequently, substituting Eq. (B.4) in Eq. (B.3) and the latter together with Eq. (B.2) in Eq. (B.1) we get
Now we gather together the terms with and solve it, obtaining
where it is to be noted that 𝒢_{i, j} and 𝒟_{i} do not appear, as they are null. We can see this point and retrieve the expressions for ℋ_{i, j} and ℱ_{i} if we write
where 𝒜_{i, j}, ℬ_{i, j}, and 𝒞_{i, j} are detailed in Table B.1. Then, taking the derivative of m_{i, j} with respect to gas pressure at constant temperature and bearing in mind Eq. (B.4) we get
where the various terms for different ionization stages are listed in Table B.1. The terms 𝒜_{i, j}, ℬ_{i, j}, and 𝒞_{i, j} for the various ionization stages (j) are either constants or functions of α_{i, 1} and/or α_{i, 2}. These terms, as can be seen in Eq. (A.5), only depend on the temperature, so their derivatives with respect to gas pressure at constant temperature vanish. In a similar way, it can be shown that by taking the derivative of d_{i} with respect to gas pressure at constant temperature, one leads to
where these derivatives vanish as α_{i, 1} and α_{i, 2} only depend on the temperature.
Actual mathematical expressions for 𝒜_{i, j}, ℬ_{i, j}, 𝒞_{i, j}, 𝒢_{i, j}, and ℋ_{i, j} terms for the various ionization stages considered.
Appendix C: Derivative of the number of electron density with respect to density at constant temperature
This derivation is very similar to the previous one (Appendix B), where one has to use instead of and taking into account that the hydrogen density number derivative is now
Thus, following the same steps as in Appendix B, one reaches a similar expression for :
where the various terms appearing in the expression can be found in Appendix B.
All Tables
Actual mathematical expressions for 𝒜_{i, j}, ℬ_{i, j}, 𝒞_{i, j}, 𝒢_{i, j}, and ℋ_{i, j} terms for the various ionization stages.
Actual mathematical expressions for 𝒜_{i, j}, ℬ_{i, j}, 𝒞_{i, j}, 𝒢_{i, j}, and ℋ_{i, j} terms for the various ionization stages considered.
All Figures
Fig. 1.
Examples of three simulated Stokes profiles as RTE is solved for a snapshot of the Rempel simulation of a sunspot (Rempel 2012). Top panel: monochromatic continuum intensity map in which three pixels are highlighted in black, red, and lightblue squares as representatives of low magnetized, penumbral, and umbral areas, respectively; their spectral profiles are presented following this same color code in panels a–d for Stokes I, Q, U, and V, respectively. The physical parameters (T(z), P_{g}(z), ρ(z), v_{LOS}(z), B_{x}(z), B_{y}(z), and B_{z}(z)) that lead to these profiles are shown in panels e–k, respectively. 

In the text 
Fig. 2.
Analytical dimensionless response functions for, from the second row to the last one, T(z), P_{g}(z), ρ(z), B_{x}(z), B_{y}(z), B_{z}(z), and v_{LOS}(z) for the penumbral pixel highlighted in Fig. 1 (red square) and whose profiles are shown on the first row. Columns, from left to right: Stokes I, Q, U, and V. 

In the text 
Fig. 3.
Results of the inversion of three thermodynamic parameters. Three different inversions are presented, in each of which a single physical parameter is inverted: the temperature, T(z), (solidpurple line), the gas pressure, P_{g}(z), (dashedred), and the density, ρ(z), (dashdottedblue). Left panels: simulated spectral lines inverted (𝕀^{syn}(λ), in black) coming from the same penumbral pixel highlighted in Fig. 1 (red square) together with the results of the three inversions. Right panels: thermodynamic parameters (synthetic X^{syn}(z) and inverted X′(z)), from top to bottom, T(z), P_{g}(z), and ρ(z), (the initialization for each inversion is shown in thin lines for each thermodynamical parameter). Lineofsight velocity and magnetic field components are taken to be constant with height. For the sake of clarity, in each panel we only show the physical parameter inverted, as the other ones are taken as known. The colored shadows represent the uncertainties as derived by the code (see Sect. 2.3). Right panels include a horizontal line for each physical parameter (following the same color code) with ticks for, from left to right, the heights at which lgτ_{5} = 0, −1, −2, −3. 

In the text 
Fig. 4.
Two different inversions of the gas pressure, P_{g}(z), for the penumbral pixel highlighted in Fig. 1 by a red square. Left column, from top to bottom: Stokes I, Q, U, and V for the simulated case (𝕀^{syn}(λ) solidblack lines) and for the two inversions ( solidred lines and dashedblue lines) with fixed and “different from the simulated case” temperatures: T_{1}(z) and T_{2}(z), respectively. Right column: thermodynamic parameters (synthetic X^{syn}(z) and inverted X′(z)) of the various models, following the same color code. Upper panel: horizontal line for each model atmosphere (following the same color code) with ticks for, from left to right, the heights at which lgτ_{5} = 0, −1, −2, −3. 

In the text 
Fig. 5.
Inversion of Stokes profiles in hydrostatic equilibrium for different values of the gas pressure top boundary condition, P_{g}(z_{max} = 2048 km). Left column: Stokes parameters, from top to bottom, I, Q, U, and V, for the simulated case, 𝕀^{syn}(λ) (solidblack lines), the inversion with dyn cm^{−2} (, solid red lines) and the inversion with dyn cm^{−2} (, dash blue lines). Middle column: thermodynamic parameters used for the simulated data (X^{syn}(z)) together with the ones obtained for each inversion ( and , following the same color code). Right column: same thermodynamic parameters (, , and ), but setting a common height origin for the three cases h(τ_{5} = 1)=1000 km. In the uppercenter and right panels we include horizontal lines with ticks for, from left to right, the heights at which lgτ_{5} = 0, −1, −2, −3 for each atmosphere (following the same color code). 

In the text 
Fig. 6.
Inversion of a penumbral pixel (highlighted in Fig. 1 with a red square). In black, the simulated Stokes parameters, 𝕀^{syn}(λ), (left column) and atmosphere model parameters inverted in geometrical scale (X^{syn}(z), middle column) and in optical depth scale (X^{syn}(lgτ_{5}), rightmost column). The results of the inversion in hydrostatic equilibrium for the third and fourth cycle (see text), ) and ), are in solid red and dashed blue, respectively. The physical parameters inverted are T(z), v_{LOS}(z), B_{x}(z), B_{y}(z), and B_{z}(z). Since the height origin is not absolute (see Sect. 3.2), each atmosphere model is forced to have h(τ_{5} = 1) = 1000 km. The atmosphere models are only shown between lg(τ_{5}) = 1.5 and lg(τ_{5}) = − 5 as the Stokes profiles, given the spectral line setup, are sensitive inside these optical depths, though the whole atmosphere is inverted. The shadowed area shows the lower limit errors for each physical parameter. 

In the text 
Fig. 7.
Difference between the Stokes profiles of the inversion results and the inverted ones for the example considered in Fig. 6 (following the same color code). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.