Understanding Heating in Active Region Cores through Machine Learning I. Numerical Modeling and Predicted Observables

TL;DR

This paper develops a pipeline for forward modeling solar active region emissions to analyze how different heating frequencies affect observable diagnostics, providing insights into coronal heating mechanisms.

Contribution

It introduces a novel modeling pipeline that predicts emission diagnostics from active regions based on magnetic and hydrodynamic models, linking heating frequency to observable signatures.

Findings

Signatures of heating frequency persist in diagnostics.

Emission measure slope narrows with decreasing heating frequency.

Time lag becomes more spatially coherent at lower heating frequencies.

Abstract

To adequately constrain the frequency of energy deposition in active region cores in the solar corona, systematic comparisons between detailed models and observational data are needed. In this paper, we describe a pipeline for forward modeling active region emission using magnetic field extrapolations and field-aligned hydrodynamic models. We use this pipeline to predict time-dependent emission from active region NOAA 1158 as observed by SDO/AIA for low-, intermediate-, and high-frequency nanoflares. In each pixel of our predicted multi-wavelength, time-dependent images, we compute two commonly-used diagnostics: the emission measure slope and the time lag. We find that signatures of the heating frequency persist in both of these diagnostics. In particular, our results show that the distribution of emission measure slopes narrows and the mean decreases with decreasing heating frequency and that the range of emission measure slopes is consistent with past observational and modeling work. Furthermore, we find that the time lag becomes increasingly spatially coherent with decreasing heating frequency while the distribution of time lags across the whole active region becomes more broad with increasing heating frequency. In a follow up paper, we train a random forest classifier on these predicted diagnostics and use this model to classify real AIA observations of NOAA 1158 in terms of the underlying heating frequency.

Tables

Table 1. Table 1: All three heating models plus the two single-event control models. In the single-event models, the energy flux is not constrained by Equation 3 .

Name	$\epsilon$ (see Eq.1)	Energy Constrained?
high	0.1	yes
intermediate	1	yes
low	5	yes
cooling	1 event per strand	no
random	1 event per strand	no

Table 2. Table 2: Elements included in the calculation of Equation 5 . For each element, we include all ions for which CHIANTI provides sufficient data for computing the emissivity.

Element	Number of Ions	Number of Transitions
O	8	11892
Mg	11	31965
Si	13	30047
S	16	33091
Ca	17	42823
Fe	25	553541
Ni	19	83517

Table 3. Table 3: σ I / I ¯ subscript 𝜎 𝐼 ¯ 𝐼 \sigma_{I}/\bar{I} as defined by Equation 11 of Guarrasi et al. ( 2010 ) computed on a single image at t = 15 × 10 3 𝑡 15 superscript 10 3 t=15\times 10^{3} s for each channel and heating frequency. A larger value denotes a greater degree of contrast.

Channel [Å]	High	Intermediate	Low
94	3.06	4.63	4.19
131	5.56	3.61	6.13
171	2.79	2.81	3.25
193	2.69	2.80	2.79
211	2.73	2.83	2.84
335	2.63	3.08	3.21

Equations

$$\epsilon =\frac{⟨t_{\text{wait}}⟩}{{\tau }_{\text{cool}}}\{\begin{array}{ll}<1,& \text{high frequency},\\ \sim 1,& \text{intermediate frequency},\\ >1,& \text{low frequency}.\end{array}$$

$$F_{AR}=\frac{{\tau }_{\text{event}}}{2}\frac{{\sum }_l^{N_{\text{strands}}}{\sum }_i^{N_l}{q}_i{L}_l}{t_{\text{total}}}$$

$$\frac{|F_{AR}/N_{\text{strands}}-F_{\ast }|}{F_{\ast }}<\delta ,$$

$$P({\lambda }_{ij})=\frac{n_H}{n_e}\mathrm{Ab}(X)N_j(X,k)f_{X,k}{A}_{ji}\Delta E_{ji}{n}_e,$$

$$I_c=\frac{1}{4\pi }\sum _{\{ij\}}{\int }_{\text{LOS}}\mathrm{d}hP({\lambda }_{ij})R_c({\lambda }_{ij})$$

$$\frac{\partial f_k}{\partial t}=n_e(R_{k+1}{f}_{k+1}+I_{k-1}{f}_{k-1}-I_k{f}_k-R_k{f}_k)$$

$$\frac{\partial }{\partial t}\mathbf{F}=\mathbf{AF},$$

$$\mathbf{A}=n_e\left(\begin{array}{ccccc}-I_1& R_2& 0& \dots & 0\\ I_1& -(I_2+R_2)& R_3& \dots & 0\\ & \ddots & \ddots & & \\ ⋮& I_{k-1}& -(I_k+R_k)& R_{k+1}& ⋮\\ & & \ddots & \ddots & \\ 0& \dots & 0& I_Z& -R_{Z+1}\end{array}\right).$$

$${\mathbf{F}}_{j+1}={\mathbf{F}}_j+\frac{\Delta t}{2}\left(\frac{\partial }{\partial t}{\mathbf{F}}_{j+1}+\frac{\partial }{\partial t}{\mathbf{F}}_j\right)+\mathcal{O}(\Delta t^2),$$

$${\mathbf{F}}_{j+1}\approx {\left(\mathbb{I}-\frac{\Delta t}{2}{\mathbf{A}}_{j+1}\right)}^{-1}\left(\mathbb{I}+\frac{\Delta t}{2}{\mathbf{A}}_j\right){\mathbf{F}}_j,$$

$${\mathcal{C}}_{AB}(\tau )={\mathcal{I}}_A(t)\star {\mathcal{I}}_B(t)={\mathcal{I}}_A(-t)\ast {\mathcal{I}}_B(t)$$

$${\mathcal{I}}_c(t)=\frac{I_c(t)-{\bar{I}}_c}{{\sigma }_c},$$

$$\mathcal{F}\left\{{\mathcal{C}}_{AB}(\tau )\right\}$$

$$=\mathcal{F}\left\{{\mathcal{I}}_A(-t)\right\}\mathcal{F}\left\{{\mathcal{I}}_B(t)\right\},$$

$${\mathcal{C}}_{AB}(\tau )$$

$${\tau }_{AB}=\underset{\tau }{\mathrm{argmax}}{\mathcal{C}}_{AB}(\tau ).$$

Full text

Sun:corona, Sun:UV radiation, methods:numerical, hydrodynamics

††journal: The Astrophysical Journal††software: Astropy (The Astropy Collaboration et al., 2018), Dask(Rocklin, 2015), IPython(Pérez & Granger, 2007), matplotlib(Hunter, 2007), Numba(Lam et al., 2015), NumPy(Oliphant, 2006), PlasmaPy(PlasmaPy Community et al., 2018), PythonTeX(Poore, 2015), seaborn(Waskom et al., 2018), scipy(Jones et al., 2001), SolarSoftware(Freeland & Handy, 1998), SunPy(Mumford et al., 2018), yt(Turk et al., 2011)

1 Introduction

Nanoflares have long been used to explain the observed million-degree temperatures in the non-flaring solar corona. Due to their faint, transient nature, direct observations of nanoflares are made difficult by several factors, including inadequate spectral coverage of instruments, the efficiency of thermal conduction, and non-equilibrium ionization (Cargill, 1994; Winebarger et al., 2012; Barnes et al., 2016a). However, recent observations of “very hot” 8-10 MK plasma, the so-called “smoking gun” of nanoflares, have provided compelling evidence for their existence (e.g. Brosius et al., 2014; Caspi et al., 2015; Parenti et al., 2017; Ishikawa et al., 2017).

Critical to understanding the underlying heating mechanism is knowing whether the corona in non-flaring active regions is heated steadily or impulsively. However, because both waves and reconnection can lead to impulsive heating (Klimchuk, 2015), it is better to ask at what frequency do nanoflares repeat on a given magnetic strand. In the case of low-frequency nanoflares, the time between consecutive events on a strand is long relative to its characteristic cooling time, giving the strand time to fully cool and drain before it is re-energized. In the high-frequency scenario, the time between events is short relative to the cooling time such that the strand is not allowed to fully cool before being heated again. Steady heating may be regarded as nanoflare heating in the very high-frequency limit.

Before proceeding, we note that a magnetic strand, the fundamental unit of the low-$\beta$ corona, is a flux tube oriented parallel to the magnetic field that is isothermal in the direction perpendicular to magnetic field. We make the distinction that a coronal loop is an observationally-defined feature representing a field-aligned intensity enhancement relative to the surrounding diffuse emission, such that a single coronal loop may be composed of many thermally-isolated strands. Furthermore, we define the active region core as the area near the center of the active region whose X-ray and EUV emission is dominated by closed loops with both footpoints rooted in the photosphere .

In lieu of a direct observable signature of nanoflare heating, two parameters in particular have been used to diagnose the heating frequency in active region cores: the emission measure slope and the time lag. These diagnostics provide indirect signatures of the energy deposition via observations of the plasma cooling by thermal conduction, enthalpy, and radiation. We now discuss each of these observables in detail.

The emission measure distribution, $\mathrm{EM}(T)=\int\mathrm{d}h\,n_{e}^{2}$, where $n_{e}$ is the electron density and the integration is taken along the line of sight, is a useful diagnostic for parameterizing the frequency of energy deposition. Many observational and theoretical studies have suggested that the “cool” portion of the $\mathrm{EM}(T)$ (i.e. leftward of the peak, $10^{5.5}\lesssim T\lesssim 10^{6.5}$ K), can be described by $\mathrm{EM}(T)\sim T^{a}$ (Jordan, 1976; Cargill, 1994; Cargill & Klimchuk, 2004). The so-called emission measure slope, $a$, is an important diagnostic for assessing how often a single strand may be reheated and has been used by several researchers to interpret active region core observations in terms of both high- and low-frequency heating (see Table 3 of Bradshaw et al., 2012, and references therein). The “cool” emission measure slope typically falls in the range $2<a<5$, with shallower slopes indicative of low-frequency heating and steeper slopes associated with high-frequency heating. Many observational studies of active region cores have used the emission measure slope to make conclusions about the heating frequency (e.g. Tripathi et al., 2011; Warren et al., 2011; Winebarger et al., 2011; Schmelz & Pathak, 2012; Warren et al., 2012; Del Zanna et al., 2015).

To better understand observable properties of nanoflare heating, several researchers have used hydrodynamic models of coronal loops to examine how the emission measure slope varies with heating frequency (e.g. Mulu-Moore et al., 2011; Bradshaw et al., 2012; Reep et al., 2013). Most recently, Cargill (2014) found that varying the time between consecutive heating events from 250 s (high-frequency heating) to 5000 s (low-frequency heating) could account for the wide observed distribution of emission measure slopes, with higher values of $a$ corresponding to higher heating frequency due to the $\mathrm{EM}(T)$ distribution becoming increasingly isothermal (see also Barnes et al., 2016b).

In addition to the emission measure slope, the time lag analysis of Viall & Klimchuk (2012) has also been used by several workers to understand the frequency of energy release in active region cores. The time lag is the temporal delay which maximizes the cross-correlation between two time series, and, qualitatively, can be thought of as the amount of time which one signal must be shifted relative to another in order to achieve the best “match” between the two signals. As the plasma cools through the six EUV channels of the Atmospheric Imaging Assembly instrument (AIA, Lemen et al., 2012) onboard the Solar Dynamics Observatory spacecraft (SDO, Pesnell et al., 2012), we expect to see the intensity peak in successively cooler passbands of AIA according to the sensitivity of each channel in temperature space (Viall & Klimchuk, 2011). Computing the time lag between light curves in different channels provides a proxy for the cooling time between channels and insight into the thermal evolution of the plasma. Calculating the time lag in each pixel of an AIA image can reveal large scale cooling patterns in coronal loops as well as the diffuse emission between loops across an entire active region.

Viall & Klimchuk (2012) computed time lags for all possible AIA EUV channel pairs in every pixel of active region NOAA 11082 and found positive time lags across the entire active region core, indicative of cooling plasma. They interpreted these observations as being inconsistent with a steady heating model. Viall & Klimchuk (2017) extended this analysis to the 15 active regions catalogued by Warren et al. (2012) and found overwhelmingly positive time lags, or cooling plasma, in all cases, with only a few isolated instances of negative time lags, or heating plasma. These observations are consistent with an impulsive heating scenario in which little emission is produced during the heating phase because of the time needed to fill the corona by chromospheric evaporation and the efficiency of thermal conduction. Bradshaw & Viall (2016) predicted AIA intensities for a range of nanoflare heating frequencies in a model active region and applied the time lag analysis to their simulated images. They found that aspects of both high- and intermediate-frequency nanoflares reproduced the observed time lag patterns, but neither model could fully account for the observational constraints, suggestive of a range of heating frequencies across the active region. Additionally, Lionello et al. (2016) used a field-aligned hydrodynamic model to compute time lags for several loops in NOAA 11082 and concluded that an impulsive heating model could not account for the long ($>5000$ s) time lags calculated from observations by Viall & Klimchuk (2012).

Any successful heating model must be able to reproduce the observed distribution of emission measure slopes and time lags. In order to carry out such a test, both advanced forward modeling and sophisticated comparisons to data are required. In this paper, we carry out a series of nanoflare heating simulations in order to better understand how the frequency of impulsive heating events on a given strand is related to observable properties of the plasma, notably the emission measure slope and the time lag as derived from AIA observations. To do this, we use a combination of magnetic field extrapolations, hydrodynamic models, and atomic data to produce simulated AIA emission which can be treated in the same manner as real observations. We then apply the emission measure and time lag analyses to this simulated data. Section 2 provides a detailed description of both our forward modeling pipeline and the nanoflare heating model. In Section 3, we show the predicted intensities for each heating model and AIA channel (Section 3.1), the resulting emission measure slopes (Section 3.2.1) and the time lags (Section 3.3). Section 4 provides some discussion of our results and Section 5 includes a summary and concluding remarks.

This paper is the first in a series concerned with constraining nanoflare heating properties through forward modeled observables and serves to describe our forward modeling pipeline and lay out the results of our nanoflare simulations. In Barnes et al. (2019, Paper II hereafter), we use machine learning to make detailed comparisons to AIA observations of active region NOAA 1158. We train a random forest classifier using the predicted emission measure slopes and time lags presented here over the entire heating frequency parameter space in order to classify the heating frequency in each pixel of the observed active region. In contrast to past studies which have relied on a single diagnostic, this approach allows us to simultaneously account for an arbitrarily large number of observables in deciding which model fits the data “best.” The ability to quantitatively compare models with large quantities of data is crucial for progress in the current era where the amount of solar coronal data is orders of magnitude larger than in the past. Combined, these two papers demonstrate a novel method for using real and simulated observations to systematically predict heating properties in active region cores.

2 Modeling

In order to understand how signatures of the heating frequency are manifested in the emission measure slope and time lag, we predict the emission over the entire active region as observed by SDO/AIA for a range of nanoflare heating frequencies. To do this, we have constructed an advanced forward modeling pipeline through a combination of magnetic field extrapolations, field-aligned hydrodynamic simulations, and atomic data111Our forward modeling pipeline, called synthesizAR, is modular and flexible and written entirely in Python. The complete source code, along with installation instructions and documentation, are available here: github.com/wtbarnes/synthesizAR. In the following section, we discuss each step of our pipeline in detail.

2.1 Magnetic Field Extrapolation

We choose active region NOAA 1158, as observed by the Helioseismic Magnetic Imager (HMI, Scherrer et al., 2012) on 12 February 2011 15:32:42 UTC, from the list of active regions studied by Warren et al. (2012). The line-of-sight (LOS) magnetogram is shown in the left panel of Figure 1. We model the geometry of active region NOAA 1158 by computing the three-dimensional magnetic field using the oblique potential field extrapolation method of Schmidt (1964) as outlined in Sakurai (1982, Section 3). The extrapolation technique of Schmidt is well-suited for our purposes due to its simplicity and efficiency though we note it is only applicable on the scale of an active region. We include the oblique correction to account for the fact that the active region is off of disk-center.

The HMI LOS magnetogram provides the lower boundary condition of the vector magnetic field (i.e. $B_{z}(x,y,z=0)$) for our field extrapolation. We crop the magnetogram to an area of 300″-by-300″centered on $(-288.26\arcsec,-223.21\arcsec)$ and resample the image to 100-by-100 pixels to reduce the computational cost of the field extrapolation. Additionally, we define our extrapolated field to have a dimension of 100 pixels and spatial extent of $0.3R_{\sun}$ in the $z-$direction such that each component of our extrapolated vector magnetic field $\vec{B}$ has dimensions $(100,100,100)$.

After computing the three-dimensional vector field from the observed magnetogram, we trace $5\times 10^{3}$ field lines through the extrapolated volume using the streamline tracing functionality in the yt software package (Turk et al., 2011). Furthermore, we keep only closed field lines in the range $20<L<300$ Mm, where $L$ is the full length of the field line. The right panel of Figure 1 shows a subset of the traced field lines overlaid on the observed AIA 171 Å image of NOAA 1158. Contours from the observed HMI LOS magnetogram are shown in red (positive polarity) and blue (negative polarity). A qualitative comparison between the extrapolated field lines and the loops visible in the AIA 171 Å image reveals that the field extrapolation and line tracing adequately capture the three-dimensional geometry of the active region. Figure 2 shows the distribution of footpoint-to-footpoint lengths for all of the traced field lines.

2.2 Hydrodynamic Modeling

Due to the low-$\beta$ nature of the corona, we can treat each field line traced from the field extrapolation as a thermally-isolated strand. We use the Enthalpy-based Thermal Evolution of Loops model (EBTEL, Klimchuk et al., 2008; Cargill et al., 2012a, b), specifically the two-fluid version of EBTEL (Barnes et al., 2016a), to model the thermodynamic response of each strand. The two-fluid EBTEL code solves the time-dependent, two-fluid hydrodynamic equations spatially-integrated over the corona for the electron pressure and temperature, ion pressure and temperature, and density. The two-fluid EBTEL model accounts for radiative losses in both the transition region and corona, thermal conduction (including flux limiting), and binary Coulomb collisions between electrons and ions. The time-dependent heating input is configurable and can be deposited in the electrons and/or ions. A detailed description of the model and a complete derivation of the two-fluid EBTEL equations can be found in Appendix B of Barnes et al. (2016a).

For each of the $5\times 10^{3}$ strands, we run a separate instance of the two-fluid EBTEL code for $3\times 10^{4}$ s of simulation time to model the time-dependent, spatially-averaged coronal temperature and density. For each simulation, the loop length is determined from the field extrapolation. We include flux limiting in the heat flux calculation and use a flux limiter constant of 1 (see Equations 21 and 22 of Klimchuk et al., 2008). Additionally, we choose to deposit all of the energy into the electrons . To map the results back to the extrapolated field lines, we assign a single temperature and density to every point along the strand at each time step. Though EBTEL only computes spatially-averaged quantities in the corona, its efficiency allows us to calculate time-dependent solutions for many thousands of strands in a few minutes.

2.3 Heating Model

We parameterize the heating input in terms of discrete heating pulses on a single strand with triangular profiles of duration $\tau_{\textup{event}}=200$ s. For each event $i$, there are two parameters: the peak heating rate $q_{i}$ and the waiting time prior to the event $t_{\textup{wait},i}$. We define the waiting time such that $t_{\textup{wait},i}$ is the amount of time between when event $i-1$ ends and event $i$ begins. Following the approach of Cargill (2014), we relate the waiting time and the event energy such that $t_{\textup{wait},i}\propto q_{i}$. The physical motivation for this scaling is as follows. In the nanoflare model of Parker (1988), random convective motions continually stress the magnetic field rooted in the photosphere, leading to the buildup and eventual release of energy. If the field is stressed for a long amount of time without relaxation, large discontinuities will have time to develop in the field, leading to a dramatic release of energy. Conversely, if the field relaxes quickly, there is not enough time for the field to become sufficiently stressed and the resulting energy release will be relatively small.

In this work we explore three different heating scenarios: low-, intermediate-, and high-frequency nanoflares. We define the heating frequency in terms of the ratio between the fundamental cooling timescale due to thermal conduction and radiation, $\tau_{\textup{cool}}$, and the average waiting time of all events on a given strand, $\langle t_{\textup{wait}}\rangle$,

[TABLE]

We choose to parameterize the heating in terms of the cooling time rather than an absolute waiting time as $\tau_{\textup{cool}}\sim L$ (see appendix of Cargill, 2014). While a waiting time of 2000 s might correspond to low-frequency heating for a 20 Mm strand, it would correspond to high-frequency heating in the case of a 150 Mm strand. By parameterizing the heating in this way, we ensure that all strands in the active region are heated at the same frequency relative to their cooling time. Figure 3 shows the heating rate, electron temperature, and density as a function of time, for a single strand, for the three heating scenarios listed above.

For a single impulsive event $i$ with a triangular temporal profile of duration $\tau_{\textup{event}}$, the energy density is $E_{i}=\tau_{\textup{event}}q_{i}/2$. Summing over all events on all strands that comprise the active region gives the total energy flux injected into the active region,

[TABLE]

where $t_{\textup{total}}$ is the total simulation time, $N_{\textup{strands}}$ is the total number of strands comprising the active region, and $N_{l}=(t_{\textup{total}}+\langle t_{\textup{wait}}\rangle)/(\tau+\langle t_{\textup{wait}}\rangle)$ is the total number of events occurring on each strand over the whole simulation. Note that the number of events per strand is a function of both $\varepsilon$ and $\tau_{\textup{cool}}$.

For each heating frequency, we constrain the total flux into the active region to be $F_{\ast}=10^{7}$ erg cm*-2* s*-1* (Withbroe & Noyes, 1977) such that $F_{AR}$ must satisfy the condition,

[TABLE]

where $\delta\ll 1$. For each strand, we choose $N_{l}$ events each with energy $E_{i}$ from a power-law distribution with slope $-2.5$ and fix the upper bound of the distribution to be $\bar{B}_{l}^{2}/8\pi$, where $\bar{B}_{l}$ is the spatially-averaged field strength along the strand $l$ as derived from the field extrapolation. This is the maximum amount of energy made available by the field to heat the strand. We then iteratively adjust the lower bound on the power-law distribution for $E_{i}$ until we have satisfied Equation 3 within some numerical tolerance. We note that the set of $E_{i}$ we choose for each strand may not uniquely satisfy Equation 3.

We use the field strength derived from the potential field extrapolation to constrain the energy input to our hydrodynamic model for each strand. While the derived potential field is already in its lowest energy state and thus has no energy to give up, our goal here is only to understand how the distribution of field strength may be related to the properties of the heating. In this way, we use the potential field as a proxy for the non-potential component of the coronal field, with the understanding that we cannot make any quantitative conclusions regarding the amount of available energy or the stability of the field itself.

In addition to these three multi-event heating models, we also run two single-event control models. In both control models every strand in the active region is heated exactly once by an event with energy $\bar{B}_{l}^{2}/8\pi$. In our first control model, the start time of every event is $t=0$ s such that all strands are allowed to cool uninterrupted for $t_{\textup{total}}=10^{4}$ s. In the second control model, the start time of the event on each strand is chosen from a uniform distribution over the interval $[0,3\times 10^{4}]$ s, such that the heating is likely to be out of phase across all strands. In these two models, the energy has not been constrained according to Equation 3 and the total flux into the active region is $(\sum_{l}\bar{B}_{l}^{2}L_{l})/8\pi t_{\textup{total}}$. From here on, we will refer to these two models as the “cooling” and “random” models, respectively. All five heating scenarios are summarized in Table 1.

2.4 Forward Modeling

2.4.1 Atomic Physics

For an optically-thin, high-temperature, low-density plasma, the radiated power per unit volume, or emissivity, of a transition $\lambda_{ij}$ of an electron in ion $k$ of element $X$ is given by,

[TABLE]

where $N_{j}$ is the fractional energy level population of excited state $j$, $f_{X,k}$ is the fractional population of ion $k$, $\mathrm{Ab}(X)$ is the abundance of element $X$ relative to hydrogen, $n_{H}/n_{e}\approx 0.83$ is the ratio of hydrogen and electron number densities, $A_{ji}$ is the Einstein coefficient, and $\Delta E_{ji}=hc/\lambda_{ij}$ is the energy of the emitted photon (see Mason & Fossi, 1994; Del Zanna & Mason, 2018). To compute Equation 4, we use version 8.0.6 of the CHIANTI atomic database (Dere et al., 1997; Young et al., 2016). We use the abundances of Feldman et al. (1992) as provided by CHIANTI. For each atomic transition, $A_{ji}$ and $\lambda_{ji}$ can be looked up in the database. To find $N_{j}$, we solve the level-balance equations for ion $k$, including the relevant excitation and de-excitation processes as provided by CHIANTI (see Section 3.3 of Del Zanna & Mason, 2018).

The ion population fractions, $f_{X,k}$, provided by CHIANTI assume ionization equilibrium (i.e. the ionization and recombination rates are always in balance). However, in the rarefied solar corona, where the plasma is likely heated impulsively, it is not guaranteed that the ionization timescale is less than the heating timescale such that the ionization state may not be representative of the electron temperature (Bradshaw & Cargill, 2006; Reale & Orlando, 2008; Bradshaw, 2009). To properly account for this effect, we compute $f_{X,k}$ by solving the time-dependent ion population equations for each element using the ionization and recombination rates provided by CHIANTI. The details of this calculation are provided in Appendix A.

2.4.2 Instrument Effects

We combine Equation 4 with the wavelength response function of the instrument to model the intensity as it would be observed by AIA,

[TABLE]

where $I_{c}$ is the intensity for a given pixel in channel $c$, $P(\lambda_{ij})$ is the emissivity as given by Equation 4, $R_{c}$ is the wavelength response function of the instrument for channel $c$ (see Boerner et al., 2012), $\{ij\}$ is the set of all atomic transitions listed in Table 2, and the integration is along the LOS. Note that when computing the intensity in each channel of AIA, we do not rely on the temperature response functions computed by SolarSoft (SSW, Freeland & Handy, 1998) and instead use the wavelength response functions directly. This is because the response functions returned by aia_get_response.pro assume both ionization equilibrium and constant pressure. Appendix B provides further details on our motivation for recomputing the temperature response functions.

We compute the emissivity according to Equation 4 for all of the transitions in Table 2 using the temperatures and densities from from our hydrodynamic models for all $5\times 10^{3}$ strands. We then compute the LOS integral in Equation 5 by first converting the coordinates of each strand to a helioprojective (HPC) coordinate frame (see Thompson, 2006) using the coordinate transformation functionality in Astropy (The Astropy Collaboration et al., 2018) combined with the solar coordinate frames provided by SunPy (SunPy Community et al., 2015). This enables us to easily project our simulated active region along any arbitrary LOS simply by changing the location of the observer that defines the HPC frame. Here, our HPC frame is defined by an observer at the position of the SDO spacecraft on 12 February 2011 15:32:42 UTC (i.e. the time of the HMI observation of NOAA 1158 shown in Figure 1).

Next, we use these transformed coordinates to compute a weighted two-dimensional histogram, using the integrand of Equation 5 at each coordinate as the weights. We construct the histogram such that the bin widths are consistent with the spatial resolution of the instrument. For AIA, a single bin, representing a single pixel, has a width of 0.6″-per-pixel. Finally, we apply a gaussian filter to the resulting histogram to emulate the point spread function of the instrument. We do this for each time step, using a cadence of 10 s, and for each channel. For every heating scenario, this produces approximately $6(3\times 10^{4})/10\approx 2\times 10^{4}$ separate images.

3 Results

We forward model time-dependent AIA intensities using the method outlined in Section 2.4 for the heating scenarios discussed in Section 2.3. We discuss the predicted intensities in Section 3.1 for all six EUV channels of AIA and all five heating models. In Section 3.2 and Section 3.3, we show the results of the emission measure and time lag analyses, respectively, as applied to our simulated data. In Paper II, we use these simulated observables to train a machine learning classification model to understand with which heating scenario the real data are most consistent.

3.1 Intensities

We compute the intensities for the 94, 131, 171, 193, 211, and 335 Å channels of SDO/AIA using the procedure described in Section 2.4. We compute the intensity in each pixel of the model active region over a total simulation period of $3\times 10^{4}\,\mathrm{s}\,\approx 8.3$ hours with the exception of the cooling case which is only run for $10^{4}$ s. For the high-, intermediate-, low-frequency and “random” models, we discard the first and last $5\times 10^{3}$ s of evolution to avoid any transient effects in the strand evolution associated with the initial conditions and the constraints placed on the energy, respectively. We complete this procedure for each of the five heating scenarios in Table 1.

Figure 4 shows a snapshot of the intensity map at $t=15\times 10^{3}$ s for each channel and for the high-, intermediate-, and low-frequency nanoflare heating cases. The rows correspond to the different heating scenarios while the columns show the six AIA EUV channels. In each column, the intensities are normalized to the maximum intensity in the low-frequency case and are on a square-root scale. In general, we find that in the high-frequency intensity maps, individual loops are difficult to distinguish while in the low-frequency case individual loops appear bright relative to the surrounding emission. This distinguishability or “fuzziness” can be measured quantitatively as $\sigma_{I}/\bar{I}$, where $\sigma_{I}$ is the standard deviation taken over all pixels and $\bar{I}$ is the mean intensity (Guarrasi et al., 2010, Equation 11). A larger value of $\sigma_{I}/\bar{I}$ indicates a greater degree of contrast and vice versa. $\sigma_{I}/\bar{I}$ for each channel and heating frequency is shown in Table 3. With the exception of 131 Å, for every channel, the high-frequency case is the most “fuzzy”. The low-frequency case shows the most contrast in each channel except 94 Å though the margin between the low and intermediate cases is quite small in some cases.

Looking at the first two columns of Figure 4, we see that the intensity in the 94 and 131 Å channels increases as the heating frequency decreases. Both channels are double peaked and have “hot” ($\approx 7$ MK for 94 Å, $\approx 12$ MK for 131 Å) and “warm” ($\approx 1$ MK for 94 Å, $\approx 0.5$ MK for 131 Å) components. In the case of high-frequency heating, less energy is available per event such that few strands are heated to $>4$ MK. There is little emission in the 131 Å channel as strands are not often permitted to cool to $\leq 0.5$ MK either. However, in the low- and intermediate-frequency cases, we see several individual bright loops in both the 94 and 131 Å channels as the heating rate is sufficient to produce “hot” (i.e. 8-10 MK) loops. We see only a few of these loops as the lifetime of this hot plasma is short due to the efficiency of thermal conduction. In contrast, the faint, diffuse component of the 94 Å emission that is present in all three cases is due to the contribution of the “warm” component.

Additionally, we find that the 171 Å channel is dimmer for high-frequency heating as the peak sensitivity of this channel is $<1$ MK and in the case of high-frequency heating, strands are rarely allowed to cool below $1$ MK. In contrast, we note that the overall intensity in the 193, 211, and 335 Å channels is relatively constant over heating frequency as compared to the three previous channels though individual loops do become more visible with decreasing heating frequency. This relative insensitivity is because the temperature response functions of these three channels all peak in between 1.5 MK and 2.5 MK. In the case of high-frequency heating, strands are being sustained near these temperatures while in the low-frequency case, strands are cooling through this temperature range. This is illustrated for a single strand in Figure 3.

While there are clear differences in the AIA intensities between all three heating frequencies, quantifying these differences is difficult due in part to the multidimensional nature of the intensity data. To better understand how observational signatures differ as a function of heating frequency, we need to find a reduced representation of our data set that retains signatures of the underlying energy deposition. To this end, we compute two common observables: the emission measure slope (Section 3.2) and the time lag (Section 3.3).

3.2 Emission Measure Distributions

As discussed in Section 1, the emission measure slope is a useful quantity for understanding how frequently strands are re-energized. We compute emission measure distributions from our forward-modeled intensities using the regularized inversion method of Hannah & Kontar (2012). This method was designed to work with the narrowband coverage provided by AIA and so is well-suited to our needs. We choose our temperature bins such that the leftmost edge is at $10^{5.5}$ K and the rightmost edge at $10^{7.2}$ K with bin widths of $\Delta\log T=0.1$. Rather than computing $\mathrm{EM}(T)$ at each time step, we compute the time-averaged intensity in each pixel of each channel and compute $\mathrm{EM}(T)$ only once. We compute the uncertainties on the intensities in each channel using the aia_bp_estimate_error.pro procedure in SSW which incorporates uncertainties due to shot noise, read noise, dark subtraction, quantization, photometric calibration, and onboard compression.

3.2.1 Emission Measure Slopes

After computing $\mathrm{EM}(T)$ in each pixel using the regularized inversion procedure, we fit a first-order polynomial to the log-transformed emission measure and temperature bin centers, $\log_{10}\mathrm{EM}\sim a\log_{10}T$, to calculate the emission measure slope, $a$. In Paper II, we compare our modeled emission measure slopes to those derived from real AIA observations of NOAA 1158 using this same method.

Figure 6 shows the resulting emission measure slope, $a$, in each pixel of our forward-modeled active region for the high-, intermediate-, and low-frequency cases. We fit $\mathrm{EM}(T)$ over bins in the temperature range $8\times 10^{5}\,\textup{K}\leq T<T_{\textup{peak}}$. To assess the “goodness-of-fit” we use $r^{2}$, the correlation coefficient for the first-order polynomial fit, and mask pixels with $r^{2}<0.75$. Looking at the three panels in Figure 6, we find that overall, $a$ tends to decrease with decreasing frequency, consistent with previous modeling work (see Section 1). The low-frequency map (right panel) shows many values close to 2.

As the heating frequency increases, the slopes become larger, indicating an increasingly isothermal emission measure distribution. The intermediate-frequency map (center panel) shows predominantly higher slopes, with most pixels in the range $2\lesssim a\lesssim 3.5$ while the high-frequency map (left panel) shows much steeper slopes, with many $a\geq 3.5$, and a much broader range of slopes, $3\lesssim a\lesssim 8$. Note that in the high-frequency case, the slope varies considerably across the active region while the distribution of $a$ appears more spatially uniform in the intermediate- and low-frequency cases.

Below $T_{\textup{peak}}$, Cargill (1994) noted that the $\mathrm{EM}(T)$ could be described by $\mathrm{EM}(T)\sim n^{2}\tau_{rad}$, where $\tau_{rad}\sim T^{1-\alpha}n^{-1}$ is the radiative cooling time . We find that emission measure slopes produced by low-frequency nanoflares as shown in the right panel of Figure 6 are consistent with analytical results for single nanoflares .

3.2.2 Histograms

Figure 7 shows histograms of the emission measure slopes for the high-, intermediate-, and low-frequency cases. We find that the low-frequency distribution peaks at $a\approx 2.3$, inside the range expected from analytical results as noted above. The intermediate- and high-frequency distributions peak at successively higher values, $\approx 2.8$ and $\approx 4.0$, respectively. While the low- and intermediate-frequency distributions are more narrowly distributed around their peak values, the distribution of slopes in the high-frequency case is relatively broad and has a positive skew towards steeper slopes.

Looking at the distribution of slopes across the entire active region in Figure 7, we find that when the strands are heated infrequently such that each strand is allowed to cool fully prior to the next event, the distribution of slopes “saturates” in the range expected for single nanoflares. However, as the strands are reheated more often, the value of the slope becomes unsaturated and is subject to a range of cooling times due to the dependence of each waiting time on the power-law heating rate (see Section 2.3). These results are consistent with Cargill (2014) who computed $\mathrm{EM}(T)=n^{2}L$ for a single strand for a range of heating frequencies and found $a$ converged to $\approx 2$ for low frequency nanoflares and increased slowly with increasing heating frequency.

Our modeled emission measure slopes show that, even when accounting for the LOS integration, atomic physics, and information lost in the $\mathrm{EM}(T)$ inversion, signatures of the heating frequency still persist in the emission measure slope. However, while this quantity retains information about the frequency of energy deposition, drawing conclusions about the heating based solely on the observed emission measure slope, particularly for a small number of pixels may be misleading. As shown here and in Del Zanna et al. (2015), the slope may vary significantly across a given active region. Additionally, calculating $\mathrm{EM}(T)$ from observations is non-trivial due to several factors, including the mathematical difficulties of the ill-posed inversion (Craig & Brown, 1976; Judge et al., 1995, 1997), uncertainties in the atomic data (Guennou et al., 2013), and insufficient constraints from spectroscopic observations (e.g. Landi & Klimchuk, 2010; Winebarger et al., 2012), among others.

3.3 Time Lags

Next, we apply the time lag method of Viall & Klimchuk (2012) to our synthetic intensities for all of the heating scenarios discussed in Section 2.3. For each pixel in the active region, we compute the cross-correlation (Appendix C) for all pairs of the EUV channels of AIA (15 in total) and find the temporal offset which maximizes the cross-correlation according to Equation C3. The details of the cross-correlation and time lag calculation are given in Appendix C. Using the convention of Viall & Klimchuk (2012), we order the channel pairs such that the “hot” channel is listed first, meaning that a positive time lag corresponds to variability in the hotter channel followed by variability in the cooler channel. In other words, a positive time lag indicates cooling plasma. For the 94 Å and 131 Å channels, both of which have a bimodal temperature response function (see Figure 11), the order is determined by the component which is most dominant such that 94 Å is ordered first while 131 Å is ordered second. Thus, it is possible for cooling plasma to produce negative time lags in these channel pairs and the ambiguity can be resolved in the context of the time lags in other channel pairs.

3.3.1 Time Lag Maps

Figure 8 shows $\tau_{AB}$ (Equation C3) in each pixel of our simulated active region for all heating scenarios listed in Table 1 and three selected channel pairs: 94-335 Å, 211-131 Å, and 193-171 Å. Blacks, blues, and greens correspond to negative time lags; reds, oranges, and yellows correspond to positive time lags; and olive green indicates zero time lag. The range of the colorbar is $\pm 6000$ s. Note that the heating frequency decreases from left to right across each row. If the correlation in a given pixel is too low ($\max{\mathcal{C}_{AB}}<0.1$), the pixel is masked and colored white.

Looking at the first two rows of Figure 8, we find that the positive time lags in the 211-131 Å channel pair are significantly longer than those in the 94-335 Å pair. In the temperature range $2.5<T<7.3$ MK (94-335 Å), the dominant cooling mechanism is field-aligned thermal conduction while radiative cooling dominates in the range $0.6<T<2.5$ MK (211-131 Å). Because thermal conduction is far more efficient at high temperatures, the plasma spends less time in the $[T_{335},T_{94}]$ temperature range than in $[T_{131},T_{211}]$. The 193-171 Å time lags for the cooling case fall in the middle as radiative cooling also tends to dominate in this temperature interval ($0.9<T<1.5$ MK), but the separation in temperature space is smaller than the 211-131 Å pair. In all cases, these differences in the magnitude of the positive time lags become more apparent at lower heating frequencies.

In the 94-335 Å pair, we find negative time lags in the longest loops near the edge of the active region, inconsistent with our previous assertion that longer loops lead to longer, positive time lags. These loops are rooted in areas of weaker magnetic field (compared to the center) and thus do not have sufficient energy to evolve significantly into the temperature range of the “hot” component of the 94 Å channel (see Section 2.3). Thus, cooling from 335 Å to the cooler part of 94 Å dominates the cross-correlation. These negative time lags become more prominent as the heating frequency decreases. Consistent with the observations by Viall & Klimchuk (2017) of this same active region, our simulation results for the “cooling” scenario (rightmost column of Figure 8) show negative 94-335 Å time lags of $\approx-2000$ s in the longer loops on the lower edge of the active region. Viall & Klimchuk found far fewer positive 94-335 Å time lags .

We also find negative 211-131 Å time lags in the center of the active region for the high-, intermediate-, and low-frequency cases, indicative of plasma cooling from the hot part of the 131 Å channel through the 211 Å channel. Though we have not shown them here, similar negative time lag signatures are present in nearly all of the other 131 Å channel pairs as well. These results are consistent with that of Cadavid et al. (2014) who found that in inter-moss regions of active region NOAA 11250, intensity variations in the 131 Å channel preceded brightenings in all other EUV channels. In the two control cases, we do not find any negative time lags as the cross-correlations in the core are dominated by uninterrupted cooling from 211 Å to the cool part of 131 Å.

For the 193-171 Å channel pair, we find very few negative time lags because, unlike the 94 and 131 Å channels, the 193 and 171 Å channels are strongly peaked about a single temperature. Along with 211 Å, these channels which are strongly-peaked about a single temperature are important for disambiguating the signals in channels with a bimodal temperature response function (see Figure 11).

our simulated time lags show far fewer zero time lags than the observations of Viall & Klimchuk (2012, 2017) and the modeling work of Bradshaw & Viall (2016) due to the lack of transition region emission in our model. Transition region emission shows near-zero time lag because every layer (or temperature) of the transition region responds in unison. However, for the 193-171 Å channel pair, we find that zero time lags still dominate the inner core of the active region for all five heating scenarios, suggesting that the plasma is cooling into, but not through the 171 Å temperature bandpass (Viall & Klimchuk, 2017). This underscores the point that zero time lags do not imply steady heating (Viall & Klimchuk, 2015, 2016).

3.3.2 Cross-Correlation Maps

Figure 9 shows the peak cross-correlation value, $\max\mathcal{C}_{AB}$, for each selected channel pair. Looking first at all three channel pairs, we see that, on average, the cross-correlation increases as the heating frequency decreases. Additionally, we find that the highest cross-correlations tend to be in the center of the active region while the lowest tend to be on the outer edge. Comparing Figure 9 with the time lags in Figure 8 also reveals that negative time lags are correlated with lower peak cross-correlation values. Furthermore, other than the “cooling” scenario, we find that there are large variations from one loop to the next for all heating frequencies such that the spatial coherence of these peak cross-correlation values is low. In Paper II, we will use the peak cross-correlation value, in addition to the time lag, to classify the heating frequency in each observed pixel.

3.3.3 Histograms

Figure 10 shows histograms of time lags for every channel pair and all five heating scenarios. The time lags are binned between $-10^{4}$ s and $+10^{4}$ s in 60 s bins. Each histogram is colored according the corresponding heating scenario, consistent with Figure 3 and Figure 7. The columns are arranged such that heating frequency decreases from left to right. We show each channel pair for all heating models to demonstrate how the distribution of time lags evolves as the heating frequency varies.

Note that as the heating frequency decreases (from left to right), the number of negative time lags decreases. In the “cooling” case, there are very few negative time lags except for channel pairs which include one or both of the double-peaked channels (94 Å and 131 Å). For those channel pairs which include 94 Å and/or 131 Å, we expect to find negative time lags, even in the single-nanoflare cooling case as our convention of ordering the “hot” channel first has been violated such that cooling plasma can lead to negative time lags. For the remaining channel pairs, negative time lags are associated with the heating and cooling cycle being interrupted by repeated events on a given strand.

4 Discussion

For all of our heating models, we find negative time lags in at least one of the three channel pairs as shown in Figure 8. Negative time lags can be used to disambiguate the temperature sensitivity of the AIA passbands and can be produced in one of two ways: high-frequency heating in which the time lag is dominated by many frequent reheatings or a channel pair in which one channel is bimodal. While intensity in the 131 Å channel can correspond to either $<0.4$ MK or $>10$ MK plasma (see Figure 11), negative time lags in the 211-131 Å channel pair provide a possible signature of $\geq 10$ MK plasma because a negative time lag implies the plasma is cooling from 131 Å to 211 Å. This also holds for the 171-131 and 193-131 Å channel pairs as well while the 94-131 and 335-131 Å channels are more ambiguous due to the first channel in the pair being bimodal as well.

As noted in Section 3.3.1, the high-, intermediate-, and low-frequency maps for the 211-131 Å channel pair all show coherent negative time lags in the core. Because these strands are rooted in areas of strong field, enough energy is made available by the field (see Section 2.3) to heat them well into (and likely above) the hot component of the 131 Å passband. Since these strands are relatively short, the density increases rapidly enough for this hot plasma to be visible before it is washed out by thermal conduction.

Plasma undergoing pure cooling by radiation and thermal conduction produces a predictable and well-understood time lag signature. However, complicated heating scenarios and LOS geometries are likely to make it more difficult to interpret observed time lag signatures. Consider the case of a single cooling strand such that the maximum allowed time lag for a given channel pair $AB$ is the amount of time it takes to cool from $T_{A}$ to $T_{B}$ by thermal conduction and radiation. We may regard the “cooling” case in Figure 10 as the baseline time lag distribution given that all strands were heated only once at $t=0$ s. Because the time lag is primarily determined by the cooling phase of the strand, the time lag becomes primarily a function of the loop length $L$ since $\tau_{\textup{cool}}\propto L$. Two effects are likely to increase the decoherence of the baseline time lag distribution: multiple structures evolving out-of-phase along a given LOS (the “random” heating scenario) and multiple reheatings before the end of the cooling and draining cycle on a given strand. We note that multiple polluting structures along the LOS seem to primarily add negative time lags to the distribution (the “random” case) while increasing the frequency of events on a given strand extends the distribution in the positive direction. The latter effect also produces more negative time lags. Since steady heating can be thought of as nanoflare heating in the high-frequency limit ($\langle t_{\textup{wait}}\rangle\to 0$), we expect the distribution of time lags to approach a uniform distribution as the heating frequency increases (Viall & Klimchuk, 2016). steadily-heated loops have no preferred time lag.

While our model for the energy deposition (see Section 2.3) does not assume any specific physical heating mechanism, the parameterization of the heating frequency in Equation 1 has an interesting implication in the context of the Parker (1988) nanoflare model. Rearranging Equation 1 and recalling that $\tau_{\textup{cool}}\propto L$ gives $\langle t_{\textup{wait}}\rangle\propto L$, i.e. longer strands have longer absolute waiting times between heating events. Given that longer field lines tend to be rooted in regions of weaker magnetic field, this further implies that, where the field is stronger, energy is more quickly dissipated. According to Parker (1988), this dissipation is due to small-scale reconnection of flux tubes that are continually stressed by the convective motion of the underlying photosphere. Thus, in this context, our heating model implies that the reconnection and the underlying driver are more efficient in areas where the field is strongest.

Though we have not addressed it here, another possible mechanism for producing time-varying intensity in active regions is thermal non-equilibrium (TNE) wherein condensation cycles driven by highly-stratified, but steady footpoint heating lead to long-period intensity pulsations (Kuin & Martens, 1982). Though originally used to explain coronal rain (Antolin et al., 2010, 2015; Auchère et al., 2018) and prominences (Antiochos & Klimchuk, 1991), several workers (Mok et al., 2016; Winebarger et al., 2016; Froment et al., 2017; Winebarger et al., 2018; Froment et al., 2018) have recently claimed that TNE can produce time lag signatures similar to those of impulsive heating models, suggesting that observed time lags may be consistent with both impulsive and steady heating. However, it is not yet clear whether TNE is consistent with observed signatures of very hot (8-10 MK) plasma. Detailed comparisons between TNE and nanoflare simulations and observations are desperately needed.

5 Summary

We have carried out a series of numerical simulations in an effort to understand how signatures of the nanoflare heating frequency are manifested in two observables: the emission measure slope and the time lag. Additionally, we described each component of our pipeline for forward modeling active region emission. For a given magnetogram observation of the relevant active region (in this case, NOAA 1158), we compute a potential field extrapolation and trace a large number of field lines through the extrapolated vector field. For each traced field line, we run an EBTEL hydrodynamic simulation and use the resulting temperatures and densities, combined with data from CHIANTI and the instrument response function, to compute the time-dependent intensity. These intensities are then mapped back to the magnetic skeleton and integrated along the LOS in each pixel to create time-dependent images of the active region.

Using our novel and efficient forward modeling pipeline, we produced AIA images for all six EUV channels for $\approx 8$ hours of simulation time. From these results, we computed both the emission measure slope and the time lag for all possible channel pairs. We carried out these steps for three different nanoflare heating frequencies, high, intermediate, and low, (see Equation 1) in addition to two control models, for a total of five different heating scenarios (see Table 1).

Our results can be summarized in the following points:

As the heating frequency decreases, the emission measure slope, $a$, becomes increasingly shallow, saturating at $a\approx 2$. As the heating frequency increases, the distribution of slopes over the active region is shifted to higher values and broadens. 2. 2.

The time lag becomes increasingly spatially coherent with decreasing heating frequency. When strands are allowed to cool without being re-energized, the spatial distribution of time lags is largely determined by the distribution of loop lengths over the active region. 3. 3.

The distribution of time lags becomes increasingly broad as the heating frequency increases, consistent with the results of Viall & Klimchuk (2016). 4. 4.

Negative time lags in channel pairs where the second (“cool”) channel is 131 Å provide a possible diagnostic for $\geq 10$ MK plasma

In this paper, we have used our advanced forward modeling pipeline to systematically examine how the emission measure slope and time lag are affected by the nanoflare heating frequency. In Paper II, we use the model results presented here to train a random forest classifier and apply it to emission measure slopes and time lags derived from real AIA observations of NOAA 1158. The 15 channel pairs for the time lag and cross-correlation combined with the emission measure slope represent a 31-dimensional feature space and a single 500-by-500 pixel active region amounts to $2.5\times 10^{5}$ sample points. Performing an accurate assessment over this amount of data manually or “by eye” is at least impractical and likely impossible. Thus, the application of machine learning to the problem of assessing models in the context of real data is a critical step in understanding the underlying energy deposition in active region cores and, to our knowledge, has not yet been applied in this context.

CHIANTI is a collaborative project involving George Mason University (USA), the University of Michigan (USA), University of Cambridge (UK), and NASA Goddard Space Flight Center (USA). This research makes use of SunPy, an open-source and free community-developed solar data analysis package written in Python (SunPy Community et al., 2015) and PlasmaPy, a community-developed open source core Python package for plasma physics (PlasmaPy Community et al., 2018). SJB and WTB were supported by the NSF through CAREER award AGS-1450230. The work of NMV was supported by the NASA Supporting Research program. The complete source of this paper, including the data, code, and instructions for running the forward modeling code, can be found at github.com/rice-solar-physics/synthetic-observables-paper-models. SDO(AIA,HMI)

Appendix A Non-equilibrium Ion Populations

In order to account for ionization non-equilibrium, we compute $f_{X,k}$ as a function of time $t$ by solving the time-dependent ion population equations for each ion $k$ in each element $X$,

[TABLE]

where $n_{e}$ is the electron density and $R_{k}$ and $I_{k}$ are the temperature-dependent recombination and ionization rates of ion $k$, respectively. The ionization and recombination rates are computed using the data provided in CHIANTI. The ionization rates include both direct ionization and excitation autoionization and the recombination rates include both radiative and dielectronic recombination (see section 6 of Young et al., 2016). Setting the left hand side of Equation A1 to zero gives the equation of ionization equilibrium.

Note that for an element with atomic number $Z$, we must solve $Z+1$ coupled differential equations to find the non-equilibrium level populations. Following the approaches of Masai (1984); Hughes & Helfand (1985); Shen et al. (2015), we can write Equation A1 in matrix form,

[TABLE]

where $\mathbf{F}=(f_{1},f_{2},\ldots,f_{k},\ldots,f_{Z+1})$ and $\mathbf{A}$ is a ${Z+1\times Z+1}$ tridiagonal matrix containing the ionization and recombination rates, multiplied by the electron density,

[TABLE]

Due to drastic changes in the ionization and recombination rates with temperature, the above system of equations is very “stiff,” making explicit schemes extremely sensitive to the choice of time step (MacNeice et al., 1984; Bradshaw, 2009). To solve Equation A2, we use the “deferred correction” method of NPL (1961), as pointed out by MacNeice et al. (1984),

[TABLE]

where $\Delta t$ is the time step, $j$ indexes time, and $\mathcal{O}(\Delta t^{2})$ denotes terms of second order or higher in $\Delta t$. Dropping the higher-order terms and using Equation A2 yields an expression for $\mathbf{F}_{j+1}$,

[TABLE]

where $\mathbb{I}$ is the identity matrix. To solve Equation A2, we need only compute $\mathbf{A}_{j}$ for each $T(t_{j})$, set $\mathbf{F}_{0}$ to the equilibrium ion populations, and iteratively compute Equation A4 for all $j$. We solve Equation A2 for all elements in Table 2 and for each strand in the active region.

Though this method is unconditionally stable, $\Delta t$ should still be chosen carefully as $\mathbf{F}$ will relax to the equilibrium charge states for very long time steps. In general, smaller time steps should be chosen when the electron temperature varies rapidly and the electron density is high. In this paper, we exploit the adaptive time step provided by the two-fluid EBTEL code which accounts for changes in the temperature and electron density. MacNeice et al. (1984) provide two rules for adaptively adjusting the time step to ensure $f_{k}$ changes sufficiently slowly. Additionally, Shen et al. (2015) provide an alternate scheme for choosing the time step a priori based on the input electron density and the eigenvalues of Equation A3.

Appendix B Effective AIA Response Functions

As discussed in Section 2.4.1, the assumption of ionization equilibrium is likely to be violated in the impulsive heating cases considered here. Thus, we must recompute the contributions of each ion to the total channel response, using the result of Equation A2 in place of the equilibrium ion population fractions. Figure 11 shows the effective temperature response functions for the six EUV channels on AIA compared to those calculated from aia_get_response.pro in SSW. Even though we include a limited number of transitions from the CHIANTI database (see Table 2), we recover nearly all of the response from each channel. The high-temperature contributions in the SSW functions are due to continuum emission which we do not include in our model. In all cases, the continuum contribution is several orders of magnitude below peak of the channel response. Additionally, we do not account for the time variation in the wavelength response functions due to the degradation of the detector (see Section 2.1.6 of Boerner et al., 2012).

Furthermore, during the evolution of a strand, the pressure is not constant for any of our heating scenarios as evidenced by Figure 12. The black line of constant pressure $p=10^{15}$ K cm*-3* shows the default pressure at which the SSW AIA response functions are evaluated. The other lines show the temperature-density phase space evolution for the high-, intermediate-, and low-frequency cases for a single strand, none of which is well described by the assumption of constant pressure. By recomputing and interpolating the emissivity to the temperatures and densities as defined by our hydrodynamic simulation, we ensure that we are evaluating all quantities in Equation 4 at the correct temperature and density.

Appendix C Computing Time Lags

To find the associated time lag for a channel pair in a given pixel, we compute the cross-correlation between the time series associated with each channel and find the delay which maximizes this cross-correlation. We can express the cross-correlation $\mathcal{C}$ between two channels $A$ and $B$ as,

[TABLE]

where $\star$ and $\ast$ represent the correlation and convolution operators, respectively, $\tau$ is the lag and

[TABLE]

is the intensity of channel $c$ as a function of time, $I(t)$, centered to zero mean and scaled to unit standard deviation. Taking the Fourier transform, $\mathscr{F}$, of both sides of Equation C1, using the convolution theorem, and then taking the inverse Fourier transform, $\mathscr{F}^{-1}$, gives,

[TABLE]

Scaling Appendix C by the length of the intensity time series $I(t)$ yields the same result as that of the correlation defined in Section 2 of Viall & Klimchuk (2012). Furthermore, the time lag between channels $A$ and $B$ is defined as,

[TABLE]

By exploiting the definition of the cross-correlation as given in Appendix C, we can leverage existing Fourier transform algorithms in order to compute $\mathcal{C}_{AB}$ in a scalable and efficient manner. For a 500-by-500 pixel active region observation and 15 possible channel pairs, we need to compute $\tau_{AB}$ nearly $4\times 10^{6}$ times. We use the highly-optimized and thoroughly tested Fourier transform algorithms in the NumPy library for array computations (Oliphant, 2006) combined with the Dask library for parallel and distributed computing (Rocklin, 2015). Using Dask, we are able to parallelize this operation over many cores such that, on a 64-core machine, we can compute time lags for all 15 channel pairs in every pixel of the active region in less than ten minutes. For comparison, doing the same calculation by calling the IDL function c_correlate.pro on each pixel in serial would take $\approx 14$ hours for all 15 channel pairs.