text:
\section{Introduction}\label{sec:intro}

EUV lines emitted by atomic species with a variety of characteristic formation temperatures provide a powerful probe of evolving physical conditions in the active solar atmosphere. The observed count rates (s$^{-1}$) from a specified pixel in an image, measured in a set of of $m$ specified wavelength channels, are

\begin{equation}\label{eq:DEM-def}
    I_i = \int K_i (T) \, \xi(T) \, dT \, ; \qquad i = 1, \cdots, m \,\,\,
\end{equation}
where $K_i(T)$ (cm$^5$~pixel$^{-1}$~s$^{-1}$) is the emissivity function corresponding to the one or more spectral lines that fall within the bounds of the observed wavelength channel and $\xi(T) = n^2 \, ds/dT$ (cm$^{-5}$~K$^{-1}$) is the \textit{differential emission measure}, a quantity that provides information on the physical state of the emitting region. Here $n$ (cm$^{-3}$) is the ambient density, $T$ (K) the temperature, and $s$ (cm) the position along the line of sight.

Early observations made with the Skylab Apollo Telescope Mount (ATM) captured the emission in a set of spectral lines formed at temperatures ranging from a few $\times 10^4$~K to more than $10^6$~K, within a set of $5.5^{\prime \prime} \times 5.5^{\prime \prime}$ pixels, revealing \citep{1978SoPh...57..373E} significant temporal correlations between the emissions corresponding to different temperature ranges during several small flares. \cite{1978SoPh...57..373E} also computed the weighted differential-emission-measure-related quantity ($\mu(T)$, in their notation)

\begin{equation}\label{eq:mu-T}
\left \langle \xi(T) \right \rangle_i = \frac{\int K_i(T) \, \xi(T) \, dT}{\int K_i(T) \, dT} = \frac{I_i}{\int K_i(T) \, dT}
\end{equation}
for each observed spectral line $i$, thus allowing a relatively straightforward estimate of the evolution of the $\xi(T)$ profile throughout the events studied. Relatively modest increases in $\langle \xi(T) \rangle$ over the durations of the various flares were found, which were interpreted as being due to a combination of a significant increase in density in the emitting volume, coupled with a steepening of the temperature gradient $dT/ds$ (e.g., from an enhanced thermal conductive flux) and hence reduced thickness of the emitting volume.

Subsequent measurements by the Hinode X-Ray Telescope \citep[XRT;][]{2007SoPh..243...63G} and EUV Imaging Spectrometer \citep[EIS;][]{2007SoPh..243...19C}, and by the Solar Dynamics Observatory's Atmospheric Imaging Assembly \citep[AIA;][]{2012SoPh..275...17L} and EUV Variability Experiment \citep[EVE;][]{2012SoPh..275..115W} have provided a very large amount of flare-event data from atomic species formed over a broad range of temperatures. Given the superior nature of this data, it is manifestly appropriate to go beyond the construction of average $\langle \xi(T) \rangle$ values corresponding to each wavelength channel using the prescription~(\ref{eq:mu-T}), and instead treat Eq.~(\ref{eq:DEM-def}) as an integral equation to be solved for the desired \textit{source function} $\xi(T)$, given the set of observed \textit{data vector} values $I_i$ and the prescribed \textit{kernel functions} $K_i(T)$ (shown in Fig.~\ref{fig:aia_t_resp} for the six EUV channels in the SDO/AIA instrument).

\begin{figure}[t]
\centering
\includegraphics[height=.47\textwidth,
angle=90]{tresp.pdf}
\caption{Temperature response functions for the SDO/AIA EUV channels \citep{2012SoPh..275...41B}.}
\label{fig:aia_t_resp}
\end{figure}

In principle this appears straightforward: Eq.~(\ref{eq:DEM-def}) can be discretized as an $m \times n$ matrix equation

\begin{equation}\label{eq:DEM-matrix}
    I_i = \sum_{j=1}^n K_{ij} \, \xi_j \, ; \qquad i = 1, \cdots, m \,\,\, ,
\end{equation}
where the $\xi_j$ are components of a vector of $\xi(T)$ values at a set of chosen temperatures $T_j; j = 1, \cdots n$ and the $m \times n$ matrix $\mathbf{K}$ has coefficients $K_{ij}$ that represent the emissivity functions integrated over the various temperature bins $[T_j, T_j + \Delta T_j]$. Naively, this has a straightforward solution in terms of the generalized inverse of the matrix $\mathbf{K}$:

\begin{equation}\label{eq:DEM-inverse-matrix}
    \xi_j = \sum_{i=1}^m ([\mathbf{K}^{T} \mathbf{K}]^{-1} \mathbf{K}^T)_{ji} \, I_i \, ; \qquad j = 1, \cdots, n \,\,\, .
\end{equation}
However, it is well understood \citep[e.g.,][]{1986ipag.book.....C} that, due to a significant inter-dependence between the rows of the matrix $\mathbf{K}$, this is an ill-posed inverse problem \citep{1963SvPhU...5..594T,1985InvPr...1..301B,1996ApJ...457..882S}. Because of the ill--conditioning of $\mathbf{K}$, the solution given by Eq.~(\ref{eq:DEM-inverse-matrix}) is usually corrupted by oscillations due to noise amplification effects, and it often yields (unphysical) negative values at one or more points. Hence straightforward application of Eq.~(\ref{eq:DEM-inverse-matrix}) typically does not yield physically useful solutions. In practice, viable solutions for $\xi(T)$ must instead be found by imposing additional ``reasonableness'' requirements (e.g., smoothness, boundedness, positivity) on the solution.

Given the considerable importance of the recovered $\xi(T)$ profile in probing the physics of the emitting plasma \citep[see, e.g.,][and references therein]{2022ApJ...939...19E}, numerous authors have proposed a variety of methods for generating physically acceptable $\xi(T)$ profiles. For a general overview of such methods we refer the reader to, for instance, Chapter~5 of \cite{2008uxss.book.....P}. Whatever the method chosen, quantitative uncertainties in the $\xi(T)$ profiles derived can straightforwardly be determined by constructing multiple realizations of the data, with each realization drawn from a range of values consistent with the level of noise in the data; synthesis of the $\xi(T)$ profiles obtained from multiple realizations of the data produces a ``confidence strip'' of $\xi(T)$ profiles \citep[cf.][]{2006ApJ...643..523B}.

Perhaps the most basic approach to solving Eq.~\eqref{eq:DEM-def} is the so-called ``EM Loci'' method, originally applied to observations of UV lines observed in stellar spectra by \cite{1987MNRAS.225..903J} and subsequently applied to solar observations from the SoHO Coronal Diagnostic Spectrometer \citep[CDS;][]{1995SoPh..162..233H} by \cite{2007ApJ...655..598C} and \cite{2007ApJ...658L.119S}, and to data from the Hinode EUV Imaging Spectrometer \citep[EIS;][]{2007SoPh..243...19C} by \cite{2009ApJ...694.1256T}. The basis of this method is to initially assume that the emitting plasma is isothermal with temperature $T_o$, so that the differential emission measure can be approximated by a Dirac delta-function $\xi(T) = EM \, \delta(T-T_o)$, where $EM$ is the total source emission measure (cm$^{-5}$). The problem then reduces to determining the location and strength of the Dirac delta-function, i.e., the correct temperature $T_o$ and corresponding emission measure $EM$. For such a delta-function form of $\xi(T)$, Eq.~(\ref{eq:DEM-def}) reduces to $I_i = EM \, K_i(T_o)$, so that dividing an observed intensity $I_i$ by the corresponding response function value $K_i(T_o)$ gives the emission measure $EM$. Using different assumed temperatures $T_o$ yields a (generally, concave upward) set of $EM(T_o)$ points corresponding to that single observation, and repeating this for different observing channels gives a set of such ``EM Loci.'' If the source were truly isothermal, then all the EM Loci would intersect at a single point, which would then inform both the source temperature $T_o$ and the emission measure $EM$. However, for a non-isothermal emitting region the various EM Loci will intersect at different $(T_o, EM)$ points. Further, since the observed emission in a given observing channel will have contributions from material at temperatures other than the assumed $T_o$, the actual differential emission measure $\xi(T)$ at any chosen temperature is necessarily less than that deduced for an isothermal plasma at that temperature. It follows that the set of local minima of the overlapping EM Loci curves forms an upper bound (``curtain'') to the $\xi(T)$ profile \citep[see, e.g., Fig.~17 of][]{2012A&A...539A.146H}. This can be useful in determining a starting point for iterative methods but, as we shall show below, is not necessary for our proposed Regularized Maximum Likelihood (RML) method.

A more physical approach, that allows \textit{a priori} for the non--isothermal nature of the region under observation, is to model $\xi(T)$ as a locally smooth function, using, for example, a discretized cubic spline fitting procedure. Such an approach was originally applied to data from the NRL High Resolution Telescope and Spectrograph \citep[HRTS;][]{1977ApOpt..16..870T} by \cite{1992MmSAI..63..767M}, from the Solar EUV Rocket Telescope and Spectrograph \citep[SERTS; see][]{1992SoPh..137...87N} by \cite{1996ApJS..106..143B}, and from the SoHO CDS by \cite{2000A&A...363..800P}. \cite{1998ApJ...503..450K} have presented a Bayesian approach to derive $\xi(T)$ using a Metropolis Markov Chain Monte Carlo (MCMC) method, which, similar to the cubic spline fitting of the above references, applies a local smoothing, but over a scale determined by correlating the line emissivity function $K_i(T)$ \citep[$G_{u\ell} (T)$ in the notation of][]{1998ApJ...503..450K} with a parametric Gaussian-based ``Mexican Hat'' wavelet function (their Eq.~(8)). Even for simulated data generated from rather complicated forms of assumed $\xi(T)$~functions (their Figs.~7 and~8), the recovered $\xi(T)$ have a generally good overall agreement with the initially assumed forms, but unfortunately with rather large uncertainties. \cite{1998ApJ...503..450K} argue that ``a careful selection of the spectral lines used to infer the DEM [i.e., $\xi(T)$] is needed in order to avoid `artificially' generating structure in the DEM.''

A fundamentally different approach to recovering the differential emission measure $\xi(T)$ involves assuming a parametric form and determining the best values of the associated parameters by forward-fitting to the observations. Such an approach, using a parametric form constructed from multiple Gaussian functions of $\log T$ (see their Eq.~(7)), was employed in the analysis of SDO/AIA data by \cite{2011ApJ...732...81A}. A similar parametric ``basis function'' approach was employed by \cite{2015ApJ...807..143C}, using a combination of finite-width Gaussians of $\log T$ (see their Eq.~(13)) and (isothermal) delta-functions. However, comparison of the $\xi(T)$ profiles derived from simulated data generated even from relatively simple functional forms of $\xi(T)$ showed on occasion spurious features at high temperatures \citep[e.g., Fig.~4 in][]{2015ApJ...807..143C}. This ``Basis Pursuit'' technique seeks the minimum number of non--zero basis coefficients needed to fit the data within plausible uncertainties and, to avoid (unphysical) negative solutions, it applies a positivity constraint to each of the non--zero coefficients returned. However, the underlying simplex optimization procedure sometimes fails (cf. remarks in Sect.~\ref{sec:double_gaussian}) to generate a solution that satisfies all the required constraints, leading to reconstructed differential emission measure maps (see Sect.~\ref{sec:dem-maps} below)
that contain a number of ``empty'' pixels where the method has been unable to find a solution.

The SITES (Solar Iterative Temperature Emission Solver) algorithm \citep{2019SoPh..294..135M} first normalizes the temperature response curves $K_i(T)$ for each AIA observing channel to produce a set of ``relative response curves'' that sum to unity in each of several temperature bins. Count rates obtained by forward processing of an initially assumed $\xi(T)$ profile with the normalized temperature response curves are then compared to the observed counts and the differences used to generate a correction term to $\xi(T)$, and this step incorporates a positivity constraint. This iterative process is then repeated until the $\xi(T)$ profile converges to acceptable limits. \cite{2019SoPh..294..136P} have expanded this concept into a ``Gridded SITES'' method in which pixels with similar intensities in all six EUV channels of the AIA are grouped together and a $\xi(T)$ for each group of pixels is found; this saves considerable computational time by avoiding the calculation of $\xi(T)$ profiles that are similar to those that have already been calculated for another pixel with similar data properties. Such a grouping method could obviously also be used to expedite other algorithms used to generate $\xi(T)$ profiles from datasets that comprise a set of spectral line intensities.

\cite{2010A&A...523A..44G} employed a Maximum Likelihood (ML; see  Sect.~\ref{sec:method-description}) approach, specifically termed a Bayesian Iterative Method \citep[BIM;][]{richardson1972bayesian}. They compared the results using this method with previous analyses of spectral line data from the SoHO SUMER \citep{1995SoPh..162..189W} instrument in order to validate the conclusions of \cite{2008ApJ...672..674L}, with the results of \cite{2010AstL...36...44S}, who had applied the BIM method (with fewer iterative steps) to spectral line data from the SPIRIT/CORONAS-F instrument, and lastly, to broadband soft X-ray data from the Hinode XRT to validate the conclusions of \cite{2009ApJ...698..756R} that were based on a parametric basis function technique similar to that of \cite{2015ApJ...807..143C} but using ``top-hat,'' rather than Gaussian, basis functions. Figs.~2 and~3 of \cite{2010A&A...523A..44G} show that the BIM method can recover rather complicated $\xi(T)$ functions (e.g., those with multiple peaks) with a reasonable degree of accuracy (note that, in their figures, a logarithmic scale is applied to the $y$-axes). However, we will show in Sect.~\ref{section:simulated_data} that ML approaches without additional regularization constraints, in general produce unreliable results when applied to data consisting of only a few spectral lines, as in the case of AIA data. Another ML method proposed by \cite{1975SoPh...45..301W} and \cite{1980SoPh...67..285S} was utilized in the interpretation of Hinode XRT data by \cite{2008AnGeo..26.2999S}, and of EIS and AIA data by \cite{2014A&A...564A.130M}. Certain similarities between the performances of this ML~method and those of a genetic algorithm have been noted \citep{2008AnGeo..26.2999S}, and the use of a genetic algorithm has also proven effective \citep{2000ApJ...529.1115M} as a preconditioning step to determine the subset of spectral lines for which the corresponding~$\mathbf{K}$~matrix has minimum condition number. After this preconditioning step, the profile of $\xi(T)$ is then determined by applying the \cite{1963SvPhU...5..594T} regularization method to this selected subset of data.

Several authors have utilized \textit{regularized} inversion methods; such methods impose \textit{a priori} information on the solution in order to suppress the amplification of data uncertainties and thus obtain a a relatively stable ``smooth'' recovery of the $\xi(T)$ profile \citep[see][]{1986ipag.book.....C}. Several forms of regularized inversion, such as zeroth-order regularization, second-order regularization and maximum entropy regularization have been employed and tested using simulated EUV spectral line emission data \citep[e.g.,][]{1997ApJ...475..275J}. For example, the second-order regularized inversion method of \cite{2012A&A...539A.146H} seeks to minimize the quantity

\begin{equation}\label{eq:regularized-equation}
\vert \vert \, \mathbf{I} - \mathbf{K} \, \mathbf{\xi} \, \vert \vert^2 + \lambda \, \vert \vert \, \xi \, \vert \vert^2 ~,
\end{equation}
where $\mathbf{I} = (I_1, \dots, I_m)$, $\mathbf{\xi} = (\xi(T_1), \dots, \xi(T_n))$, and $\lambda$ is the regularization parameter. High values of $\lambda$ smooth the solution by penalizing large deviations in $\xi(T)$, while lower values of $\lambda$ place greater emphasis on accurate fitting of the reconstructed $\xi(T)$ profiles to the observed data. Although these regularization approaches are generally superior to simple methods such as EM Loci, (unphysical) negative solutions can still result from over-fitting of (noisy) data; accordingly, \cite{2012A&A...539A.146H} chose the lowest value of $\lambda$ that is consistent with a non-negative $\xi(T)$ at all temperatures $T$. The addition of the regularization term resolves many of the problems with $\xi(T)$ profiles corrupted by unphysical artifacts and generally recovers the underlying $\xi(T)$ (with uncertainties) quite well. However, analysis of simulated SDO/AIA data constructed from an assumed Gaussian $\xi(T)$ profile showed that the method can underestimate the actual $\xi(T)$ at high temperatures and overestimate it at low temperatures \citep[e.g., Figs.~11 and~13 of][]{2012A&A...539A.146H}. \cite{2013A&A...553A..10H} applied this method to a number of ``interesting'' pixels within the SDO/AIA images for a solar eruptive event on 2010~November~3, and we shall return to these results in Sect.~\ref{sec:AIA-data-application} below.

\cite{2013ApJ...771....2P} have presented a ``Fast Iterative Regularized'' (FIR) method based on the $L^2$ norm regularization of Eq.~(\ref{eq:regularized-equation}), starting with construction of a $\xi(T)$ profile in terms of a number of prescribed basis functions $B_l(T)$, the coefficients $e_l$ of which are to be determined. \cite{2013ApJ...771....2P} also impose a positivity constraint on the recovered $\xi(T)$ values, although their six-step iterative procedure is rather involved. Their $\xi(T)$ reconstructions compare favorably with those of \cite{2012A&A...539A.146H} (their Fig.~13, noting that this comparison does not [as it possibly could] impose positivity on the reconstructions produced by the \cite{2012A&A...539A.146H} method). However, recovery of $\xi(T)$ profiles from simulated data constructed from simple Gaussian functions of $\log T$ shows (e.g., their Figs.~1 through~4) that the method tends to generate spurious features, especially at high temperatures.

In summary, it is apparent that both ML and regularized approaches each bring their own set of strengths (and also some weaknesses) to the overarching problem of constructing $\xi(T)$ profiles from a discrete set of optically-thin spectral line intensities. Strengths include a high degree of mathematical rigor, the absence of any need to assume an \textit{a priori} (parametric) mathematical form for $\xi(T)$, the ability to impose global (rather than local) smoothness constraints, the ability to authentically reconstruct model $\xi(T)$ profiles from simulated data and to determine uncertainties in the solution. However, these approaches, taken individually, also bring a number of weaknesses and/or concerns into play, including the possible generation of spurious features in the solution (particularly at high temperatures) and the inability to find a solution which satisfies the imposed constraints.

We therefore here present (Sect.~\ref{sec:method-description}) a new \textit{Regularized Maximum Likelihood} (RML) method for the generation of $\xi(T)$ profiles, which combines the favorable elements of the Maximum Likelihood approaches of \cite{1975SoPh...45..301W}, \cite{1980SoPh...67..285S}, \cite{2010A&A...523A..44G}, and the regularization approaches of \cite{2012A&A...539A.146H} and \cite{2013ApJ...771....2P}. We shall demonstrate, using simulated data generated from idealized single-Gaussian (Sect.~\ref{sec:single_gaussian}) and double-Gaussian (Sect.~\ref{sec:double_gaussian}) differential emission measure forms, that this mathematically rigorous method, which does not require an \textit{a priori} choice of a functional form for the solution, is nevertheless robust and \textit{always} generates a (physically required) non-negative solution. In Sect.~\ref{sec:AIA-data-application} we apply the method to the same set of representative pixels in a solar eruptive event previously studied by \cite{2013A&A...553A..10H}, and we compare the $\xi(T)$ reconstructions obtained for each of these representative pixels by several methods, noting the similarities and differences between the reconstructions obtained. Overall, we find that the RML method produces results that are generally consistent with those obtained using other methods and, unlike other methods, never generates ``outlier'' results. The method also regularly produces acceptably low reduced $\chi^2$ values, showing a fidelity to the data without generating unphysical features caused by unnecessary overfitting. In Sect.~\ref{sec:summary} we summarize the results obtained and offer prospects for the implementation of this powerful differential emission measure reconstruction method in the near-real-time prediction of solar activity.

\section{The Regularized Maximum Likelihood (RML) Method}\label{sec:method-description}

Let $g_i(x,y)$ denote the count rate per unit time detected in the $i$-th channel of the SDO/AIA telescope corresponding to emission that originates within an $0.6^{\prime \prime} \times 0.6^{\prime \prime}$ pixel centered at the location $(x,y)$ on the solar disk, $K_i(T)$ (cm$^5$~pixel$^{-1}$~s$^{-1}$) the temperature response function of the $i$-th AIA channel (see Fig.~\ref{fig:aia_t_resp}), and $\xi(T;x,y)$ (cm$^{-5}$~K$^{-1}$) the line-of-sight differential emission measure at location $(x,y)$. These quantities are related by the integral equation

\begin{equation}\label{eq:eq1}
g_i(x,y) \approx \int_{T} K_i(T) \, \xi(T;x,y)  \, dT ~,
\end{equation}
where the approximation is due to the fact that experimental data are inevitably corrupted by statistical (and possibly other sources of) noise. Once discretized, Eq.~(\ref{eq:eq1}) becomes

\begin{equation}\label{eq:discretized}
\mathbf{g} = \mathbf{K} \, \mathbf{\xi} ~,
\end{equation}
where we have denoted $\mathbf{g}_i \equiv g_i(x,y)$, $\mathbf{K}_{ij} \equiv K_i(T_j)$, and $\mathbf{\xi}_{j} \equiv \xi(T_j;x,y)$, $i=1, \dots, m$; $j=1, \dots, n$. The essence of solving the differential emission measure inverse problem involves finding a (physically acceptable) source vector $\mathbf{\xi}$ that, given an observed data vector $\mathbf{g}$, satisfies Eq.~\eqref{eq:discretized}) within the bounds of observational uncertainty.

We assume that data are affected by statistical noise only, i.e.,

\begin{equation}
\mathbf{g}_i \simeq \mathcal{P}((\mathbf{K}\mathbf{\xi})_i) ~,
\end{equation}
where $\mathcal{P}(\eta)$ denotes a Poisson random variable with a mean $\eta > 0$. The \textit{Maximum Likelihood problem} is to determine $\mathbf{\xi}^\ast$ such that

\begin{equation}\label{eq:ml-prob}
\mathbf{\xi}^\ast = \mathop{\mathrm{arg\:max}}_{\mathbf{\xi}\geq 0} P(\mathbf{g} \, | \, \mathbf{K} \mathbf{\xi}) ~,
\end{equation}
where $P$ is the likelihood function associated with a Poisson distribution, defined as

\begin{equation}\label{eq:poisson-distribution}
P(\mathbf{g} \, | \, \mathbf{K} \mathbf{\xi}) = \prod_i \frac{\exp(-(\mathbf{K} \mathbf{\xi})_i) \, (\mathbf{K} \mathbf{\xi})_i^{\mathbf{g}_i}}{\mathbf{g}_i !} ~.
\end{equation}

The condition~(\ref{eq:ml-prob}) is equivalent to minimizing the negative logarithm of the likelihood function, viz.
\citep[cf. the C-statistic;][]{1979ApJ...228..939C}

\begin{equation}\label{eq:ml_eq}
\mathbf{\xi}^\ast = \mathop{\mathrm{arg\:min}}_{\mathbf{\xi}\geq 0} \left [ -\ln(P(\mathbf{g} \, | \, \mathbf{K} \mathbf{\xi})) \right ] = \mathop{\mathrm{arg\:min}}_{\mathbf{\xi}\geq 0} D_{\mathrm{KL}}(\mathbf{g}, \mathbf{K} \mathbf{\xi}) ~,
\end{equation}
where

\begin{equation}\label{dkl}
D_{\mathrm{KL}}(\mathbf{g}, \mathbf{K} \mathbf{\xi}) = \sum_i \, (\mathbf{K} \, \mathbf{\xi})_i - g_i \log (\mathbf{K} \, \mathbf{\xi})_i
\end{equation}
is the Kullback--Leibler divergence \citep{bertero2008iterative} and we have ignored terms in Eq.~(\ref{eq:poisson-distribution}) that depend only on the (known) quantities $\mathbf{g}_i$. Now, the partial derivative of $D_{\mathrm{KL}}$ with respect to one of the unknowns $\xi_k$ is

\begin{equation}\label{dkxi-gk}
\begin{split}
\frac{\partial D_{\mathrm{KL}}}{\partial \xi_k} & = \sum_i \, \left [ \, \frac{\partial}{\partial \xi_k} \sum_j K_{ij} \xi_j - \frac{g_i}{(K \xi)_i} \frac{\partial}{\partial \xi_k} \sum_j K_{ij} \xi_j \, \right ]
= \\
&= \sum_i \left [ K_{ik} - \frac{g_i}{(K \xi)_i} K_{ik} \right ] = \sum_i K^{T}_{ki} \left [ {\bf{1}}_i -  \frac{g_i}{(K \mathbf{\xi})_i} \right ]~,
\end{split}
\end{equation}
where ${\bf{1}}$ is an $n$-vector with components all equal to unity, and hence the gradient

\begin{equation}\label{eq:dkxi-grad}
\nabla_\xi D_{\mathrm{KL}}(\xi) = \mathbf{K}^T \left( {\bf{1}} - \frac{\mathbf{g}}{\mathbf{K} \mathbf{\xi}} \right)~,
\end{equation}
where $\mathbf{g}/\mathbf{K} \xi$ is the vector with $i$-th component equal to the ratio of the respective components $\mathbf{g}_i$ and $(\mathbf{K} \mathbf{\xi})_i$. By using the Karush-Kuhn-Tucker (KKT) conditions \citep{kuhn2014nonlinear1,kuhn2014nonlinear2}, it can be shown that the minimization problem~(\ref{eq:ml_eq}) is equivalent to finding a solution of

\begin{equation}\label{eq:kkt}
\mathbf{\xi} \cdot \nabla_\xi D_{\mathrm{KL}}(\xi) = 0 ~, \quad \xi \geq 0 ~,
\end{equation}
where the multiplication and the inequality between arrays are to be interpreted component--wise. This equation can be solved by means of the iterative scheme

\begin{equation}\label{eq:ml-solve}
\xi^{(l+1)} = \frac{\xi^{(l)}}{\mathbf{K}^T{\bf{1}}} \mathbf{K}^T \left(\frac{\mathbf{g}}{\mathbf{K} \xi^{(l)}}\right) \,; \qquad \xi^{(0)} = {\bf{1}}~,
\end{equation}
starting with a ``gray'' trial vector with all unit values.
This iterative algorithm, also known as Expectation Maximization \citep{dempster1977maximum} or, in the context of astronomical image deconvolution, Richardson--Lucy \citep{richardson1972bayesian,lucy1974iterative}, has been applied to the solution of inverse problems in several different fields, including medical imaging \citep{shepp1982maximum}, microscopy \citep{sarder2006deconvolution}, and most recently hard X-ray imaging of solar sources \citep{2013A&A...555A..61B,massa2019count,pianabook}. We refer the interested reader to \cite{bertero2008iterative}, \cite{10.1088/2053-2563/aae109} and \cite{massa2021predictive} for a comprehensive overview of the ML method and its applications. We note that the ML method defined in Eq.~(\ref{eq:ml-solve}) is the same as the BIM~method proposed by \cite{2010A&A...523A..44G} (cf. their Eq.~(22)). Also, the ML~strategy is similar to the method proposed by \cite{1975SoPh...45..301W} and \cite{1980SoPh...67..285S}, the main difference being the adoption by the latter method of empirically--defined weight functions in the iterative process. Contrary to many of the methods discussed in Sect.~\ref{sec:intro}, the ML method naturally includes a positivity constraint on the solution, which prevents the generation of unphysical negative $\xi(T)$ values. However, as will be evident from the simulation results below, the ML strategy is not well suited for addressing the $\xi(T)$~reconstruction problem from EUV AIA~data only, since the low number of available data channels does not permit the application of a sufficient number of constraints; this generally leads to the presence of spurious features in the reconstructions. This suggests that it would be highly advantageous to incorporate \citep{green1990use} a penalty term to allow the inclusion of \textit{a priori} information on the solution; specifically, we incorporate a term which penalizes solutions with high $\xi(T)$ values at high temperatures in order to avoid the appearance of artifacts at the upper end of the temperature range under consideration. With this term added, problem~(\ref{eq:ml_eq}) becomes finding the solution for

\begin{equation}\label{eq:rml}
\mathbf{\xi}^\ast = \mathop{\mathrm{arg\:min}}_{\mathbf{\xi}\geq 0} \left [ D_{\mathrm{KL}}(\mathbf{g}, \mathbf{K} \mathbf{\xi}) + \lambda \sum_j \, T_j \, \xi_j \right ] ~,
\end{equation}
where $\lambda > 0$ is a regularization parameter and the penalty term represents the differential-emission-measure-weighted mean temperature. The rationale for this choice of the penalty term is to select the solution, among all the ones with similar accuracy in fitting the data, that results in a $\xi(T)$ profile that is most concentrated at low temperatures. Indeed, the penalty term can be written as $\overline{n} \times \sum_j n_j T_j \, (ds/dT)$, which is the mean source density $\overline{n}$ multiplied by a term proportional to the total thermal energy per cm$^2$ along the line of sight $\int 3 \, n k_B T \, ds$ (where $k_B$ is the Boltzmann constant). Hence the regularization constraint effectively seeks, among solutions that adequately fit the data, the one with the lowest thermal energy content.

For solving problem (\ref{eq:rml}), we note that the derivative $\partial/\partial \xi_k$ of the objective function now includes an additional term $\lambda \, T_k$, which accordingly modifies the gradient expression~(\ref{eq:dkxi-grad}). Incorporating this additional term, the iteration formula~(\ref{eq:ml-solve}) becomes

\begin{equation}\label{eq:rml_solve}
\xi^{(l+1)} = \frac{\xi^{(l)}}{\mathbf{K}^T{\bf{1}} + \lambda \, \mathbf{T}} \, \mathbf{K}^T \left(\frac{\mathbf{g}}{\mathbf{K} \xi^{(l)}}\right) \,; \qquad \xi^{(0)} = {\bf{1}} ~,
\end{equation}
where $\mathbf{T} = (T_1, \dots, T_m)$ is the array of temperatures selected for the discretization of Eq.~(\ref{eq:DEM-def}). Since the array $\lambda \, \mathbf{T}$ has only positive entries, at each iteration the estimated $\xi^{(l)}$ profile is multiplied component-wise by a non--negative array to yield a revised estimate $\xi^{(l+1)}$. Therefore, since the iterative scheme is initialized with an array $\xi^{(0)}$ with all positive entries, each iteration, and hence the final solution, automatically satisfies the (physically required) non--negativity requirement. The value of the regularization parameter $\lambda$ is determined according to the \cite{morozov1966solution} discrepancy principle: starting from an intentionally high value of $\lambda$ (which leads to an over--regularized solution with a correspondingly high value of the reduced $\chi^2$), we then progressively decrease $\lambda$ by a constant multiplicative factor (typically $2/3$), until we reach a value $\tilde{\lambda}$ at which the solution $\tilde{\xi}$ has a reduced $\chi^2$ lower than 1. This criterion proves to be almost always effective and leads to solutions which do not overfit the data. We term the method expressed by Eq.~(\ref{eq:rml_solve}) the \textit{Regularized Maximum Likelihood}~(RML) method.

We conclude this section with a few remarks.

\begin{enumerate}

\item First, in rare instances, especially in weak pixels with relatively poor statistics, it is possible that the Morozov discrepancy principle cannot be satisfied for \textit{any} value of the regularization parameter $\lambda \,$; in such cases we simply adopt the solution corresponding to the initially chosen value of $\lambda$. Further, the main drawback of the progressively-decreasing-estimate approach for determining $\lambda$ is that it involves performing several reconstructions in sequence and hence substantially increases the computational cost. In a future work we will explore techniques for determining \textit{a priori} the optimum value of $\lambda$, possibly exploiting the use of neural networks \citep[see, e.g.,][]{alberti2021learning}.

\item Second, we can straightforwardly apply the iterative formula (\ref{eq:rml_solve}) for each of the $N$ pixels in an AIA image in parallel by casting $\mathbf{g}$ and $\xi$ as matrices of dimension $m \times N$ and $n \times N$, respectively. We start by applying Eq.~(\ref{eq:rml_solve}) to every pixel and, for a given pixel, we terminate the iterative process once that pixel satisfies the discrepancy principle constraint $\chi^2 < 1$. Pixels for which acceptable solutions have been thus obtained are then discarded when we recompute the next solution with a decreased value of $\lambda$. This approach significantly reduces the total amount of computational time required to produce acceptable $\xi(T)$ profiles for each pixel in the image.

\item Third, similar to the ML method of \cite{2010A&A...523A..44G} and the regularized approach of \cite{2012A&A...539A.146H}, we can provide a quantitative estimate of the uncertainty of the RML~solution by means of the confidence strip approach. Specifically, we repeat the optimization procedure defined in Eq.~(\ref{eq:rml_solve}) 25~times, each using a Poisson noise perturbation of the data. Finally, we estimate the uncertainty of the $\xi(T)$ profile provided by the RML method by using these 25~reconstructions to compute the standard deviations of each temperature bin. For each AIA pixel, the selection of the regularization parameter value $\lambda$ is performed just once (according to the Morozov discrepancy principle) during the reconstruction process; this value of $\lambda$ is then used for each of the 25~reconstructions, thus speeding up the uncertainty estimation process.

\item Finally, to the best of our knowledge, this is the first time that a maximum--likelihood--type algorithm has been applied to the reconstruction of $\xi(T)$ profiles from AIA data exclusively.

\end{enumerate}

\section{Results}\label{sec:results}

Throughout this section, the $\xi(T)$ profiles will be generated as a function of the temperature $T$~(K). However, in order to highlight features in different temperature ranges, they will be plotted as a function of the logarithm (to base $10$) of the temperature.

\subsection{Simulated data}\label{section:simulated_data}

As a test of the ability of the various inversion algorithms to accurately recover the source differential emission measure function $\xi(T)$, several ``ground truth'' analytic forms of $\xi(T)$ were used to generate simulated SDO/AIA data in the five SDO/AIA EUV channels 94~\AA, 131~\AA, 171~\AA, 193~\AA, and 211~\AA. The 335~\AA\ EUV channel was excluded because of its relatively weak temperature response function (Fig.~\ref{fig:aia_t_resp}); it is the only AIA channel that does not dominate\footnote{the 94 \AA\ channel is never strictly dominant, but does have a response roughly equal to that of the 131~\AA\ and 193~\AA\ channels at $\log_{10} T \simeq 6.9$
(Fig.~\ref{fig:aia_t_resp}) and so is included in the analysis.} the instrument response at \textit{any} temperature (Fig.~\ref{fig:aia_t_resp}), and we have found that attempting to fit the data in this channel often introduces spurious features in the recovered solutions. This simulated data in these five channels was then inverted, using five of the various algorithms described above, namely:

\begin{enumerate}
    \item the Basis Pursuit (BP) technique of \cite{2015ApJ...807..143C};
    \item the Fast Iterative Regularized methodology of \cite{2013ApJ...771....2P} ;
    \item the regularized (REG) approach of \cite{2012A&A...539A.146H};
    \item the (un-regularized) Maximum Likelihood (ML) method of Eq.~(\ref{eq:ml-solve}); and
    \item the Regularized Maximum Likelihood (RML) approach of Eq.~(\ref{eq:rml_solve}) ,
\end{enumerate}
to produce recovered $\xi(T)$ forms, to be compared with the original input form as a measure of the accuracy, fidelity, and robustness of each method. The reconstructed $\xi(T)$ profiles for all methods involved the discretization of temperature into 12 bins; this number was chosen as being more than the number of data channels (five), but not so large as to unnecessarily exacerbate the non-uniqueness aspect of the solution.
The chosen temperatures (at the [logarithmic] center of each bin) are given by $\log_{10} T$~(K) = $5.85 \, (0.15 ) \, 7.50$.  For each method we calculated $\xi(T)$ profiles using 25 different realizations of the data, established by perturbing the data with Poisson noise; the similarities and/or differences between the $\xi(T)$ profiles reconstructed from each data realization allows a characterization of the robustness of the method.

We will find in Sect.~\ref{sec:AIA-data-application} that many AIA pixels during flares and/or solar eruptive events contain data that are
consistent with EUV emission generated by a single source containing material at a relatively small range of temperatures corresponding to quiet-Sun coronal temperatures ($\simeq 1.5 \times 10^6$~K), or by a two-temperature source with both quiet Sun and enhanced ($\simeq 10^7$~K) components. We will not dwell extensively on the various physical reasons for this (although we will offer plausible explanations for the appearance of both components in due course), but we will lean heavily on this observed two-component structure to inform the types of simulated data to be used as a test of the various differential emission measure reconstruction methods. Specifically, the methods should be able to clearly identify and adequately characterize both low-temperature and (if they exist) high-temperature components in the differential emission measure profile.  We therefore test the methods using simulated data generated from two main types of ``ground truth'' differential emission measure profiles: a single Gaussian function of $\log T$ with a centroid temperature around $10^{6.1}$~K -- $10^{6.4}$~K (Sect.~\ref{sec:single_gaussian}), and a double Gaussian (Sect.~\ref{sec:double_gaussian}) with both a $10^{6.1}$~K component and a higher temperature component with a centroid temperature ranging from $10^{6.6}$~K to $10^{7.2}$~K.

\begin{figure*}[ht]
\centering
\includegraphics[height=\textwidth, angle=90]{single_gauss_linear.pdf}
\caption{\textit{From top to bottom:} reconstructions of simulated Gaussian $\xi(T)$~profiles by means of the Basis Pursuit technique by \cite{2015ApJ...807..143C} (BP), the inversion method by \cite{2013ApJ...771....2P} (FIR), the regularization technique by \cite{2012A&A...539A.146H} (REG), the ML method defined in Eq.~\eqref{eq:ml-solve}, and finally our proposed RML technique. The black solid lines represent the ``ground truth'' configurations, which are Gaussian functions of $\log T$, all with total Emission Measure (EM) equal to $4\times 10^{21}$ cm$^{-5}$, standard deviation equal to 0.15 dex, and centroid temperatures given by $\log_{10}(T_c) = 6.1, 6.2, 6.3, 6.4$, from left to right, respectively. The gray lines represent the reconstructions from 25~different realizations of Poisson-noise-perturbed data; the red lines represent the mean values of these 25~reconstructions. The $\xi(T)$ profiles are plotted as a function of the (base $10$) logarithm of the temperature~(K).}

\label{fig:method-comparison-single-gaussian-linear}
\end{figure*}

\begin{figure*}[ht]
\centering
\includegraphics[height=\textwidth, angle=90]{double_gauss_linear.pdf}
\caption{\textit{From top to bottom:} reconstructions of a simulated double Gaussian $\xi(T)$ profile by means of the Basis Pursuit technique by \cite{2015ApJ...807..143C} (BP), the inversion method by \cite{2013ApJ...771....2P} (FIR), the regularization technique by \cite{2012A&A...539A.146H} (REG), the ML method defined in Eq.~(\ref{eq:ml-solve}), and finally our proposed RML technique. The black solid lines represent the ``ground truth'' configurations, which consist of double Gaussian profiles. The left peak is centered in $\log_{10}(T_c) = 6.1$, while the right peak is centered in $\log_{10}(T_c) = 6.6, 6.8, 7.0$ and $7.2$, from left to right, respectively. The left Gaussian has total EM equal to $4\times 10^{21}$ cm$^{-5}$ and standard deviation 0.1 dex, while the right one has total EM equal to $2\times 10^{21}$~cm$^{-5}$ and standard deviation 0.15 dex. The gray lines represent 25 different reconstructions from as many realizations of Poisson noise affecting the data; the red lines represent the mean values of the ten reconstructions. The $\xi(T)$ profiles are plotted as a function of the (base~$10$) logarithm of the temperature~(K).}

\label{fig:method-comparison-double-gaussian-linear}
\end{figure*}

\subsubsection{Single Gaussian forms}\label{sec:single_gaussian}

Fig.~\ref{fig:method-comparison-single-gaussian-linear} shows the results using simulated ``ground truth'' $\xi(T)$ that take the form of a Gaussian function of $\log T$. For each method we show the results for 25 different realizations of the data, established by perturbing the data with Poisson noise; the red line is the mean of these 25 reconstructions, while the gray lines show the individual reconstructions for each realization. From left to right, the various columns represent ``ground truth'' $\xi(T)$ forms with centroid temperatures given by $\log_{10} T_c = 6.1, 6.2, 6.3,$ and~$6.4$; each has a standard deviation equal to 0.15~dex.

\begin{enumerate}

\item The Basis Pursuit (BP) approach of \cite{2015ApJ...807..143C} (first row of Fig.~\ref{fig:method-comparison-single-gaussian-linear}) is very robust, as evident from the near-congruence of the gray curves for the 25 different data realizations. The method also reproduces the input $\xi(T)$ functions rather well for the low centroid temperature cases (first and second columns of Fig.~\ref{fig:method-comparison-single-gaussian-linear}); however, the fidelity of the reconstructed $\xi(T)$ forms becomes poorer for higher centroid temperatures (third and [especially] fourth columns of Fig.~\ref{fig:method-comparison-single-gaussian-linear}). Compared to the other inversion methods, BP appears to be the most accurate in reconstructing the input $\xi(T)$ with $\log_{10}T_c = 6.1, 6.2$ and $6.3$ (first three columns of Fig.~\ref{fig:method-comparison-single-gaussian-linear}). This is likely, in large part, due to the strong similarity between the shapes of the simulated ``ground truth'' configurations and the (Gaussian) basis functions used in this method.

\item The Fast Iterative Regularized (FIR) methodology of \cite{2013ApJ...771....2P} (second row of Fig.~\ref{fig:method-comparison-single-gaussian-linear}) is also very robust, but significantly overestimates the peak in $\xi(T)$ for the low centroid temperature cases, while underestimating it (and shifting the peak to lower temperatures) for cases with higher centroid temperatures.

\item For the regularized approach of \cite{2012A&A...539A.146H} (REG; third row of Fig.~\ref{fig:method-comparison-single-gaussian-linear}), there are no clear trends; however the method typically underestimates the value of the peak and/or shifts it to lower temperatures. This method is also very robust, i.e., insensitive to different realizations of the data, perhaps demonstrating an over-reliance on the regularization term in the minimization of equation~(\ref{eq:regularized-equation}).

\item The (unregularized) Maximum Likelihood (ML) method based on Eq.~\ref{eq:ml-solve}) (fourth row of Fig.~\ref{fig:method-comparison-single-gaussian-linear}) underestimates the peak $\xi(T)$ for the two lowest centroid temperature cases. It does much better for cases with a higher centroid temperature, and in all cases it well reproduces the value of that centroid temperature. However, it systematically creates a spurious artifact at temperatures $\log T \gtrsim 7.1$. The method is less robust than the others, with significant differences between the reconstructed $\xi(T)$ profiles for the 25 different data realizations, particularly on the high-temperature wing of the profile.

\item The Regularized Maximum Likelihood (RML) method (final row of Fig.~\ref{fig:method-comparison-single-gaussian-linear}) reproduces the peak intensities and centroid temperatures of the input $\xi(T)$ with $\log_{10} T_c = 6.1$~and~$6.2$ (first two columns of Fig.~\ref{fig:method-comparison-single-gaussian-linear}) with greater accuracy compared to FIR, REG and ML. However, the RML method does progressively underestimate the high-temperature wings (at $\log_{10} T \simeq 6.5$) in the third column, and significantly underestimates the $\log_{10} T = 6.4$ centroid temperature in the last column. Indeed, underestimation of $\xi(T)$ at such temperatures is a common element of all of the reconstruction methods tested (except ML). Reference to Fig.~\ref{fig:aia_t_resp} reveals why this may be the case: temperatures $\log_{10} T \simeq 6.3 - 6.7$ correspond to rather low values of the temperature response curves $K(T)$ in all channels, so that fitting the data in any observed channel is more straightforwardly accomplished by adding a relatively small amount of emission measure at other temperatures. Finally, the RML method seems to be less robust (i.e., more sensitive to perturbations in the data) compared to the other inversion techniques; however, the variance shown by the reconstructions (gray lines) associated with different data realizations still appears to be well within reasonable tolerances, with the exception of a single realization of the data, the reconstruction of which shows a small artifact at $\log_{10} T \simeq 6.7$).

\end{enumerate}

\begin{figure*}[t]
\centering
\includegraphics[height=\textwidth, angle=90]{aia_maps.pdf}
\caption{Maps recorded by AIA on November 3 2010 around 12:15:02 in the 94 \AA, 131 \AA, 171 \AA, 193 \AA, and 211 \AA\ channels, from left to right, respectively.}
\label{fig:aia_maps}
\end{figure*}

\subsubsection{Double Gaussian forms}\label{sec:double_gaussian}

Fig.~\ref{fig:method-comparison-double-gaussian-linear} shows the results obtained with simulated ``ground truth'' $\xi(T)$ that take the form of \textit{two} Gaussian functions of $\log T$. In application of the methods to real data (see Sect.~\ref{sec:AIA-data-application} below), it will become apparent that ascertaining the veracity of recovered high-temperature ($\log_{10} T_c \gtrsim 7.0$) components in $\xi(T)$ is important; hence, in the simulation test shown in Fig.~\ref{fig:method-comparison-double-gaussian-linear}, the low-temperature component is centered at $\log_{10} T_c = 6.1$, while the high-temperature component is centered at four different values of $\log_{10} T_c = 6.6, 6.8, 7.0$ and $7.2$ (from the first to fourth column, respectively). The low-temperature Gaussian has a total EM equal to $4\times 10^{21}$~cm$^{-5}$ and a standard deviation of 0.1 dex, while each of the high-temperature components has EM$ = 2 \times 10^{21}$~cm$^{-5}$ and a standard deviation of 0.15 dex. Again, the gray lines represent the reconstructions from 25 different data realizations, and the red lines represent the mean values of these 25 reconstructions.

\begin{enumerate}

\item The Basis Pursuit (BP) approach of \cite{2015ApJ...807..143C} generally underestimates the intensity of both low-temperature and high-temperature components, the latter being especially obvious in the first column, corresponding to a high-temperature component that has the lowest centroid temperature of the four cases considered. We note that, in the test case of the first column, the BP method reconstructed a \underline{zero} $\xi(T)$ profile for one of the 25 realizations of the noise in the data, probably because (as noted in Sect.~\ref{sec:intro}) occasionally the method cannot find a solution which both fits the data within uncertainties and satisfies the imposed constraints; the value of the average reconstruction plotted in red reflects only the 24 non-zero reconstructions obtained. Nevertheless, with the exception of the test presented in the first column, the BP method does correctly identify the approximate centroid temperatures of both low- and high-temperature components, and it is the most accurate in retrieving the peak value, width and centroid of the low--temperature component. The method is not as robust as it is for the single Gaussian case of Fig.~\ref{fig:method-comparison-single-gaussian-linear}, as is expected, given the more complex configuration that has to be reconstructed in this case). Indeed, there is a considerable range in the $\xi(T)$ values for the high-temperature component for the different data realizations.

\item The FIR approach \citep{2013ApJ...771....2P} typically underestimates the high-temperature component, but apparently compensates for this by significantly \textit{over}estimating the intensity of the low-temperature component. With regard to robustness, there is considerable variation in the reconstructions for different data realizations, particularly with reference to the high-temperature component.

\item The regularized inversion technique REG of \cite{2012A&A...539A.146H} significantly underestimates both low- and high-temperature components, especially for lower values of the centroid temperature of the high-temperature component (e.g., first column of Fig.~\ref{fig:method-comparison-double-gaussian-linear}), although it is the most robust, recreating very similar $\xi(T)$ profiles for the 25 different data realizations. Further, the method consistently reconstructs spurious high-temperature material at $\log_{10} T \gtrsim 7.2$. Similar artifacts are also present in the reconstructions shown in the first two columns of Fig.~\ref{fig:method-comparison-single-gaussian-linear}), although, in this case, they are less prominent.

\item The (unregularized) Maximum Likelihood (ML) approach significantly underestimates the magnitude of the low-temperature component, and \textit{severely} overestimates both the magnitude and width of the high-temperature component, in all cases creating features at temperatures well in excess of any that are actually present in the ``ground truth'' $\xi(T)$ profile. Interestingly, it also deduces these (erroneous) $\xi(T)$ profiles rather robustly.

\item The Regularized Maximum Likelihood approach (RML) tends to slightly overestimate the centroid temperature of the low--temperature component. However, it is the most accurate in reconstructing both the centroid and peak value of the high temperature components for the tests presented in the second and third columns. Further, only RML and ML suggest the presence of the (actual) high-temperature component in the test case considered in the first column, although neither the centroid temperature nor the peak value of this component are accurately reproduced. A comparison with the results obtained by means of the (un--regularized) ML method shows the effectiveness of the adopted regularization, which suppresses spurious high-temperature features present in the ML reconstructions. The RML method also proves to be quite robust, showing deviations larger than those of REG or FIR only in the case of the high temperature component of the test presented in the fourth column. For this test case, most of the RML reconstructions associated with the different data realizations show a high--temperature component which is narrower than the ground truth and the peak of which is shifted towards lower temperatures. Most importantly, comparing the single-Gaussian results of Fig.~\ref{fig:method-comparison-single-gaussian-linear} and the double-Gaussian results of Fig.~\ref{fig:method-comparison-double-gaussian-linear}, it is apparent that the RML method recovers a (suitably scaled and positioned) high-temperature component if and only if one is actually present.

\end{enumerate}

\begin{figure*}[h!]
\centering
\includegraphics[width=0.69\textwidth]{Figure_aia.pdf}
\caption{\textit{Top row:} AIA images recorded on 2010 November~3 around 12:15:02~UT in the 131~\AA~and 211~\AA~channels (left and right panel, respectively). The numbered crosses denote the location of the pixels selected for a comparison of the $\xi(T)$~profiles reconstructed by the different methods. \textit{Second to fifth row:} $\xi(T)$~profiles associated with Pixels~1 through~8. In each panel, the profiles reconstructed by the BP, FIR, REG and~RML methods are plotted in magenta, blue, green and orange, respectively. The $\pm 1 \sigma$ uncertainties associated with the REG and~RML reconstructions have been added as error bars at each temperature point in the pertinent reconstruction. The intensities of the reconstructions of Pixels~1, 3 and~5 have been multiplied by a factor of~$2$. The intensity axes are linear, while the temperature axes are in terms of $\log_{10} T$~(K).}
\label{fig:aia_dems}
\end{figure*}

\begin{table*}[pht]
\centering
\begin{tabular}{rcccccccc}
\toprule
Method  &Pixel 1 &Pixel 2 &Pixel 3 &Pixel 4 &Pixel 5 &Pixel 6 &Pixel 7 &Pixel 8 \\
\midrule
BP &2.13 &1.96 &1.05 &1.62 &2.13 &2.45 &1.49 &1.44\\
FIR  &1.11 &2.82 &0.57 &0.94 &3.98 &0.67 &1.11 &0.83\\
REG &1.03 &1.27 &1.15 &1.18 &0.90 &2.35 &1.16 &11.36\\
RML &0.55 &0.65 &0.35 &0.72 &0.99 &0.70 &0.86 &0.94\\
\bottomrule
\end{tabular}
\caption{Reduced $\chi^2$ values for the various pixels in the 2010~November~3 solar eruptive event (Fig.~\ref{fig:aia_dems}), for four different differential emission measure reconstruction algorithms.}
\label{tab:chi2}
\end{table*}

\begin{figure*}[ht]
\centering
\includegraphics[width=\textwidth]{DEM_BP.pdf}
\caption{DEM maps reconstructed from the AIA images shown in Fig.~\ref{fig:aia_maps} by the BP method \citep{2015ApJ...807..143C} . The temperature value corresponding to each DEM map is reported in the top--left corner of each panel. Further, the same (logarithmic) color map is shared by each of the panels in the same row.}
\label{fig:dem_bp}
\end{figure*}

\begin{figure*}[ht]
\centering
\includegraphics[width=\textwidth]{DEM_FIR.pdf}
\caption{Same information as in Fig.~\ref{fig:dem_bp}, but using the FIR method \citep{2013ApJ...771....2P}}
\label{fig:dem_fir}
\end{figure*}

\begin{figure*}[ht]
\centering
\includegraphics[width=\textwidth]{DEM_REG.pdf}
\caption{Same information as in Fig.~\ref{fig:dem_bp}, but using the REG method of \cite{2012A&A...539A.146H}.}
\label{fig:dem_reg}
\end{figure*}

\begin{figure*}[ht]
\centering
\includegraphics[width=\textwidth]{DEM_RML.pdf}
\caption{Same information as in Fig.~\ref{fig:dem_bp}, but using our proposed RML method.}
\label{fig:dem_rml}
\end{figure*}

\subsection{Application to AIA data}\label{sec:AIA-data-application}

Here we compare\footnote{Small differences between the results presented in this paper and those in \cite{2013A&A...553A..10H} result from the fact that we do not consider the AIA 335~\AA\ channel in reconstructing the $\xi(T)$ profiles. Although, as discussed in Sect.~\ref{section:simulated_data} above, the 335~\AA\ channel has a small emissivity function over the temperature range of interest (Fig.~\ref{fig:aia_t_resp}) and hence can lead to the introduction of spurious features in the recovered $\xi(T)$ profiles, this very consideration means that its influence on the recovered $\xi(T)$ profiles is not entirely negligible.} the $\xi(T)$ reconstructions for a selected set of pixels in the 2010~November~3 12:15~UT event previously studied by \cite{2013A&A...553A..10H}. Fig.~\ref{fig:aia_maps} shows the general morphology of the region, which produced both a flare and an erupting flux rope, at 12:15:02~UT in several SDO/AIA channels. Eight pixels (see Fig.~\ref{fig:aia_dems}) were selected by \cite{2013A&A...553A..10H} as representative of different features in the field of view. In Fig.~\ref{fig:aia_dems} we show the $\xi(T)$ reconstructions for each of these pixels by the various reconstruction methods (excluding the problematic unregularized ML method) considered above. Below we discuss and compare the general features of the results obtained.

\begin{itemize}

    \item \textit{Pixel 1 -- Core of the erupting plasmoid:} The various reconstructions are quite similar, and clearly show two components to the emission, a relatively cool one at around $1.5 \times 10^6$~K and a hotter one at around $10^7$~K. We agree with the interpretation of \cite{2013A&A...553A..10H} that these components probably represent the background corona along the line of sight, and the plasmoid emission, respectively. The methods provide consistent reconstructions of the low--temperature component. For the high--temperature component, the centroid is consistently retrieved by the different techniques, while its reconstructed peak value, width and skewness vary somewhat from method to method.

    \item \textit{Pixel 2 -- Filament/stem behind the plasmoid:} Here the relatively cool background line-of-sight emission is about $3 \times $ more intense than for Pixel~1, while the hot component is about $3 \times$~less intense. Although all methods agree on the temperature of the two components, the FIR method appears to broaden the low-temperature component and shift its centroid temperature downward slightly.

    \item \textit{Pixel 3 -- High corona away from the event:} This pixel contains a relatively small amount of emission (peak value of $\xi(T) \simeq 2 \times 10^{20}$~cm$^{-5}$~K$^{-1}$) at a temperature around $\log_{10} T = 6.2$, similar to the temperature of the background corona components in Pixels~1 and~2. The BP, REG and RML methods agree on the centroid temperature of this component, although BP produces an enhanced high--temperature wing and RML returns an higher estimate of the peak value of $\xi(T)$. However, the reconstructions by BP and REG both fall within the $\pm 1 \sigma$ error bars associated with the RML reconstruction, suggesting that this is simply due to statistical uncertainty in the data. The FIR method again (cf. Pixels~1 and~2) produces a broader peak at a lower temperature than the other methods.

    \item \textit{Pixel 4 -- Corona away from the event:} All methods agree very well as to the presence and intensity of a relatively cool $\log_{10} T \simeq 6.2$ component; however, considerations similar to those for Pixel 3 apply to the enhanced high--temperature wing produced by BP and the higher estimate of the peak flux returned by RML.

    \item \textit{Pixel 5 -- Corona near the event:} Again, all methods agree very well as to the presence and intensity of a relatively cool $\log_{10} T \simeq 6.2$K component; however, the FIR method substantially overestimates the low-temperature ``wing'' of this component relative to the other three methods. Further, similar to Pixels~3 and~4, the BP method reconstructs a $\xi(T)$ profile that is more skewed towards higher temperature values.

    \item \textit{Pixel 6 -- Low corona flare emission:} All methods agree very well with regard to both the relatively cool ($\log_{10} T \simeq 6.2$) and hot ($\log_{10} T \simeq 7.1$) components present. Again, the FIR method creates a low-temperature $\xi(T)$ component that is broader (and skewed toward lower temperatures) than those of the other methods. Further, the BP reconstruction of the high--temperature component shows an enhancement in the high--temperature wing, while the other methods agree quite well in terms of the centroid value, peak intensity and width of the reconstructed high--temperature component.

    \item \textit{Pixel 7 -- Envelope just ahead of the plasmoid:} Here the FIR method once again produces a centroid of the low--temperature component that is broader and slightly shifted towards lower temperatures compared to the other reconstructions. Also, the BP method retrieves an enhanced high-temperature ``wing'' for this component. There is also evidence of a weak high-temperature ($\log_{10} T \simeq 7.0$) component. Given the location of this material, it is not unreasonable to expect additional heating there \citep[consistent with, e.g., the findings of][that ``the CME is in the heat-releasing state (i.e., entropy loss) throughout its journey from the Sun to Earth.'']{2020FrASS...7....1M}.

    \item \textit{Pixel 8 -- Further ahead of the plasmoid:} Here all methods agree that there is a low-temperature component with an intensity similar to that of Pixel~7. The BP method again shows an enhancement in the high--temperature wing of the $\xi(T)$ profile. The 131 \AA\ image shows that this pixel is at the leading edge of the erupting material, providing (similar to Pixel~7) a plausible explanation for such enhanced heating of the low-temperature component.

\end{itemize}

In any ill-posed inversion problem, there is a necessary tradeoff between accurately fitting the data and producing a plausibly smooth solution. The fit of the reconstructed $\xi(T)$ profiles to the AIA data is determined by the reduced $\chi^2$ values shown in Table~\ref{tab:chi2}. Most of the $\chi^2$ values for the eight pixels are of order unity, showing that the data is well-fit but not overfitted at the expense of the plausibility of the $\xi(T)$ profile. Notable exceptions occur in application of the FIR method to Pixel~5, where FIR creates a more intense peak compared to the reconstructions by the other methods, and in application of the REG method to Pixel~8, where REG provides a reconstruction with a significantly lower peak value compared to that of the other methods.

\subsection{Differential emission measure maps}\label{sec:dem-maps}

In Figs.~\ref{fig:dem_bp} through~\ref{fig:dem_rml} we show differential emission measure maps at different temperatures, constructed by assigning each pixel the value of $\xi(T)$ calculated for that pixel at the temperature in question. As pointed out by \cite{2013A&A...553A..10H}, for some temperatures the differential emission measure maps closely resemble the images from one of the AIA channels: for example, the $T = 1.0$~MK map closely resembles the 171 \AA\ image, and the $T = 11$~MK map closely resembles the 131 \AA\ image. These close matches reflect the well-defined peaks in the response curves for those channels at those respective temperatures (Fig.~\ref{fig:aia_t_resp}). However, the differential emission measure maps for other temperatures reveal features that are not as apparent in the individual SDO/AIA images, such as the large region of emission extending out to coordinates $(x \simeq - 1100, y \simeq [-450 : -350])$ above the Eastern limb in the $T = 4$~MK map.
This feature is reproduced consistently by all the reconstruction methods (see Figs.~\ref{fig:dem_bp} through~\ref{fig:dem_rml}), confirming its reality and highlighting the generally complicated relationship between temperature and emitted wavelength, particularly for SDO/AIA channels with a relatively broad temperature response. We also note (cf. remarks in Sect.~\ref{sec:intro}) that the DEM~maps produced by the BP method (see Fig.~\ref{fig:dem_bp}) contain $253$ ``null'' pixels, showing points where the simplex optimization method adopted by BP has not been able to find a solution that satisfies all the required constraints.

\section{Summary}\label{sec:summary}

The results of Sects.~\ref{section:simulated_data} and~\ref{sec:AIA-data-application} show that, both for simulated and actual AIA data, the differential emission measure $\xi(T)$ reconstructions obtained using the Regularized Maximum Likelihood (RML) method described in Sect.~\ref{sec:method-description} are broadly compatible with those of other methods; in no case does the RML~method create a $\xi(T)$ profile that is an ``outlier,'' and, with the possible exception of the very weak (and hence statistically uncertain) Pixel~3, in no case is its $\chi^2$ value unacceptably large (corresponding to over-smoothing) or unacceptably small (corresponding to over-fitting of data).

As evidenced by the reconstructed $\xi(T)$ profiles constructed from different (Poisson-noise) realizations of the data (Fig.~\ref{fig:method-comparison-single-gaussian-linear}), the RML~method\footnote{Both IDL and Python codes that implement the RML method are available at
\url{https://github.com/paolomassa/WAFFLE.git}.}
is robust. Further,  it is straightforward to implement and computationally efficient, taking\footnote{When a more efficient rule for the selection of the regularization parameter (which does not involve performing multiple reconstructions as is the case when the Morozov discrepancy principle is employed) is implemented, the computational time for reconstructing $500\times 500$ $\xi(T)$~profiles (without uncertainty estimation) should decrease to less than $10$~s.} about $35$~s to reconstruct $\xi(T)$ profiles for $500\times 500$ pixels ($\sim 7000$ $\xi(T)$~reconstructions per s) on an Apple MacBook Pro M1 (Chip Apple M1, CPU 8-core) processor, without considering need to reconstruct solutions for multiple data realizations in order to estimate the uncertainty on the solution. Further, it does not require an \textit{a priori} choice of a parametric functional form and, very importantly for this particular application, always generates a non-negative solution without the need to impose \textit{a posteriori} adjustments on the solutions obtained \citep[cf.][]{2013ApJ...771....2P}. We conclude that it is an appropriate  method to use in the construction of $\xi(T)$ profiles from pixel-by-pixel AIA data.

In future work we will develop a more general version of RML that is applicable to other data sets, such as those from EIS or XRT, in order to better constrain \citep[cf.][]{2013A&A...553A..10H} the overall solution of the differential emission measure reconstruction problem. Given the general features of the RML~reconstructions, in particular existence, fidelity, robustness, and non-negativity, this differential emission measure reconstruction method is particularly suitable for the application of machine learning tools to solar data.
Furthermore, the computational efficiency of the method will allow us to generate, from AIA data in near real-time, four-dimensional data hypercubes $\xi(x,y; T; t)$ \citep[see Sect.~5 of][]{10.3389/fspas.2022.1040099} that represent a time series of differential emission measure maps. Application of machine learning tools to such hypercubes, each labeled with the eventual level of activity (e.g., maximum GOES level, duration of high levels of emission, total energy released) that results in the event represented by the data hypercube in question, can be used to identify both morphological and thermodynamic precursors of solar activity \citep[see, e.g.,][]{2020arXiv201106433G}, thus providing a promising tool for predicting the timing, location, and characteristics of solar eruptive events.

\section*{Acknowledgements}
We thank Lars Hebenstiel for several useful discussions on the features, strengths and weaknesses of various differential emission measure reconstruction algorithms. We also thank the authors of the BP \citep{2015ApJ...807..143C} and FIR \citep{2013ApJ...771....2P} methods for making their codes freely available, thus allowing us to us them in the comparison of various methods. This work was supported by NASA Kentucky under NASA award number 80NSSC21M0362.

\bibliographystyle{aa}

redpajama_set_name: RedPajamaArXiv


text:
NRGI Impact: The Resource Governance Index Drives Reform
The Resource Governance Index Drives Reform (PDF 744.06 KB)
The Resource Governance Index Drives Reform (Card) (PDF 173.66 KB)
NRGI staff and partners have continued to use the Resource Governance Index (RGI) since its launch in 2017, with important steps taken to close the governance gaps identified in several NRGI focus countries.
The RGI has compelled governments to improve their performance. Guinea's Ministry of Mines and Geology engaged closely with NRGI to change the governance of the mining sector, leading to the publication of the first RGI interim evaluation, which highlights areas in which Guinea has improved since the 2017 edition was released. The Guinean government made many disclosures as a direct result of the RGI's findings.
In Tunisia, RGI findings led the Tunisian state-owned enterprise (SOE), la Compagnie des Phosphates de Gafsa (CPG), to engage with NRGI staff to understand its 2017 RGI score, which is notable given the traditional reluctance of Tunisian SOEs to engage with civil society actors. As a result of this engagement CPG has begun to disclose phosphate revenue information on its website.
In Azerbaijan, NRGI was asked to provide analysis to the ministry of finance, the state oil company and the state oil fund on how to improve the country's rating for the next iteration of the RGI. NRGI recommendations helped inform the basis of a transparency report prepared by auditing firm Baker Tilly for Azerbaijan's Commission on Transparency in the Extractive Industry.
Qatar Petroleum, the national oil company of Qatar, approached NRGI regarding its RGI score. The company's general counsel met with NRGI to discuss ways to improve performance and NRGI advised the SOE to publicly disclose annual and sustainability reports. Subsequent to the meeting, Qatar Petroleum published annual and sustainability plans for the past 14 and eight years respectively, for the first time, on its website.
Civil society organizations and journalists in dozens of resource-rich countries have used RGI data to diagnose governance successes and shortcomings. In Nigeria, BudgIt used the RGI to support advocacy around contract transparency and the governance of the SOE NNPC by creating infographics and sharing them through social media channels. In Ukraine, Dixi Group used the RGI to support its advocacy around the development of the "Law on Ensuring Transparency in the Extractive Industries," which was then passed. Newspapers from Mexico to Myanmar have referenced the RGI in articles about resource governance. In these instances, the most relevant findings are translated into pertinent, accessible messages for a country's citizens.
Using RGI data, NRGI found that in all but two African countries, there is an "implementation gap" between what laws say and how resource governance works in practice. This insight, and others uncovered in the report, have provided a further foundation from which to engage our stakeholders in Africa, and globally. The index data have been leveraged in research by the United Nations Economic Commission for Africa on the state of African governance, referenced in a Swiss federal council report on gold trade and human rights violations and used in a study on the state of extractives disclosures by the European Commission.
The RGI is also a resource for the private sector and investors. Morgan Stanley's emerging markets bonds group has begun using the RGI to assess potential investments in resource-rich countries. In at least one instance, the group used RGI findings to inform discussions with a resource-rich government in which they are considering investment. Further, multinational investment bank Citi reported using RGI data extensively for the bank's environmental, social and governance work.
MENA Resource Governance Course a Life-Changer
Bechara Youssef Samneh
Declining Investment Challenges Pro-Extraction Government in Colombia
Carlos Monge
El Contexto Externo y Nacional Pone en Cuestión la Estrategia Extractivista del Gobierno de Colombia
Debt Sustainability Challenges in Resource-Rich Developing Countries
David Mihalyi
From Contract Transparency to Open Contracting in Natural Resources
Gavin HaymanRobert Pitman

redpajama_set_name: RedPajamaCommonCrawl


large galvanized planter galvanized garden planters originally used to wash the clothes these galvanised dolly tubs now make simple large galvanized planters large galvanized planter box.
lesro industries reception seating waiting room wood lesro industries ct.
water mister home depot inspirational water mister nozzles home depot water bottle mister home depot.
dandelion killer home depot poison ivy and tough brush killer ready to use dandelion killer home depot canada dandelion killer spray home depot.
zapp electric guitar zapp electrical group zap electric bike.
grouting shower walls grouting a shower shower re grouting grouting shower wall tiles grouting shower wall tips grouting shower walls or floor first.
wayfair outdoor dining sets outdoor furniture patio coupon code basic prestigious 5 wayfair outdoor dining chairs wayfair outdoor wicker dining chairs.
chemlawn near me green nozzle for turf gun lesco chemlawn spray gun chemlawn spray gun nozzles.
charles thomas homes avenue photo 1 st charles robert thomas homes charles thomas homes elkhorn ne.
20 20 20 liquid fertilizer liquid fertilizer 20 20 20 liquid soluble fertilizer liquid all purpose fertilizer 20 20 20.