Recovering Thermodynamics from Spectral Profiles observed by IRIS (II): improved calculation of the uncertainties based on Monte Carlo experiments

The study of the chromosphere is critical to understand the solar atmosphere (Carlsson et al., 2019). Although this region of the solar atmosphere has been observed for decades, understanding it is still a challenge. This is due to several major issues: 1. the complex coupling between the radiation field and the local magnetic field and thermodynamic conditions, which means interpretation of radiation must consider non-local thermodynamic equilibrium; 2. the transitions from fully ionized plasma to partially ionized and back to fully ionized, and from dominated by the plasma to domination by the magnetic field; 3. the highly dynamic and highly structured nature of the chromosphere on small spatio-temporal scales, necessitating sub-arcsecond high-quality observations on timescales of seconds. During the last few decades, both theoretical and observational improvements have allowed us to gain a better knowledge of the chromosphere and the events that occur in this region (e.g. Scharmer et al., 2008; Vissers et al., 2015; Leenaarts et al., 2011, 2013a, 2013b; De Pontieu et al., 2014; Quintero Noda et al., 2016; de la Cruz Rodríguez et al., 2019; Carlsson et al., 2019; Centeno et al., 2021; De Pontieu et al., 2021; Ishikawa et al., 2021; Vissers et al., 2022; Trujillo Bueno & del Pino Alemán, 2022).

The Interface Region Imaging Spectrograph (IRIS, De Pontieu et al. 2014) has been providing high-resolution observations (free from seeing effects introduced by the Earth’s atmosphere) of the chromosphere through the near-ultra violet spectral range around the Mg II h & k lines since 2013. The IRIS wavelength range also contains the Mg II UV triplet lines (hereafter denoted as Mg II UV2&3). The Mg II h & k lines are optically thick lines, being sensitive to the conditions at the high- and mid- chromosphere (Leenaarts et al., 2013a, b; Pereira et al., 2013), while the Mg II UV2&3 lines typically form lower in the chromosphere (although under flaring conditions the line formation may be different Kerr et al. 2016; Rubio da Costa et al. 2016). The most reliable method to derive physical information along the optical depth from these lines is by the "inversion" of these lines. This involves an iterative process in which, at first, an initial atmosphere is assumed and the radiative transfer equations are solved considering the non-local thermodynamic equilibrium and the partial frequency redistribution of the radiation from scattered photons, leading to a refinement in the underlying atmosphere, followed by further iterations.

The state-of-the-art Stockholm Inversion Code (STiC, de la Cruz Rodríguez et al. 2016, 2019) is the only available code capable of inverting these lines under these conditions. However, the inversion of a single Mg II h & k profile is computationally expensive (2.5 CPU-hour per observed profile). To minimize this burden, we created the IRIS Inversion based on Representative profiles Inverted by STiC ($\mathbf{IRIS^{2}}\$, Sainz Dalda et al. 2019). This technique is based on the inversion of the representative profiles (RP) of a broad selection of observations taken by IRIS in the Mg II h & k lines. A RP is the averaged profile of those profiles belonging to a data set that share the same shape, i.e. a similar distribution of the intensity over a given spectral range. This shape - or profile - is the signature of the conditions in the solar atmosphere where the radiation comes from. Therefore, the RP is the average of those profiles sharing similar conditions in the Sun. It is natural to define a Representative Model Atmosphere for the atmospheric conditions associated with the RP, which is obtained from the inversion of the RP. The core of $\mathbf{IRIS^{2}}\$is the $\mathbf{IRIS^{2}}\$database, which has 3 components: i) the synthetic RPs ($RP^{syn}$); ii) their corresponding RMAs obtained from the inversion of RPs; and iii) the uncertainty of the thermodynamics variable $p$, $\sigma_{p}$, of the RMA associated with the observed $RP^{syn}$. The $\mathbf{IRIS^{2}}\$database consists of $\approx$ 50,000 items, obtained from: 1. clustering 312 data sets on different targets (observed by IRIS) by using the $k-means$ technique (Steinhaus, 1957; MacQueen, 1967); 2. inverting each RP with STiC; 3. obtaining the $RP^{syn}$, RMAs and $\sigma_{p}$. The physical information relies on the relationship between the $RP^{syn}$, the RMA, and the uncertainties of the latter ($\sigma_{p}$), while the statistical significance of $\mathbf{IRIS^{2}}\$is given by the selection of the datasets considered in the database, which takes into account different solar features, exposure times, and locations on the solar disk.

In the first publicly released version of $\mathbf{IRIS^{2}}\$, the uncertainty of a physical variable $\sigma_{p}$ was obtained using the expression¹¹1A formal derivation of this expression can be found by using the equations of section 2.3 in Sánchez Almeida 1997 and of sections 6.2 and 6.3 in Bellot Rubio 1998. (see and del Toro Iniesta & Ruiz Cobo 2016):

with $i=0,...,q$ the sampled positions in the wavelength $\lambda_{i}$, $w_{i}$ their weights, $\sigma_{i}$ the uncertainties of the observation (e.g. photon noise), $m$ the number of physical quantities in the model $\mathbf{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. In this study, we only consider the Intensity Stokes parameter, $I$.

The expression above is valid to calculate $\sigma_{p}$, however, practical cases using $\mathbf{IRIS^{2}}\$ show an underestimation in $\sigma_{p}$ in those regions where the line is sensitive to the changes in physical variable $p$, and an overestimation of $\sigma_{p}$ where the line is not so sensitive to those changes. This is due to the fact that the $R_{p}$ is calculated considering all the optical depths (or nodes) of the model $\mathbf{M}$, while for $I(\lambda_{i};\mathbf{M})^{syn\ RP@STiC}$ only variations in a selected number of nodes in the model $\mathbf{M}$ are considered. That means, $R_{p}$ (RF) encodes the information in all the optical depths, while the $RP^{syn}$comes from a model evaluated in selected optical depths. Thus, in the particular nodes where the line is more sensitive to changes in $p$, the $R_{p}$ will be larger than in those nodes where it is less sensitive, making the $\sigma_{p}$ ($\sim R_{p}^{-1}$) smaller in regions where the line is more sensitive and larger in the regions where the line is less sensitive. This is the expected behavior, but in practice, in many cases, the obtained $\sigma_{p}$ is very low (high) for the optical depths where the line is (not) sensitive to changes in the physical variable $p$.

In this paper, we present a new approach to calculate the uncertainties of the RMAs in the $\mathbf{IRIS^{2}}\$database initially presented by Sainz Dalda et al. (2019). In Section 2, we explain how these new uncertainties have been calculated using a Monte Carlo simulation approach. The criteria used to determine the uncertainties are presented in Section 2.2. In Section 3, we evaluate the results obtained with the new version of $\mathbf{IRIS^{2}}\$with those obtained from inversion using STiC. Finally, in Section 4 we present the main conclusions and limitations of the new $\mathbf{IRIS^{2}}\$database.

When we invert an observed profile there are several factors that introduce a randomness to some key elements in the inversion. First, the noise inherent to an observation, both the one associated to the distribution of photons detected by the instrument (i.e., Poisson noise for our NUV photons), and the one associated to the readout or other electronic variations in our detector.

In addition, the initialization of the iterative inversion process is usually randomized. Thus, the path started and followed during an inversion of an observed profile may be different from another independent inversion for the same profile, which may yield different results. To better understand the impact of this randomness in the initialization of the inversion, we can invert the same profile several times with different initializations. This Monte Carlo inversion approach to quantify the uncertainty was used for the first time, to the best of our knowledge, by Westendorp Plaza 1999. Another source of possible variability in the inversion results comes from the initial atmosphere model. To minimize this, the inversion code DeSIRe (Ruiz Cobo et al., 2022) uses several initial guess models to independently invert the same profile, selecting the best fit of all the fits prociced by each inversion. We have not considered this case in our study since the $\mathbf{IRIS^{2}}\$database was built with the results from the inversion of the RPs using an unique initial guess model (FALC, Fontenla et al. 1993), and that is the one we only consider in our Monte Carlo approach.

One other aspect to consider when estimating uncertainties is that the inversion technique is based on the minimization of

with $i=0,...,q$ the sampled wavelengths,$w_{i}$ their weights, $\sigma_{i}$ the uncertainties of the observation (e.g. photon noise)²²2Formally, the Equation 2 consider a weight and a noise per spectral position per profile. However, for computational reasons only one weight and noise level per spectral profile is provided for all the profiles inverted in a batch. In this study, a batch is all the $\widetilde{RP}^{syn}_{t_{exp},noi}$ at a given $t_{exp}$ at a given $\mu$., and $\nu$ the number of observables, i.e., the spectral samples. This is the weighted Euclidean distance between the observed (input) and the synthetic (output) profile, with the weight higher for those wavelengths that we are more interested in. However, as we will see below, this metric is not optimal for those cases where the dimension of the observation, i.e. the number of observed wavelengths, is high.

The method that we have used to estimate the uncertainties associated with the $RP^{syn}$-RMA takes into account all these issues.

Let us now discuss the impact of these new calculations on the uncertainties on the thermodynamic parameters from $\mathbf{IRIS^{2}}\$, and in particular some cases that highlight the difference between the previous and new approach, and the limitations of any uncertainty calculation.

The first row of Figures 4 and 5 shows the $RP^{syn}$ (first column) and the associated RMAs for $T,v_{los},v_{turb},$ and $n_{e}$ (in the second, thirth, fourth, and fifth column respectively) with the uncertainties calculated using Eq. 3 in Sainz Dalda et al. (2019), i.e., using the response functions. In the following rows, the first column shows $\widetilde{RP}^{syn}_{t_{exp},noi}$ for the $RP^{syn}$#1619 with $t_{exp}=1,4,8,$ and $30s$, with the 5 noise randomizations over-plotted; from the second to the fifth columns the same RMA thermodynamic variables as in the first row, but now showing two types of uncertainties. In blue, we show the uncertainties calculated using the standard deviation of those ${}^{MC}\widetilde{RMA}^{syn}_{t_{exp},noi}$ associated with the inverted profiles for $\widetilde{RP}^{syn}_{t_{exp},noi}$ that satisfy the condition (5). In grey we show the uncertainties derived from all 25 Monte Carlo experiments, i.e., ${}^{[25]}\widetilde{RMA}^{syn}_{t_{exp},noi}$. In each panel of $\widetilde{RP}^{syn}_{t_{exp},noi}$ the number of profiles used to calculate the uncertainty is indicated in black, and, as a reference, the number of profiles that satisfies $\chi^{2}\leq 3$ when that number is less than 10 is indicated in green.

The uncertainties in $T$ are relatively small between $-6\leq log(\tau)\leq-3$ for all the $t_{exp}$ in both examples (#1619 and #6395). When all the 25 Monte Carlo experiments are considered (in grey), we see some differences, with the largest difference at for $-7\leq log(\tau)$ and to a lesser extent around $log(\tau)=-5$. The former location is the region in the optical depth where neither the Mg II h & k nor the Mg II UV2&3 are sensitive to the variations in the thermodynamic variables. The latter is where the Mg II h & k lines are more sensitive to changes in the atmosphere. Therefore, we should expect some uncertainty in the atmosphere for inversion cases in which the $\widetilde{RP}^{syn}_{t_{exp},noi}$ are not well fitted, and also where the Mg II h & k lines are actually sensitive to variations in the thermodynamic parameters. For $-3\leq log(\tau)$, the uncertainties are usually larger, which makes sense since the IRIS Mg II h & k profiles barely encode photospheric information, i.e., these lines are not sensitive to variations in the thermodynamics at this optical depth range.

It is important to distinguish how the uncertainties are calculated in the considered methods. In the method using the RFs, a small variation in the atmospheric parameter is introduced at given optical depth, then the response function is obtained as the difference between the synthetic profile from the atmosphere with the slightly modified parameter with respect to the profile corresponding to the unperturbed model atmosphere (i.e., without variation of any physical parameter). This process is repeated for all the optical depths considered in the model atmosphere. Let us now consider how uncertainties are determined in our new method. First, we note that during the inversion of the profiles only some optical depths (nodes³³3The cycles and nodes used in this study are the same as the ones used in Sainz Dalda et al. 2019: the first cycle considers four nodes in temperature, and three nodes both $v_{turb}$ and $v_{los}$. The second cycle uses seven nodes in temperature, and four nodes both in $v_{turb}$ and $v_{los}$.) are considered. In the Monte Carlo approach, five full inversions for the five $\widetilde{RP}^{syn}_{t_{exp},noi}$ are executed to evaluate the reproducibility of the results, using the standard deviation of the resulting models as the uncertainties for the original model. In the first case (using RFs), the synthetic profiles come from a model atmosphere evaluated in all the optical depths with a small variation, while this is not the case in our new calculations: the variation of the model atmosphere during the inversion only occurs in the nodes. In some cases, the inversion code may find a good fit generating some variations in some nodes, and none (or negligible) in other nodes because the code is able to fit the input profile without variation in these nodes. This is why in some cases the uncertainties at $-2\leq log(\tau)$ are very small. This effect can be seen for $v_{los}$ for ${}^{MC}\widetilde{RMA}^{syn}_{t_{exp},noi}$#6395 with $t_{exp}=1s$ in comparison with the other $t_{exp}$. For the latter, the Mg II UV2&3lines are more well defined (less noisy) and the inversion code may be trying to introduce a variation in the nodes at $-3\leq log(\tau)$. This effect is also noted in the $v_{turb}$, and in the $n_{e}$ at $-1\leq log(\tau)$. The conclusion then is that when assessing the uncertainties, we have to be aware of the optical depths where the observed lines are mostly sensitive to different parameters. These regions are slightly different for different solar features (e.g. umbra, penumbra or plage), and different between the physical parameters (e.g., see Fig. 2 in de la Cruz Rodríguez et al. 2016).

The first row of Figure 6 shows maps of the uncertainty calculated using the RFs ($\sigma^{RF}$) of $T,v_{los},v_{turb},$ and $n_{e}$ at $\log(\tau)=-4$. The second and the third rows show respectively the uncertainties calculated using the selective Monte Carlo experiments ($\sigma^{selMC}$) , i.e. ${}^{[N]}\widetilde{RMA}^{syn}_{t_{exp},noi}$, and all 25 Monte Carlo experiments ($\sigma^{all25MC}$), i.e. ${}^{[25]}\widetilde{RMA}^{syn}_{t_{exp},noi}$. At this optical depth, in the plage and the umbra and extended penumbra or canopy the $\sigma^{RF}_{T}<\sigma^{selMC}_{T}<<\sigma^{all25MC}_{T}$, while for the $v_{los}$, $v_{turb}$ , and $n_{e}$ the $\sigma^{RF}>\sigma^{all25MC}>>\sigma^{selMC}$. This situation is however different at $\log(\tau)=-2$ (see Figure 7), where $\sigma^{RF}_{T}>\sigma^{all25MC}_{T}>>\sigma^{selMC}_{T}$, and for the $v_{los}$, $v_{turb}$ , and $n_{e}$ the $\sigma^{RF}>>\sigma^{all25MC}>>\sigma^{selMC}$. These two figures illustrate what we mentioned above. When we calculate uncertainties from the response functions (as in Sainz Dalda et al. 2019), the uncertainties may be too low for those optical depths where the lines are sensitive to changes in the thermodynamics (large RFs), while they may be unrealistically high for those optical depths where the line is barely sensitive to changes in the thermodynamics (small RFs). We find that, for the Monte Carlo approach, the variation with optical depth of the uncertainties is more moderated, and typically smaller for the selective criterion than when considering all 25 Monte Carlo experiments.

In this paper, we present and discuss a novel methodology and the results of applying a selective Monte Carlo approach to determine the uncertainties associated with the Representative Model Atmosphere (RMA) in the $\mathbf{IRIS^{2}}\$ data base. These new uncertainties represent more realistic values than the previously publicly released uncertainties (Sainz Dalda et al., 2019) which were based on response functions. This is because the uncertainties in our new approach have been calculated from the synthetic representative profiles $RP^{syn}$ considering the different sources of uncertainty in the whole process, i.e.: different exposure times, different noise randomization, and different inversion initializations. We define the uncertainty of a physical parameter associated with the pair $RP^{syn}$-RMA as the standard deviation of this parameter in the set of depth-stratified output models (from the Monte Carlo experiments) that satisfy the ad hoc selection criterion shown in Eq. (5). The latter expression is used to minimize the impact of the output models based on inversions that produce a bad fit (with the noisy synthetic profile associated with $RP^{syn}$). In general, at the optical depths where the Mg II h & k and Mg II UV2&3 lines are sensitive to variations in a thermodynamic parameter, the difference between considering all 25 Monte Carlo experiments instead of the number that satisfies expression (5) is very small. The difference is larger for optical depths where the lines are not sensitive to thermodynamic changes.

The new uncertainties will be available to the public in the $\mathbf{IRIS^{2}}\$data base, both for IDL and Python. The uncertainty calculated from the 25 Monte Carlo experiments will also be provided as an extra field in the new version of the $\mathbf{IRIS^{2}}\$data base. Therefore, the new data base will have the following elements:

• •

  $RP^{syn}$: 472 wavelength positions, from 2794 to 2806Å, with a spectral sampling of $\approx 0.025m$Å

• •

  RMA: depth-stratified $T,v_{los},v_{turb},$ and $n_{e}$, sampled at 39 optical depths (i.e., "heights" in the atmosphere) with $\delta log(\tau)=0.2$

• •

  $\sigma_{sel}$: depth-stratified $\sigma_{T},\sigma_{v_{los}},\sigma_{v_{turb}},$ and $\sigma_{n_{e}}$, sampled at 39 optical depths with $\Delta(log(\tau))=0.2$, obtained from the selected Monte Carlo experiments (selective mode). These values are given for $t_{exp}=1,4,8,$ and $30s$. Therefore, $4\times\sigma_{sel}$ values are in the database.

• •

  $\sigma_{all}$: depth-stratified $\sigma_{T},\sigma_{v_{los}},\sigma_{v_{turb}},$ and $\sigma_{n_{e}}$, sampled at 39 optical depths with $\Delta(log(\tau))=0.2$, obtained from the 25 Monte Carlo experiments (all-in mode). These values are given for $t_{exp}=1,4,8,$ and $30s$. Therefore, $4\times\sigma_{all}$ values are in the database.

• •

  $\mu$: from 0 to 1, starting from $\mu=0.05$ at steps of 0.10, as indicated in Table 1.

The different $\mathbf{IRIS^{2}}\$inversion tools that allow users to interface with this database will use these database elements for internal calculations. The inversion of an IRIS Mg II h & k data set will only return the closest $RP^{syn}$ to the observed profiles, the corresponding RMAs, and the uncertainties taking into account the $\mu$ and the $t_{exp}$ of the observation and the uncertainty mode (selective or all-in) chosen by the user.

We believe that the empirical methodology we have developed for $\mathbf{IRIS^{2}}\$will be useful for understanding the uncertainties associated with other or similar inversion approaches.