Recovering Thermodynamics from Spectral Profiles observed by IRIS: A Machine and Deep Learning Approach

To answer some of the major open questions about the solar atmosphere, it is critical to understand the physical conditions in the chromosphere. The chromosphere has been observed for decades from ground- and space-based telescopes. Particularly the Interface Region Imaging Spectrograph explorer (IRIS, De Pontieu et al. 2014) has observed more than ≈19,000 data sets at subarcsecond resolution in the Mg ii h and k spectral range, in the near-ultraviolet, since it was launched in 2013.

The formation of the Mg ii h and k lines has been studied using numerical calculations that include the effect of partial redistribution of scattered photons and 3D radiative transfer effects (Leenaarts et al. 2013a, 2013b; del Pino Alemán et al. 2016; Sukhorukov & Leenaarts 2017). Some spectral features such as the intensity and wavelength of the emission peaks and central reversal of those lines can potentially serve as proxies of the temperature, line-of-sight velocity (v_los), and their gradients in various regions of the chromosphere (Leenaarts et al. 2013b; Pereira et al. 2015). However, so far these proxies have only been studied for quiet-Sun-like conditions, and do not provide detailed height-dependent diagnostics.

One of the most successful methods to recover physical information from spectropolarimetric observations is through nonlinear fitting techniques, where the parameters of a model atmosphere are iteratively adjusted in order to match the emerging model intensities with the observed spectra. This procedure is commonly called an "inversion" even though it is not based on a formal inversion to the radiative transfer equation.

de la Cruz Rodríguez et al. (2016, 2019) have developed the STockholm Inversion Code (STiC), which assumes nonlocal thermodynamical equilibrium, plane-parallel geometry and includes partial frequency redistribution. This inversion code (IC) recovers a depth-stratified model covering the photosphere, chromosphere, and transition region from the inversion of spectropolarimetric observations. We have used STiC to invert the Mg ii h and k intensity data observed with IRIS. However, on average, the time needed to recover such information is about 2 ${CPU}-{hour}/{profile}$. Thus, to invert an IRIS map such as the one shown in Figure 1—which contains ≈220,000 spectra—takes ≈440,000 CPU–hours.

Zoom In Zoom Out Reset image size

To reduce this computationally prohibitive task and allow for the inversion of large fields of view and time series of data, we have created a framework based on the inversion results of Mg ii h and k profiles and several machine and deep learning techniques. This new approach allows us to reconstruct models (with similar accuracy as STiC) from any IRIS data set in a few minutes using a desktop machine. Accompanying this Letter, this code is publicly available.

In Section 2 we describe the foundations of the new framework. In Sections 3 and 4 we present how the novel inversion methods work. The first results and their validation are shown in Section 5. Finally, we discuss the advantage and limitations of the framework in Section 6.

We have created a database of Mg ii h and k profiles observed with IRIS using the Representative Profiles (RP) of 250 data sets of different solar features, such as quiet Sun, plage, sunspots, emerging flux regions, active regions, flares, coronal holes, and filaments. The RPs are obtained after applying a clustering technique (k-mean analysis, Steinhaus 1957; MacQueen 1967) to the spectral profiles of Mg ii h and k from the selected data sets.

Each data set in the database is clustered in 60 RPs, which we have inverted with STiC using the same inversion scheme for all RPs. The number of RPs was determined by hardware constraints. The inversion setup consists of two cycles. The first cycle considers four nodes⁷ in temperature, and three nodes both in microturbulence (v_turb) and line-of-sight velocity (v_los). The second cycle takes as input the output model from the first cycle, now using seven nodes in temperature, and four nodes both in v_turb and v_los. Each RP is inverted three times from a different set of initial parameters (randomization) in each cycle.

For each (observed) RP we obtain a synthetic RP (RP@STIC), which is the best match found by the IC between the observed and synthetic profiles, and the corresponding Representative Model Atmosphere (RMA). An RMA consists of the depth-stratifications of temperature (T in K), v_los in cm s⁻¹, v_turb in cm s⁻¹, gas pressure (p_gas in dyn cm⁻²), mass–density (ρ in g cm⁻³), electron density (n_e in cm⁻³), column mass (c_mass in g cm⁻²), and height (z in cm).

2.1. Physical Meaning of the RPs and RMAs

2.2. Building the Database

For any pixel of a given observation, e.g., the one shown in Figure 1, we look for the closest synthetic profile obtained by the StiC inversion of the RP in the IRIS Mg ii h and k database (${I}_{i}^{{syn}\ {RP}@{STiC}}$). To determine the closest profile we use the same loss function as in STiC,

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \frac{1}{\nu }\displaystyle \sum _{i=0}^{q}{\left(I{\left({\lambda }_{i}\right)}^{{obs}}-I{\left({\lambda }_{i},{\boldsymbol{M}}\right)}^{{syn}{RP}@{STiC}}\right)}^{2}\displaystyle \frac{{w}_{i}^{2}}{{\sigma }_{i}^{2}},\end{eqnarray} \tag{ 1 }$$

with i = 0, ..., q the sampled wavelengths, w_i their weights, σ_i the uncertainties of the observation (e.g., photon noise), and ν the number of observables, i.e., the spectral samples. A low value of χ² tells us whether the fit between the observed (I_obs) and synthetic profiles ${I}_{i}^{{syn}\ {RP}@{STiC}}$ is good. We explicitly denote the dependency of the synthetic RP on the parameters of the model atmosphere (${\boldsymbol{M}}$). Once the code has found the best match between the observed profile and a synthetic RP in the IRIS database, it associates the corresponding RMA of that (closest) synthetic RP to the pixel in our observation. For large data sets, this look-up table process may take a few minutes on a desktop machine. Then, the code provides a χ²-map (to indicate the goodness of the match between the observed and synthetic profiles in the database), the output model atmosphere, and the associated uncertainty of each variable of the model.⁹

The uncertainty of a physical quantity p is calculated following Equation (42) in del Toro Iniesta & Ruiz Cobo (2016):

$$\begin{eqnarray}&&{\sigma }_{p}^{2}=\displaystyle \frac{2}{{nm}+r}\displaystyle \frac{{\displaystyle \sum }_{i=1}^{q}\left[I{\left({\lambda }_{i}\right)}^{{obs}}-I{\left.{\left({\lambda }_{i};{\boldsymbol{M}}\right)}^{{syn}{RP}@{STiC}}\right)}^{2}\right]\tfrac{{w}_{i}^{2}}{{\sigma }_{i}^{2}}}{{\displaystyle \sum }_{i=1}^{q}{R}_{p}^{2}({\lambda }_{i})\tfrac{{w}_{i}^{2}}{{\sigma }_{i}^{2}}},\end{eqnarray} \tag{ 2 }$$

with m the number of physical quantities in the model ${\boldsymbol{M}}$ evaluated in n grid points along the solar atmosphere, r the number of physical quantities considered constant along that atmosphere, and R_p the Response Function (RF) of a Stokes parameter to the physical quantity p (Mein 1971; Landi Degl'Innocenti 1979; Ruiz Cobo & del Toro Iniesta 1992). The RF provides the sensitivity of a wavelength sample in a Stokes profile to (changes of) a physical quantity. Thus, we use expressions like: "the core of the Mg ii h and k lines is sensitive to the T in optical depths around log(τ)¹⁰  = −5, while the wings are sensitive to T in −5 < log(τ) < −1."

We note that the inversion code does not operate over every grid point of the atmosphere, but over the nodes (each of them usually affects several grid points simultaneously). Therefore, in the latter case the RF is usually larger as a larger section of the atmosphere is perturbed per node and the uncertainty becomes lower than the estimates obtained by perturbing each grid point.

We have named this new tool the IRIS Inversion based on Representative Profiles Inverted by STiC or just IRIS squared (${{IRIS}}^{2}$).

IRIS² relies on two fundamental concepts: (i) the relationship between the synthetic RPs and the RMAs, given by the inversion of the observed RPs by STiC, and (ii) because the IRIS database covers a large variety of solar features, the RPs and corresponding RMAs are a meaningful representation of the variety found in the chromosphere and upper photosphere.

Since the IRIS Mg ii h and k database includes synthetic profiles of the RPs and the corresponding RMAs, we have trained several deep neural networks (DNN) to reproduce this relationship. In deep learning jargon, a synthetic RP is the input layer, i.e., the intensities of the synthetic RP at the sampled wavelengths are the input nodes (473). The corresponding RMA is the output layer, i.e., the values of physical quantities along the atmosphere (39) are the output nodes. Once the DNN is trained, we are, in principle, able to predict the physical quantity through the atmosphere for a given IRIS Mg ii h and k profile.

We have considered T, v_los, v_turb, and n_e as independent variables with respect to the corresponding synthetic RP. That means, we have trained a DNN for each of these physical quantities. The DNNs have different topologies (number of hidden layers and nodes), loss functions, and dropout parameters (to avoid overfitting). All the DNNs we have built use a rectified linear unit as activation functions, and we use 80% of the IRIS Mg ii h and k database as a training set, and the remaining as a test set. Similarly, we have trained the uncertainties (along the atmosphere) of these physical quantities.

We have named this method deep IRIS squared or ${{\bf{deepIRIS}}}^{2}$. More detailed information about the used DNNs will be given in a follow-up paper.

To validate RPs@STiC, IRIS², and deepIRIS², we have inverted, with STiC, every other pixel of the IRIS Mg ii h and k observation of NOAA AR 12480 on 2016 January 14. We consider the STiC results as the ground truth, but we should note that the STiC results also depend on initialization and are not guaranteed to provide a global minimum of the loss function.

Some results of using the RPs@STiC method are shown in the bottom of Figure 1. The first row of Figure 2 shows T at log(τ) = −4 as a result of the inversion using STiC (left), IRIS² (center), and deepIRIS² (right). Figure 3 shows v_los (top), v_turb (middle), and log(n_e) (bottom). Animations showing the variation of these parameters as a function of depth in the atmosphere are available in the electronic version.

5.1. Discussion: The Reliable Uncertainty Range

One question we should answer to validate our results is, for a physical quantity at a given optical depth p(τ), how large is the unsigned difference between the value recovered by our method and the one obtained by STiC compared to the uncertainty estimated using Equation (2)? We define the uncertainty multiplication factor (UMF) as:

$$\begin{eqnarray}&&{\mathrm{UMF}}_{p}\ (\tau )=\displaystyle \frac{| {p}^{\mathrm{method}}(\tau )-{p}^{\mathrm{STiC}}(\tau )| }{{\sigma }_{p}^{\mathrm{STiC}}(\tau )}.\end{eqnarray} \tag{ 3 }$$

We have selected five regions of interest (RoIs) of 21'' × 21'' in the FoV: plage (PL), quiet Sun (QS), umbra (UM), superpenumbra (SP), and a mix of regions (Mix, see Figure 2) to help us interpret the results. The full FoV is also evaluated.

The uncertainty maps (σ-maps) of T obtained by STiC, IRIS², and deepIRIS² are shown in the second row of Figure 2. The first panel of the third row shows the χ² map for STiC. The other panels in the third row show the UMF for T using the IRIS² (center) and deepIRIS² (right) methods. We have intentionally selected the optical depth log(τ) = −4 because it illustrates several key issues.

Some regions in the σ-maps are better (lower values) for IRIS² and deepIRIS² than for STiC.¹¹ This is a direct consequence of the better fit obtained by IRIS² compared to STiC, as mentioned above. The χ² map shown in the bottom left panel of Figure 2 is normalized in such a way that the fit is bad in regions where χ² ≫ 1, indicating that p(τ) is likely suffering from large uncertainties or may be wrong; the fits are better/good in regions with χ² of order 1 or less. In the UMF_p-maps (bottom center and right panels of Figure 2) regions with UMF_p ≈ 1 have values of p(τ) as "accurate" as STiC, or even better if UMF_p ≪ 1. Our example at log(τ) = −4 indicates that care should be taken when interpreting the results in plage and the SP.

Figure 4 shows the behavior of T^method and T^STiC (in thick and dashed line, respectively, in the first row), the ratio between the uncertainties (${\sigma }_{T}^{\mathrm{method}}/{\sigma }_{T}^{\mathrm{STiC}}$, second row), and the UMF_p (third row) for the proposed methods. The numbers next to the RoIs labels in the legend are the spatially averaged normalized ${\chi }_{\mathrm{STiC}}^{2}$ and the ratio ${\chi }_{\mathrm{method}}^{2}/{\chi }_{{STiC}}^{2}$.

Zoom In Zoom Out Reset image size

The uncertainties ratios show that mostly all the methods show the same uncertainty as STiC. In some regions the other methods give better results, in other cases, the opposite is true, but on average (blue line) the behavior is very similar, or clearly better.

The UMF_p for RPs@STiC and IRIS² is ≲1 in all RoIs for −3.5 < log(τ) < −1. For plage, the UMF for −4.5 < log(τ) < −3 is ≲1, and ${\chi }_{{method}}^{2}\leqslant {\chi }_{\mathrm{STiC}}^{2}$; therefore, the values provided by these methods are more accurate than the ones from STiC. The situation is the opposite for the SP. All in all, the values provided by RPs@STiC and IRIS² are mostly valid where Mg ii h and k lines are sensitive to T, i.e., −5 < log(τ) < −1. For deepIRIS², we have to be cautious, even when the uncertainties ratio is ≲1 in all the RoIs, the difference in T is noticeable in some of the ROIs.

Figure 5 shows the UMF_p for v_los, v_turb, and log(n_e). All the methods have a similar behavior in all the RoIs for v_los and v_turb, but there are differences for log(n_e) in plage regions for the IRIS² and deepIRIS² methods. We note that while substantial differences of physical parameters between different methods can occur, these differences may often be smaller than the intrinsic uncertainty.

Zoom In Zoom Out Reset image size

In summary, any of the proposed methods provides these values as well as STiC—or even better—within the intrinsic uncertainties. However, when we interpret a physical quantity and its uncertainty provided by our methods we should be specially cautious in (i) regions showing large values in the χ²-map, i.e., where the fit between the observed profiles and the synthetic profiles is not good, and (ii) those optical depths where Mg ii h and k lines are less sensitive to variations of a physical parameter. Under those conditions, the uncertainty will be larger, e.g., two or three times larger than the uncertainties provided in the database (as the UMF values in Figures 2 and 3 suggest).

We have created and evaluated three novel methods to rapidly obtain the atmospheric physical quantities in the chromosphere and upper photosphere from the profiles of the IRIS Mg ii h and k lines. The methods presented are valid for any spectro(polarimetric) data as far as they can be inverted by a traditional inversion code. We note that IRIS² can be used for any IRIS observation that includes Mg ii h or Mg ii k (or both) lines.

We summarize the main advantages and disadvantages for the three methods:

• 1.  
  RPs@STiC: on average, it is the closest to STiC. However, we lose spatial information. This can be minimized by using a much larger number of RPs for each data set. It stills requires a proper inversion, which takes hundreds of CPU hours (e.g., 320 CPU hours for 160 RPs).
• 2.  
  IRIS²: it offers results as good as STiC on average, being slightly better or worse than the latter in some solar features. The spatial information is almost as good as the original, although some regions may show little variation if the profiles are associated with the same RP in the database. That can be minimized by including a larger variety of profiles in the database. This method is 10⁵–10⁶ times faster than STiC.
• 3.  
  deepIRIS²: it predicts values of v_vlos, v_turb, and n_e, as good as the ones obtained by STiC. The predicted values for T are not as good but acceptable. A more complex DNN architecture and larger training and test data sets can improve this. The results do not lose spatial information and they look spatially smooth. It is ≈10⁶ times faster than STiC.

As a result of our investigation, we conclude that IRIS² is currently the fastest, easy-to-use method to recover reliable information from the chromosphere and photosphere from IRIS Mg ii h and k data. While we are improving these methods with a new database that includes 160 RPs per data set, as well as more observations, we note that the current versions of the IRIS Mg ii h and k lines database, IRIS² (both in IDL and Python)and deepIRIS² are available to the community at http://iris.lmsal.com/analysis.html. We expect that our database can be applied to a wide variety of investigations that use IRIS data.

This work is supported by NASA under contract NNG09FA40C (IRIS) and the Lockheed Martin Independent Research Program. JdlCR is supported by grants from the Swedish Research Council (2015-03994), the Swedish National Space Board (128/15) and the Swedish Civil Contingencies Agency (MSB). This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (SUNMAG, grant agreement 759548). This Letter has benefited from discussions at a meeting of team 399 studying magnetic-field-regulated heating in the solar chromosphere at the International Space Science Institute (ISSI) in Switzerland. IRIS is a NASA small explorer mission developed and operated by LMSAL with mission operations executed at NASA Ames Research center and major contributions to downlink communications funded by ESA and the Norwegian Space Centre.