© ESO, 2017

1. Introduction

Response functions of spectral lines (Mein 1971; Beckers & Milkey 1975; Landi Degl’Innocenti & Landi Degl’Innocenti 1977) describe the sensitivity of an emergent Stokes spectrum to the perturbations of physical parameters as a function of depth in the atmosphere (for a detailed study of the possible diagnostics see Ruiz Cobo & del Toro Iniesta 1994). Although they are useful in forward modeling, in order to propose and analyze various diagnostics (e.g., Uitenbroek 2006) their main strength and application is in derivative-based spectral line fitting (in solar physics widely known as “inversion”; see, e.g., Ruiz Cobo & del Toro Iniesta 1992). Given the perturbation of an atmospheric parameter (temperature, velocity, magnetic field, etc.), the computation of the response function in practice reduces to the computation of the responses (derivatives with respect to atmospheric parameters) of opacity and emissivity and then propagating these responses using the formal solution of radiative transfer equation. It is intuitively clear that finding the responses of opacity and emissivity requires the responses of the atomic and molecular number densities as well as responses of the populations of individual levels.

In the approximation of local thermodynamical equilibrium (LTE), computing responses of level populations is relatively straightforward (e.g., Ruiz Cobo & del Toro Iniesta 1992; del Toro Iniesta 2003) as all the relevant number densities depend only on local quantities. In addition, level populations and their response can be computed analytically, using the Saha-Boltzmann equation. LTE is a good approximation for the spectral lines formed deep in the photosphere where transitions between the levels are dominated by collisions. However, for lines that are formed in the upper photosphere and the chromosphere (Hα, Ca II infrared triplet, Ca II H&K, Mg h&k lines, to name just a few), this approximation is far from valid as the transition rates are dominated by the radiation field. This introduces both spatial and nonlinear coupling between the level populations, resulting in a NLTE problem. The solution of such problems is extensively discussed in the literature and probably the best reference is a classical monograph by Mihalas (1978) (or a recent re-issue by Hubeny & Mihalas 2014). For polarized radiation, however, there are wealth of unsolved and interesting problems; a detailed and comprehensive introduction is given in Landi Degl’Innocenti & Landolfi (2004).

The response of a NLTE problem to atmospheric changes, however, has not been studied in great detail. The first attempt, by Socas-Navarro et al. (1998), was based on the fixed departure coefficients (FDC) approximation, which essentially computes the NLTE response functions by scaling the LTE ones with the departure coefficients (the ratio of the level population in NLTE to that in LTE). In that paper, the authors find it to be a rather crude approximation and instead resort to computing the response functions numerically (i.e., using finite differences). In the papers by Socas-Navarro & Uitenbroek (2004) and Uitenbroek (2006), response functions are computed numerically. Also, in the NLTE inversion code NICOLE (Socas-Navarro et al. 2015), responses – and thus the χ² derivatives to model parameters – are computed in a numerical way.

In this paper we demonstrate how the response functions of the level populations to atmospheric quantities can be computed in an analytical way, at a significantly reduced computational cost compared to finite differencing, from which it is straightforward to compute the response functions of the opacity and emissivity. Given the perturbations of opacity and emissivity, we show how to compute very precisely the resulting perturbations of the intensity throughout the atmosphere. In Sect. 2 we outline the method for analytical computation of NLTE functions and show the explicit derivation for a case of pure line transfer. In Sect. 3 we compare analytically computed response functions with those computed using finite differences and discuss the discrepancies. Finally, we conclude by discussing possible applications of our approach and future steps.

2. Response functions

In the following, we will assume that the emergent radiation is formed in a plane-parallel atmosphere represented by a discrete set of points in a single dimension (1D) at which the values of atmospheric parameters are given. The relevant atmospheric parameters for the formation of spectral lines are temperature, pressure, line-of-sight velocity, microturbulent velocity, and the magnetic field. Given the values of these parameters, we can obtain the opacity and the emissivity of the medium at each spatial and spectral point, and in each direction. These then enter the radiative transfer equation (RTE), which in the 1D axisymmetric, time-independent, polarization-free case reads $\mathit{\mu }\frac{\mathrm{d}\mathit{I}\mathrm{\left(}\mathit{z,\mu ,\lambda }\mathrm{\right)}}{\mathrm{d}\mathit{z}}\mathrm{=}\mathrm{-}\mathit{\chi }\mathrm{\left(}\mathit{z,\mu ,\lambda }\mathrm{\right)}\mathit{I}\mathrm{\left(}\mathit{z,\mu ,\lambda }\mathrm{\right)}\mathrm{+}\mathit{\eta }\mathrm{\left(}\mathit{z,\mu ,\lambda }\mathrm{\right)}\mathit{,}$(1)where z is the geometrical coordinate along the surface normal μ = cosθ, where θ is the angle with respect to the surface normal, and χ and η are the unpolarized (i.e., scalar) opacity and emissivity, respectively. Given the boundary conditions, the solution of the radiative transfer equation then yields the emergent intensity spectrum, which is to be compared with or fitted to the observed one. We note that the emergent intensity depends on the opacity and emissivity throughout the whole atmosphere (i.e., the solution is nonlocal). This nonlocallity is explicitly seen from the formal solution of radiative transfer equation ${\mathit{I}}_{\mathrm{0}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathit{S}\mathrm{\left(}{\mathit{\tau }}_{\mathit{\lambda }}\mathit{,\mu ,\lambda }\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathit{\lambda }}}\mathrm{d}{\mathit{\tau }}_{\mathit{\lambda }}\mathit{/}\mathit{\mu ,}$(2)where $\mathit{S}\mathrm{=}\frac{\mathit{\eta }\mathrm{\left(}{\mathit{\tau }}_{\mathit{\lambda }}\mathit{,\mu ,\lambda }\mathrm{\right)}}{\mathit{\chi }\mathrm{\left(}{\mathit{\tau }}_{\mathit{\lambda }}\mathit{,\mu ,\lambda }\mathrm{\right)}}$ is the source function and τ_λ is monochromatic optical depth defined by dτ_λ = −χ(z,μ,λ)dz.

We define the response function of the emergent intensity to the change of a given atmospheric quantity q at depth point k as ${\mathit{R}}_{{\mathit{q}}_{\mathit{k}}}\mathrm{\equiv }\frac{\mathit{\partial }{\mathit{I}}_{\mathrm{0}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\mathit{,}$(3)which is readily expressed in terms of the derivatives of the emergent intensity to opacity and emissivity, $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{I}}_{\mathrm{0}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}& \mathrm{=}& \sum _{\mathit{l}}^{{\mathit{N}}_{\mathrm{D}}}\frac{\mathit{\partial }{\mathit{I}}_{\mathrm{0}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{\chi }}_{\mathit{l}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}\frac{\mathit{\partial }{\mathit{\chi }}_{\mathit{l}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\\ & & \end{array}$(4)where N_D is the number of depth points in the atmosphere. The summation over l, the numerical equivalent of integration over the entire atmosphere, stems from the observation that a perturbation q_k in point k in general changes the opacity and emissivity not only in point k, but in the whole atmosphere. In general, $\frac{\mathit{\partial }{\mathit{I}}_{\mathrm{0}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{\chi }}_{\mathit{l}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}$ and $\frac{\mathit{\partial }{\mathit{I}}_{\mathrm{0}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{\eta }}_{\mathit{l}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}$ depend on the exact scheme used to solve the RTE, and we discuss this computation in more detail in Appendix A. For the time being, however, we assume that we are able to compute $\frac{\mathit{\partial }{\mathit{I}}_{\mathrm{0}}}{\mathit{\partial }{\mathit{\chi }}_{\mathit{l}}}$ and $\frac{\mathit{\partial }{\mathit{I}}_{\mathrm{0}}}{\mathit{\partial }{\mathit{\eta }}_{\mathit{l}}}$ in some way (we omit angle and wavelength dependence, which should be clear from the context).

To evaluate Eq. (4), we require the derivatives of the opacity and the emissivity everywhere with respect to quantity q_k. For the sake of simplicity, we restrict the discussion here to line processes only, and relax this assumption later. In the case of pure line transfer, we have* $\begin{array}{ccc}{\mathit{\chi }}_{\mathit{l}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}\mathrm{=}\frac{\mathit{h\nu }}{\mathrm{4}\mathit{\pi }}\mathrm{(}{\mathit{n}}_{\mathit{l,j}}{\mathit{B}}_{\mathit{ji}}\mathrm{-}{\mathit{n}}_{\mathit{l,i}}{\mathit{B}}_{\mathit{ij}}\mathrm{)}{\mathit{\phi }}_{\mathit{l}}^{\mathit{ij}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}\mathit{,}& & \end{array}$(5)where i and j are the indices of the upper and lower levels of each transition; n_l,i denotes the number density of atoms in state i at spatial point l; B_ji and B_ij are Einstein coefficients of absorption and stimulated emission, respectively; and ${\mathit{\phi }}_{\mathit{l}}^{\mathit{ij}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}$ is the line absorption profile. Similarly, $\begin{array}{ccc}{\mathit{\eta }}_{\mathit{l}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}\mathrm{=}\frac{\mathit{h\nu }}{\mathrm{4}\mathit{\pi }}{\mathit{n}}_{\mathit{l,i}}{\mathit{A}}_{\mathit{ij}}{\mathit{\phi }}_{\mathit{l}}^{\mathit{ij}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}\mathit{,}& & \end{array}$(6)where A_ij is the Einstein coefficient of spontaneous emission. Here we have assumed complete frequency redistribution (CRD), which is a good approximation for most spectral lines of interest, with the exception of very strong lines, that are formed very high in the atmosphere. Even for those lines, however, CRD usually describes the line core very well.

The derivative of the opacity with respect to q_k is then $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{\chi }}_{\mathit{l}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}& \mathrm{=}& \frac{\mathit{h\nu }}{\mathrm{4}\mathit{\pi }}\mathrm{(}{\mathit{n}}_{\mathit{l,j}}{\mathit{B}}_{\mathit{ji}}\mathrm{-}{\mathit{n}}_{\mathit{l,i}}{\mathit{B}}_{\mathit{ij}}\mathrm{)}\frac{\mathit{\partial }{\mathit{\phi }}_{\mathit{l}}^{\mathit{ij}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\\ & & \end{array}$(7)A similar expression can be written for the derivative of the emissivity. We thus see that the derivative (response) of the opacity/emissivity can be expressed through the derivatives of the level populations and the derivative of local absorption profile. The latter is determined by the temperature, microturbulent velocity, and the number density of the collisional partners (neutral hydrogen, electrons, etc.).

Although in the general problem, we can find the NLTE responses of electron density simultaneously with all other atomic populations; if we assume that the electron density response can be approximated by the LTE response, the calculation of the derivative of the absorption profile becomes strictly local, i.e.,$\frac{\mathit{\partial }{\mathit{\phi }}_{\mathit{l}}^{\mathit{ij}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\mathrm{=}{\mathit{\delta }}_{\mathit{lk}}\frac{\mathit{\partial }{\mathit{\phi }}_{\mathit{k}}^{\mathit{ij}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\mathit{,}$(8)so that the only remaining quantities to be computed are the derivatives of the level populations. Because the approximation of LTE does not hold, this is the most complicated part of the computation.

2.1. Responses of populations in NLTE

Where the approximation of LTE is accurate, the level populations are determined only by local quantities, and can be obtained with relative ease. However, if the density is sufficiently low, the approximation of LTE is not valid, and the level populations are determined both by the local temperature and by the nonlocal radiation field. Under such non-LTE (NLTE) conditions, the level populations are governed by the population evolution equation which, for the population of level i at depth point l, reads $\frac{\mathrm{d}{\mathit{n}}_{\mathit{l,i}}}{\mathrm{d}\mathit{t}}\mathrm{=}\sum _{\mathit{j}}\mathrm{(}{\mathit{n}}_{\mathit{l,j}}{\mathit{T}}_{\mathit{l,ji}}\mathrm{-}{\mathit{n}}_{\mathit{l,i}}{\mathit{T}}_{\mathit{l,ij}}\mathrm{)}\mathit{.}$(9)This is a nonlinear and nonlocal set of coupled differential equations, which is notoriously difficult to solve. However, if the conditions can be assumed stationary, $\frac{\mathrm{d}{\mathit{n}}_{\mathit{l,i}}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{0}$, then we can solve the considerably friendlier statistical equilibrium equation $\sum _{\mathit{j}}\mathrm{(}{\mathit{n}}_{\mathit{l,j}}{\mathit{T}}_{\mathit{l,ji}}\mathrm{-}{\mathit{n}}_{\mathit{l,i}}{\mathit{T}}_{\mathit{l,ij}}\mathrm{)}\mathrm{=}\mathrm{0}$(10)instead. Here T_ij is total rate of transitions from level i to level j. Generally, T_ij = C_ij + R_ij, where C stands for collisional and R for radiative transitions.

Deep in the photosphere, the density is so high that collisional rates dominate over the radiative ones and thus the approximation of LTE may be used. However, the situation is reversed in the chromosphere, with radiative transitions dominating over collisions. This means that the populations of the energy levels are no longer determined only by the temperature, but also by the radiation field, and, since the radiation field in turn depends on the emissivity and opacity, the radiation field thus effectively depends on itself. Many spectral lines are very optically thin in this region and this effect can be safely neglected, but for some spectral lines that are sufficiently strong to be still optically thick there, the populations must be calculated using Eq. (10), self-consistently with the radiation field using Eq. (1).

At least from the theoretical point of view, this problem can be considered “solved” mostly thanks to the application of the “Accelerated Lambda Iteration” methods (for an insightful review see Hubeny 2003). Here, we want to go one step further and compute not only the populations, but also their responses (derivatives). For brevity, we use the notation ${\mathrm{ℛ}}_{\mathit{lik}}\mathrm{\equiv }\frac{\mathit{\partial }{\mathit{n}}_{\mathit{l,i}}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}$(11)and refer to it as “response functions of level populations”. To compute them, we start by differentiating Eq. (10)with respect to q_k: $\sum _{\mathit{j}}\left({\mathrm{ℛ}}_{\mathit{ljk}}{\mathit{T}}_{\mathit{l,ji}}\mathrm{-}{\mathrm{ℛ}}_{\mathit{lik}}{\mathit{T}}_{\mathit{l,ij}}\mathrm{+}{\mathit{n}}_{\mathit{lj}}\frac{\mathit{\partial }{\mathit{T}}_{\mathit{l,ji}}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\mathrm{-}{\mathit{n}}_{\mathit{li}}\frac{\mathit{\partial }{\mathit{T}}_{\mathit{l,ji}}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\right)\mathrm{=}\mathrm{0.}$(12)The system in Eq. (10)is undetermined, so the statistical equilibrium equation for one of the levels is usually replaced with $\sum _{\mathit{j}}{\mathit{n}}_{\mathit{lj}}\mathrm{=}{\mathit{N}}_{\mathit{l}}\mathit{,}$(13)which, after taking the derivative becomes $\sum _{\mathit{j}}{\mathrm{ℛ}}_{\mathit{ljk}}\mathrm{=}\frac{\mathit{\partial }{\mathit{N}}_{\mathit{l}}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\mathrm{·}$(14)Where N_l is the total number density of species in question. For simplicity’s sake we will assume that this derivative is strictly local and that it can be computed directly from the equations of chemical equilibrium. To solve Eq. (12)for ℛ_lik, we need the derivative of the rates: $\frac{\mathit{\partial }{\mathit{T}}_{\mathit{l,ji}}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\mathrm{=}\frac{\mathit{\partial }{\mathit{C}}_{\mathit{l,ji}}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\mathrm{+}\frac{\mathit{\partial }{\mathit{R}}_{\mathit{l,ji}}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\mathrm{·}$(15)If we neglect the radiative rates in both Eqs. (10)and 12 and solve for populations and population responses, we end up with LTE values, however, in this specific case, it is actually easier to find level populations, as well as their responses, directly from the Saha-Boltzmann equations (see, e.g., Ruiz Cobo & del Toro Iniesta 1992).

If we account for the radiative rates in Eq. (10), but keep only the collisional rates (and their derivatives) in Eq. (12), we arrive at the FDC (fixed departure coefficients) approximation of Socas-Navarro et al. (1998), which uses the approximation ${\mathrm{ℛ}}_{\mathit{lik}}\mathrm{=}{\mathit{b}}_{\mathit{l,i}}{\mathrm{ℛ}}_{\mathit{lik}}^{\mathrm{\ast }}$(16)for the population responses, where ${\mathrm{ℛ}}_{\mathit{lik}}^{\mathrm{\ast }}$ denotes the LTE population response function and b_l,i is the NLTE departure coefficient of each level. Since the dependence of the departure coefficient on q_k is not taken into account, the nonlocal spatial coupling is not present, and the responses obtained in this way are therefore strictly local. While this makes the computations much easier, it does not account for the most significant NLTE effect, that of the nonlocality, which is usually most pronounced in the upper more dilute layers of the atmosphere, and becomes more important the finer the spatial grid is. This is probably why Socas-Navarro et al. (1998) find poor agreement between the FDC approximation and numerically computed responses in the upper layers of their atmospheres.

Clearly, a fully consistent approach needs to account not only for the correct radiative rates, but also for the response of those radiative rates.

2.2. Derivative of radiative rates

The radiative rates we are interested in (bound-bound and bound-free) are functionals of the specific intensity, and have the form ${\mathit{R}}_{\mathit{l,ji}}\mathrm{\propto }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathrm{\int }}_{-1}^{\mathrm{1}}{\mathrm{\Psi }}_{\mathit{l,ij}}{\mathit{I}}_{\mathit{l}}\mathit{\mu ,\lambda }\hspace{0.17em}\mathrm{d}\mathit{\mu }\mathrm{d}\mathit{\lambda ,}$(17)where Ψ_l,ij is some weighting function (e.g., the profile function in the case of pure line transfer). The derivative of the radiative rates then becomes $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{R}}_{\mathit{l,ji}}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}& \mathrm{\propto }& \int \begin{array}{c}\mathrm{\infty }\\ \mathrm{0}\end{array}{\mathrm{\int }}_{-1}^{\mathrm{1}}\left[\frac{\mathit{\partial }{\mathrm{\Psi }}_{\mathit{l,ij}}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}{\mathit{I}}_{\mathit{l}}\mathit{\mu ,\lambda }\mathrm{+}{\mathrm{\Psi }}_{\mathit{l,ij}}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{l}}\mathit{\mu ,\lambda }}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\right]\\ & & \mathrm{\times }\mathrm{d}\mathit{\mu }\mathrm{d}\mathit{\lambda }\mathit{.}\end{array}$(18)Further expanding the derivative of the specific intensity, in a similar fashion to Eq. (4), we get, omitting dependencies, which should be clear from the context: $\frac{\mathit{\partial }{\mathit{I}}_{\mathit{l}}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\mathrm{=}\sum _{{\mathit{l}}^{\mathrm{\prime }}}\left[\frac{\mathit{\partial }{\mathit{I}}_{\mathit{l}}}{\mathit{\partial }{\mathit{\chi }}_{{\mathit{l}}^{\mathrm{\prime }}}}\frac{\mathit{\partial }{\mathit{\chi }}_{{\mathit{l}}^{\mathrm{\prime }}}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\mathrm{+}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{l}}}{\mathit{\partial }{\mathit{\eta }}_{{\mathit{l}}^{\mathrm{\prime }}}}\frac{\mathit{\partial }{\mathit{\eta }}_{{\mathit{l}}^{\mathrm{\prime }}}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\right]\mathrm{·}$(19)Opacity and emissivity depend on all the populations of the relevant species, but also explicitly on the atmospheric parameters, $\frac{\mathit{\partial }{\mathit{\chi }}_{{\mathit{l}}^{\mathrm{\prime }}}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\mathrm{=}{\frac{\mathit{\partial }{\mathit{\chi }}_{{\mathit{l}}^{\mathrm{\prime }}}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}}_{\mathrm{explicit}}\mathrm{+}\sum _{\mathit{i}}\frac{\mathit{\partial }{\mathit{\chi }}_{{\mathit{l}}^{\mathrm{\prime }}}}{\mathit{\partial }{\mathit{n}}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{,i}}}\frac{\mathit{\partial }{\mathit{n}}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{,i}}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\mathit{,}$(20)where the summation over i takes into account all the relevant energy levels of all the species. A similar equation to Eq. (20)can be written for emissivity. Starting from Eq. (20)and substituting back all the way to Eq. (12), we end up with a linear system which couples all points in the atmosphere (because of the radiative transfer) and all the atomic level populations (because of their influence on opacity and emissivity). Below we give an explicit derivation for the case of pure line transfer and in Sect. 3.2 we treat the atomic model which also involves radiative bound-free and free-bound rates.

For pure line transfer, the radiative transition rates from bound level i to bound level j, have the form ${\mathit{R}}_{\mathit{l,ij}}\mathrm{=}{\mathit{A}}_{\mathit{ij}}\mathrm{+}{\mathit{B}}_{\mathit{ij}}{\mathit{J}}_{\mathit{l,ij}}\mathit{,}$(21)where A_ij ≡ 0 when i<j and J_l,ij is the scattering integral for the transition i → j at depth point l: ${\mathit{J}}_{\mathit{l,ij}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}{\mathit{\phi }}_{\mathit{l,ij}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}{\mathit{r}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{d}\mathit{\lambda }{\mathrm{\int }}_{-1}^{\mathrm{1}}{\mathit{I}}_{\mathit{l}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\mu }\mathit{.}$(22)Here, r_λ is the ratio of the line opacity to the total opacity, which is in this case identical to one. Taking the derivative of Eq. (22)yields $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{J}}_{\mathit{l,ij}}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}& \mathrm{=}& \hspace{0.17em}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{\int }}_{-1}^{\mathrm{1}}(\frac{\mathit{\partial }{\mathit{I}}_{\mathit{l}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}{\mathit{\phi }}_{\mathit{l,ij}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\\ & & \end{array}$(23)that naturally splits in two distinct terms. Because the second term involves only derivatives of the line profile, which depends strictly on local quantities, it can be computed with relative ease, and will be referred to by ${\mathrm{\Phi }}_{\mathit{l,ij}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{\int }}_{-1}^{\mathrm{1}}{\mathit{I}}_{\mathit{l}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}\frac{\mathit{\partial }{\mathit{\phi }}_{\mathit{l,ij}}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\hspace{0.17em}\mathrm{d}\mathit{\lambda }\hspace{0.17em}\mathrm{d}\mathit{\mu }\mathit{.}$(24)To evaluate the first term, we expand the derivatives of the local intensity through the derivatives of opacity and emissivity, $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{l}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}& \mathrm{=}& \sum _{{\mathit{l}}^{\mathrm{\prime }}}\sum _{{\mathit{i}}^{\mathrm{\prime }}}(\frac{\mathit{\partial }{\mathit{I}}_{\mathit{l}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{\chi }}_{\mathit{l}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}\frac{\mathit{\partial }{\mathit{\chi }}_{\mathit{l}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{n}}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{,}{\mathit{i}}^{\mathrm{\prime }}}}\\ & & \mathrm{+}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{l}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{\eta }}_{\mathit{l}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}\frac{\mathit{\partial }{\mathit{\eta }}_{\mathit{l}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{n}}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{,}{\mathit{i}}^{\mathrm{\prime }}}}){\mathrm{ℛ}}_{{\mathit{l}}^{\mathrm{\prime }}{\mathit{i}}^{\mathrm{\prime }}\mathit{k}}\\ & & \mathrm{+}\hspace{0.17em}\sum _{{\mathit{l}}^{\mathrm{\prime }}}\sum _{{\mathit{i}}^{\mathrm{\prime }}}\sum _{{\mathit{i}}^{\mathrm{\prime \prime }}\mathit{<}{\mathit{i}}^{\mathrm{\prime }}}{\mathit{p}}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{,}{\mathit{i}}^{\mathrm{\prime }}{\mathit{i}}^{\mathrm{\prime \prime }}}\mathit{,}\end{array}$(25)where $\begin{array}{ccc}{\mathit{p}}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{,}{\mathit{i}}^{\mathrm{\prime }}{\mathit{i}}^{\mathrm{\prime \prime }}}& \mathrm{=}& \left(\frac{\mathit{\partial }{\mathit{I}}_{\mathit{l}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{\chi }}_{\mathit{l}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}\frac{\mathit{\partial }{\mathit{\chi }}_{\mathit{l}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{\phi }}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{,}{\mathit{i}}^{\mathrm{\prime }}{\mathit{i}}^{\mathrm{\prime \prime }}}}\mathrm{+}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{l}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{\eta }}_{\mathit{l}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}\frac{\mathit{\partial }{\mathit{\eta }}_{\mathit{l}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{\phi }}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{,}{\mathit{i}}^{\mathrm{\prime }}{\mathit{i}}^{\mathrm{\prime \prime }}}}\right)\\ & & \end{array}$(26)describes the influence of changes of the line profile in the transition i′i′′ at depth point l′ on the intensity at depth point l, so that the response of the specific intensity to a change in q_k can thus be written in the simple form $\frac{\mathit{\partial }{\mathit{I}}_{\mathit{l}}\mathrm{\left(}\mathit{\mu ,\lambda }\mathrm{\right)}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\mathrm{=}\sum _{{\mathit{l}}^{\mathrm{\prime }}}\sum _{{\mathit{i}}^{\mathrm{\prime }}}{\mathit{a}}_{\mathit{l}{\mathit{l}}^{\mathrm{\prime }}\mathit{i}{\mathit{i}}^{\mathrm{\prime }}}{\mathrm{ℛ}}_{{\mathit{l}}^{\mathrm{\prime }}{\mathit{i}}^{\mathrm{\prime }}\mathit{k}}\mathrm{+}\sum _{{\mathit{l}}^{\mathrm{\prime }}}\sum _{{\mathit{i}}^{\mathrm{\prime }}}\sum _{{\mathit{i}}^{\mathrm{\prime \prime }}}{\mathit{p}}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{,}{\mathit{i}}^{\mathrm{\prime }}{\mathit{i}}^{\mathrm{\prime \prime }}}\mathit{.}$(27)Substituting Eq. (27)back into Eq. (23), and integrating over angles and wavelengths yields $\frac{\mathit{\partial }{\mathit{J}}_{\mathit{l,ij}}}{\mathit{\partial }{\mathit{q}}_{\mathit{k}}}\mathrm{=}\sum _{{\mathit{l}}^{\mathrm{\prime }}}\mathrm{(}{\mathit{A}}_{\mathit{l}{\mathit{l}}^{\mathrm{\prime }}\mathit{ij}}{\mathrm{ℛ}}_{\mathit{ljk}}\mathrm{+}{\mathit{A}}_{\mathit{l}{\mathit{l}}^{\mathrm{\prime }}\mathit{ii}}{\mathrm{ℛ}}_{\mathit{lik}}\mathrm{)}\mathrm{+}\sum _{{\mathit{l}}^{\mathrm{\prime }}}{\mathit{P}}_{{\mathit{l}}^{\mathrm{\prime }}\mathit{,ij}}\mathrm{+}{\mathrm{\Phi }}_{\mathit{l,ij}}\mathit{,}$(28)where A and P correspond to angle and wavelength integrated quantities a and p, respectively. Substituting this back into Eq. (12)yields for each state i and depth level l the equation $\begin{array}{ccc}& & \sum _{\mathit{j}}\mathrm{[}{\mathit{T}}_{\mathit{l,ji}}{\mathrm{ℛ}}_{\mathit{ljk}}\mathrm{-}{\mathit{T}}_{\mathit{l,ij}}{\mathrm{ℛ}}_{\mathit{lik}}\mathrm{+}\mathrm{(}{\mathit{n}}_{\mathit{l,j}}{\mathit{B}}_{\mathit{ji}}\mathrm{-}{\mathit{n}}_{\mathit{l,i}}{\mathit{B}}_{\mathit{ij}}\mathrm{)}\\ & & \mathrm{\times }\sum _{{\mathit{l}}^{\mathrm{\prime }}}\mathrm{(}{\mathit{A}}_{\mathit{l}{\mathit{l}}^{\mathrm{\prime }}\mathit{ij}}{\mathrm{ℛ}}_{\mathit{ljk}}\mathrm{+}{\mathit{A}}_{\mathit{l}{\mathit{l}}^{\mathrm{\prime }}\mathit{ii}}{\mathrm{ℛ}}_{\mathit{lik}}\mathrm{)}\\ & & \end{array}$(29)On the left-hand side of this equation, we find the level responses ℛ_ljk, together with a number of coefficients that do not depend on the quantity q_k (i.e., temperature, velocity, magnetic field ...) with respect to which we have taken the derivative. The right-hand side contains all the other terms, including those that are explicitly dependent on q_k, but no terms that depend on ℛ_ljk.

We can thus express Eq. (29)as a system of coupled linear equations of the form $\mathit{â}{x}\mathrm{=}{b}\mathit{,}$(30)where x are the level population responses ℛ_ljk, and â describes level interdependencies. This is a linear system of equations with ND × NL unknowns, where ND is number of depth points and NL is number of levels in the considered atom. Crucially, matrix â does not depend on q_k in any way, so it can be computed once, decomposed or inverted, and then used to solve the response to many different atmospheric quantities q at different heights k efficiently and quickly.

The right-hand side, b, contains the (local) responses of the atmospheric properties (collisional rates, profile functions, etc.) to q_k. In the case of pure line transfer, there are three contributors:

• the derivative of the collisional rates, which is in firstapproximation strictly local (i.e., responds only to perturbationin point l). Also, it is reasonable to assume thatcollisional rates respond only on temperature and pressure;

• the derivative of the local line profile, which influences the integration over wavelengths and angles to obtain the scattering integral. This factor is local as well and depends practically on all the atmospheric parameters (temperature, pressure, velocity, magnetic field);

• the derivative of line profiles at other points l′ (also strictly local), which influences opacity and emissivity in points l′, which then in turn influence intensity I_l through the transfer of radiation. This factor also depends on all the atmospheric parameters.

In the general case, bound-free transitions must be considered in the statistical equilibrium, and in special situations even overlapping transitions, not to mention other contributors to opacity and emissivity (electron and Rayleigh scattering, bound-free, free-free, and H-opacity, etc.). A generalization to such conditions is straightforward, and in principle only results in a number of extra terms in Eq. (26)arising from additional terms in the expressions for the opacity and emissivity. These terms, however, inevitably introduce off-diagonal terms on the left-hand side that make it intrinsically more difficult to invert. When such terms can be treated in LTE, however, these terms appear on the right-hand side instead, which makes them trivial to include. Such complications are, however, an integral part of the NLTE problem itself and are not specific to the method presented here.

Electron density deserves special mention, as it is a quantity which couples the chemical equilibrium and statistical equilibrium equations. If the electron density is assumed to be in LTE, the response is easily computed explicitly and included in the computations. A fully consistent approach would require finding the analytical derivatives of both the equations of chemical and statistical equilibrium. We postpone this discussion for the moment. We note that using finite differences for computing response functions automatically takes into account all NLTE effects on the responses of electron density.

2.3. Advantages of analytical response functions

The brute-force alternative to computing the response functions analytically is to do it numerically, that is, to compute ${\mathit{R}}_{{\mathit{q}}_{\mathit{k}}}\mathrm{\approx }\frac{{\mathit{I}}_{\mathrm{0}}\mathrm{(}\mathit{\mu ,\lambda }\mathrm{;}{q}\mathrm{-}\mathrm{\Delta }{\mathit{q}}_{\mathit{k}}\mathit{/}\mathrm{2}\mathrm{)}\mathrm{-}{\mathit{I}}_{\mathrm{0}}\mathrm{(}\mathit{\mu ,\lambda }\mathrm{;}{q}\mathrm{+}\mathrm{\Delta }{\mathit{q}}_{\mathit{k}}\mathit{/}\mathrm{2}\mathrm{)}}{\mathrm{\Delta }{\mathit{q}}_{\mathit{k}}}$(31)for a very small Δq_k, where each computation of I₀ requires one full solution of the NLTE problem. There are two major problems related to this approach.

• 1.

  If we want to compute the response to perturbations at all thepoints, we need NLTE solutions, which takesmuch longer than solving the linear system fromEq. (29). If we are somehow parameterizing ouratmosphere using a reduced number of parameters (as, forexample, when using nodes in depth-dependent inversions), thisproblem is mitigated somewhat, but it still requiresapproximately ten NLTE solutions, compared tojust one + the solution of the linear system for theanalytical approach.

• 2.

  The NLTE problem is highly nonlinear, and thus to obtain precise response functions, Δq_k must be kept small. However, the change of the NLTE solution induced by Δq_k must exceed the residual convergence error in the solution for the derivative to be accurate. This means that the NLTE solutions must be converged to a much higher precision than justified based solely on the estimated accuracy of the numerical scheme or the convergence state of the inverted atmosphere, a problem that does not affect the analytical method.

Although at first sight the first point seems to computationally strongly favor the analytic method, the use of only a small number of nodes in most inversion codes seems to reduce this advantage considerably, as thus far the method calculates the response to q_k, that is, a single quantity at a single height only. However, an important consequence of the left-hand side of Eq. (29)not depending in any way on q_k is that it allows us to calculate the combined response of different quantities at different heights by adding up the right-hand side for each of them, and then solving the linear system using the sum. This actually makes the method particularly suitable for use in inversion codes, as the problem that a perturbation of a parameter in a single node perturbs the value of that parameter over an extended height range in the atmosphere can be easily incorporated at little extra cost.

As in the case of response functions in LTE, however, a clear drawback of this method is the need to explicitly include any and all sources of opacity and emissivity in the calculation of the responses, and to follow their dependencies to the appropriate level. Failure to do so can result in significant discrepancies between the true response and calculated ones, which will not be apparent unless an explicit comparison with a response calculated using a finite difference approach is made.

3. Results

For demonstration purposes, we implemented the above procedure in a computer code that solves the NLTE radiative transfer problem using the MALI formalism of Rybicki & Hummer (1991). We chose a short characteristics-based formal solver which uses second-order Bezier splines for the interpolation of the source function, as described by de la Cruz Rodríguez & Piskunov (2013). While this choice of formal solver brings much-needed accuracy and stability to the solution of the NLTE problem, it causes some complication in the explicit computation of ∂I_l/∂χ_l′ and ∂I_l/∂η_l′, which we discuss in detail in Appendix A.

Our main aim is to compare the response functions computed with the analytical approach against those computed numerically (i.e., using Eq. (31)). We start with a simple example and then work our way up to more complicated ones.

3.1. Pure line transfer

3.1.1. Two-level atom

In the first example we consider a very simple model for the formation of the Hα line in a semi-empirical FALC model of the solar atmosphere (Fontenla et al. 1993). The line is modeled as a two level atom, with levels corresponding to the second and third energy level of hydrogen. We assume that the population of the ionized state is fixed to the LTE value, which means that the NLTE effects are only redistributing electrons between the two energy levels of the line. Emissivity and opacity are due to the line absorption and emission only. Here, we focus on the responses of the level populations and the emergent intensity to the temperature only, as our main aim is not to study the line response in detail, but only to compare the responses computed using the analytical approach with those obtained using the FDC approximation, and those computed using finite differences. In the following, we show “relative” responses, that is, the response of a certain quantity divided by the quantity itself. To show the response of the emergent intensity, we plot $\mathit{R̅}\mathrm{\left(}{\mathit{h}}_{\mathit{k}}\mathit{,\lambda }\mathrm{\right)}\mathrm{=}\frac{\mathit{R}\mathrm{\left(}{\mathit{h}}_{\mathit{k}}\mathit{,\lambda }\mathrm{\right)}}{{\mathit{I}}_{\mathrm{0}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}}\mathrm{·}$(32)The response functions to temperature computed using the three different approaches are shown in Fig. 1. For the numerical computation we used a perturbation of 1 K and converged the NLTE problem down to the greatest relative change in populations of 10^-8.

Fig. 1 Intensity response functions to temperature for a two-level atom line in FALC model atmosphere normalized with respect to the emergent profile, given in units of 10-4. Top left: numerical (finite difference) computation of intensity responses; top right: responses computed analytically using the method described in the paper; bottom left: responses computed using FDC approach; bottom right: relative differences, normalized with respect to the maximum response, between analytical and numerical computations (in log scale).

The line intensity response functions resemble those computed by Socas-Navarro & Uitenbroek (2004) (who plot absolute, normalized response functions, so negative responses cannot be seen). We see that the FDC approximation (bottom left) severely overestimates the response to temperature at heights around 1500 km. The overestimation occurs because these approximations are computed using the same assumptions as for LTE response functions, that is, the populations respond strictly locally and the response is governed by Saha-Boltzmann statistics. In extreme examples of NLTE lines (such as Hα), this is a very bad approximation as the level populations are predominantly determined by the radiation field coming from lower layers. It can be seen that a response function computed analytically using the method described above (top right) is visually indistinguishable from the one computed using finite differences. The bottom right panel of Fig. 1 shows the difference between the analytical and numerically computed response functions, divided by maximum absolute value of the numerical response function: $\mathit{r}\mathrm{\left(}\mathit{h,\lambda }\mathrm{\right)}\mathrm{=}\frac{\mathit{R̅}\mathrm{num}\mathrm{\left(}{\mathit{h}}_{\mathit{k}}\mathit{,\lambda }\mathrm{\right)}\mathrm{-}\mathit{R̅}\mathrm{an}\mathrm{\left(}{\mathit{h}}_{\mathit{k}}\mathit{,\lambda }\mathrm{\right)}}{\max \mathrm{\left(}\mathrm{\right|}\overline{\mathit{R}}\mathrm{num}\mathrm{\left(}{\mathit{h}}_{\mathit{k}}\mathit{,\lambda }\mathrm{\right)}\mathrm{\left|}\mathrm{\right)}}\mathrm{·}$(33)Here r(h,λ) is on the order of 0.001 (0.1%) everywhere in the atmosphere, both for this example and the two considered below. This is a level of agreement which suggest that the analytical response functions describe the line formation process not only qualitatively, but also quantitatively with sufficient accuracy to make use of them in an inversion code.

The time needed to set up, compute, and de-compose the matrix on the left-hand side of Eq. (29)is, in our specific implementation and for this specific example, smaller then one full NLTE solution, whereas for a model described with N_P parameters, an order of N_P full NLTE solutions are needed to numerically compute the response to all the parameters. Clearly, in this case, the use of analytically computed response functions reduces the computing time significantly.

3.1.2. Multilevel atom

We, very briefly present some results obtained using a four-level hydrogen-like atom. The atom consists of a ground level, which is assumed to be in LTE, and three NLTE levels with atomic constants corresponding to levels 2–4 of hydrogen. As before, we only consider line processes, the difference with respect to the previous example is thus only that the effects of nonlinearity are now more pronounced as we are dealing with a multilevel atom. Figure 2 shows the response function of the Hβ line, computed as in the previous example, with a similar level of agreement everywhere in the (λ,h) plane.

Fig. 2 Comparison between numerically (left) and analytically (middle) computed response functions to temperature for emergent intensity in Hβ line, normalized with respect to emergent intensity, in units of 10-4. Right: absolute difference between the analytical and numerical approach, normalized with respect to the maximum response.

3.2. Calcium II 8542 line

The infrared triplet of singly ionized calcium around 850 nm are among the most promising candidates for obtaining atmospheric diagnostics (see Quintero Noda et al. 2016, for an in-depth discussion of response function and diagnostic capabilities) in the solar chromosphere. The most frequently used member of this triplet, the line at 8542 Å, is one of the lines that is routinely inverted using the NLTE inversion code NICOLE, but the numerical cost of calculating the response functions makes such inversions very costly, a situation that has resulted in some authors resorting to approximate approaches (Beck et al. 2015). We thus investigate the potential of the analytical response functions for use in inversion codes.

We consider a FALC atmospheric model, without magnetic fields or velocities, and consider a CaII atomic model with five levels, as in Quintero Noda et al. (2016). For simplicity, we do not consider line asymmetry due to the presence of isotopes, but this effect should be straightforward to include. Although we now also consider other sources of opacity, namely H- bound-free and free-free opacity, electron scattering, Rayleigh scattering on neutral hydrogen, and bound-free opacity of the CaII atom, we simplify the problem somewhat by considering the electron density to be in LTE.

As discussed in Sect. 2.2, we need to take into account all the processes that contribute to the radiative rates, but also to the opacity and emissivity. This means that Eq. (26)has to be modified to account for continuum processes. By accounting carefully for the sources of opacity and emissivity and their derivatives, both the left-hand side of Eq. (29), which describes the coupling between all levels at all points in the atmosphere, and the right-hand side can be appropriately modified and the level responses can be obtained. This process is laborious but straightforward, after which the corresponding opacity and emissivity responses are readily computed and propagated to yield the response of the emergent intensity.

Fig. 3 Intensity response functions to temperature for a 8542 Ca line in FALC model atmosphere, normalized with respect to the emergent profile, given in units of 10-4. Top left: emergent line profile; top right: numerical (finite difference) computation of intensity responses; bottom left: responses computed analytically using the method explained in the paper; bottom right: absolute differences between analytical and numerical computations normalized to the maximum response (log scale).

Figure 3 shows the emergent line profile and response function to the temperature, computed with finite differencing and analytically. The emergent profile is somewhat different from the one shown in Quintero Noda et al. (2016), which could be due to the absence of magnetic field in our atmospheric model (the mentioned authors consider a 500 Gauss magnetic field), or because of differences in the treatment of the collisional rates. However, this difference is not critical here, as here we are predominantly interested in the response and not in the accuracy of the physics. We see that the intensity response function computed with the analytical approach and finite-differences are practically indistinguishable, and that their relative differences are fairly small. The response functions are not identical to those in Quintero Noda et al. (2016); in particular, they lack the negative “dip” very high in the atmosphere (see their Fig. 6). However, this feature is not present in the response functions calculated using either method, and thus points to differences in the atmospheric model, rather than shortcomings in the calculation of the response function itself.

4. Discussion and conclusions

In this paper we outline an analytical approach for the computation of the intensity response functions for lines formed in nonlocal thermodynamic equilibrium (NLTE). We take the derivatives of the rate equations analytically and then follow them through to get a large linear system where the matrix on the left-hand side describes the coupling between levels and points in the atmosphere (and hence has the dimension of NL × ND) and the right-hand side describes “local” dependencies (derivatives of the line profile, collisional rates, etc.). The left-hand side matrix does not depend on the quantity with respect to which the response is to be calculated, so it can be decomposed once and then used to obtain responses to different model parameters at different depths. Once responses of level populations are known, we use them to compute appropriate responses of opacity and emissivity, which are in turn propagated to get the response of emergent line profile.

In this paper we have restricted ourselves to the scalar (i.e., polarization-free) case and to the response functions to temperature. We have considered three examples: a pure line transfer example for a two- and four-level atom, and a five-level Ca II atom where we focused on the response of Ca II 8542 line. For all the cases considered we obtained excellent agreement between the intensity response function computed using the finite-difference approximation, and the new analytical approach.

Although a downside of the new method is that all interdependencies of the statistical equilibrium must be explicitly taken into account when calculating the derivatives, something that must be carefully checked using finite differences, the advantage of the analytical approach is that it is much faster. For our specific implementation, it takes less time than one NLTE solution, while a finite difference calculation of the same requires NLTE solutions where ND is number of depth points in the atmosphere.

While this work has its own interesting applications from the standpoint of theoretical radiative transfer, its main intended area of application is that of NLTE inversions. Current state-of-the-art codes use numerical response functions, which are time consuming to calculate, since typically on the order of ten NLTE solutions are needed for one set of derivatives, for each iterative step in minimization procedure. The method presented here would require a computing time similar to that of one NLTE solution and would thus offer an acceleration of one order of magnitude, and a corresponding inversion code is under construction.

Acknowledgments

We thank Smitha Narayanamurthy and Rafael Manso Sainz for stimulating and useful discussions.

References

Appendix A: Computing intensity responses to opacity and emissivity perturbations

A necessary ingredient in the computation of the response function for the intensity is the computation of the derivative of the specific monochromatic intensity in a given point with respect to opacity and emissivity everywhere in the atmosphere, that is $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{l}}}{\mathit{\partial }{\mathit{\chi }}_{\mathit{l}}^{\mathrm{\prime }}}& & \end{array}$and $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{l}}}{\mathit{\partial }{\mathit{\eta }}_{\mathit{l}}^{\mathrm{\prime }}}\mathit{,}& & \end{array}$where, for reasons of clarity, we have omitted the dependence on direction and wavelength. The specific intensity depends on the emissivity and the opacity everywhere in the atmosphere through the radiative transfer equation, which in most cases needs to be solved numerically. Thus, the value of the intensity will depend not only on the opacity and emissivity, but also on the selected numerical scheme.

It is important to realize that because different radiative transfer solvers yield different results, the use of different solvers will result in different inversion results. We can only try to ensure that our solver is sufficiently accurate for the question at hand and in light of the observational uncertainties.

After adopting a specific formal solver, the task is to compute the response of the intensity in point l to an infinitesimal perturbation of opacity and emissivity in point l′. del Toro Iniesta (2003) proposed an analytic approach, which, in the scalar case, yields $\frac{\mathrm{d}\mathit{\delta I}}{\mathrm{d}\mathit{s}}\mathrm{=}\mathrm{-}\mathit{\chi \delta I}\mathrm{+}{\mathit{\eta }}_{\mathrm{eff}}\mathit{,}$(A.1)where ${\mathit{\eta }}_{\mathrm{eff}}\mathrm{=}\mathit{\delta \eta }\mathrm{-}\mathit{\delta \chi I}$(A.2)and quantities with the prefix δ refer to perturbed quantities. Equation (A.1)is then solved by a suitable numerical method, typically the same one that was used to solve the radiative transfer equation to obtain the spectrum itself. This approach would be exact if the numerical solution were infinitely precise. This is, however, not the case, as integrating different functions on a discrete grid yields different inaccuracies.

An accurate result can be obtained by following the actual process used to obtain the numerical formal solution and systematically propagating the derivative of each expression through to the end. For this work we use a second-order Bezier solver, as in de la Cruz Rodríguez & Piskunov (2013), where the intensity in point l is computed using ${\mathit{I}}_{\mathit{l}}\mathrm{=}{\mathit{I}}_{\mathit{l}\mathrm{-}\mathrm{1}}{\mathrm{e}}^{\mathrm{-}\mathrm{\Delta }}\mathrm{+}{\mathit{w}}_{\mathit{l}\mathrm{-}\mathrm{1}}\mathrm{\left(}\mathrm{\Delta }\mathrm{\right)}{\mathit{S}}_{\mathit{l}\mathrm{-}\mathrm{1}}\mathrm{+}{\mathit{w}}_{\mathit{l}}\mathrm{\left(}\mathrm{\Delta }\mathrm{\right)}{\mathit{S}}_{\mathit{l}}\mathrm{+}{\mathit{w}}_{\mathit{C}}\mathrm{\left(}\mathrm{\Delta }\mathrm{\right)}\mathit{C,}$(A.3)where Δ = | τ_l−τ_l−1 | and C is the control point, which is computed from the derivative of the source function with respect to the optical depth.

For two given arrays containing the values of χ and η for a particular direction and wavelength, it is necessary to first compute the optical depth scale, source function, and source function derivative with respect to the optical depth before the Bezier solver can be employed. The computation of the responses of the intensity with respect to the opacity and the emissivity has a few additional steps and requires the derivatives of all quantities in Eq. (A.3)with respect to the opacity and the emissivity. The process consists of the following steps:

• 1.

  given the opacity, numerically compute its spatial derivativeneeded to perform the spatial integration of the opacity, i.e., to getthe optical depth;

• 2.

  compute the response of the opacity derivative to unit perturbations in the opacity at each point;

• 3.

  from the values of the opacity and the spatial opacity derivative, compute the optical depth scale;

• 4.

  from the perturbations of the opacity and the perturbations of the opacity derivative, compute the perturbations of the optical depth scale;

• 5.

  from the emissivity and the opacity, compute the source function (S) and the perturbations to the source function;

• 6.

  from the values of the source function, compute the derivative of the source function with respect to the optical depth scale (dS/ dτ);

• 7.

  from the perturbations of the source function and the perturbations of the optical depth, compute the perturbations of the derivative of the source function;

• 8.

  using S and dS/ dτ, following the Bezier scheme, formally solve the radiative transfer equation to obtain the specific monochromatic intensity at each point in the atmosphere;

• 9.

  finally, using the responses of S and dS/ dτ, compute the response of the intensity at each point in the atmosphere to perturbations in the opacity and the emissivity in each point of the atmosphere.

The above procedure is cumbersome, especially once translated into actual equations, but straightforward, and it yields the responses (∂I_l/∂χ_l′ and ∂I_l/∂η_l′), which agree with those computed using finite differences down to the level of the nonlinearity of the solution as exposed by the finite step size of the finite difference approximation. After these values are computed, the response of the intensity at any point to any combination of perturbations in the opacity and the emissivity can be found easily through $\mathit{\delta }{\mathit{I}}_{\mathit{l}}\mathrm{=}\sum _{{\mathit{l}}^{\mathrm{\prime }}}\left[\frac{\mathit{\partial }{\mathit{I}}_{\mathit{l}}}{\mathit{\partial }{\mathit{\chi }}_{{\mathit{l}}^{\mathrm{\prime }}}}\mathit{\delta }{\mathit{\chi }}_{{\mathit{l}}^{\mathrm{\prime }}}\mathrm{+}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{l}}}{\mathit{\partial }{\mathit{\eta }}_{{\mathit{l}}^{\mathrm{\prime }}}}\mathit{\delta }{\mathit{\eta }}_{{\mathit{l}}^{\mathrm{\prime }}}\right]\mathit{.}$(A.4)We use this method for the construction of the left-hand side matrix (i.e., when computing level responses) and for propagating opacity and emissivity perturbations once the level perturbations are known. This method provides very precise response functions in the LTE case as well (in the case of LTE, level responses can be found immediately), and it is quite possible that it will enable better converging inversion schemes for LTE as well.

All Figures

Fig. 1 Intensity response functions to temperature for a two-level atom line in FALC model atmosphere normalized with respect to the emergent profile, given in units of 10-4. Top left: numerical (finite difference) computation of intensity responses; top right: responses computed analytically using the method described in the paper; bottom left: responses computed using FDC approach; bottom right: relative differences, normalized with respect to the maximum response, between analytical and numerical computations (in log scale).

In the text

	Fig. 2 Comparison between numerically (left) and analytically (middle) computed response functions to temperature for emergent intensity in Hβ line, normalized with respect to emergent intensity, in units of 10-4. Right: absolute difference between the analytical and numerical approach, normalized with respect to the maximum response.
In the text

Fig. 3 Intensity response functions to temperature for a 8542 Ca line in FALC model atmosphere, normalized with respect to the emergent profile, given in units of 10-4. Top left: emergent line profile; top right: numerical (finite difference) computation of intensity responses; bottom left: responses computed analytically using the method explained in the paper; bottom right: absolute differences between analytical and numerical computations normalized to the maximum response (log scale).

In the text