% spell-checker: disable %
\begin{pycode}[manager_appendix]
manager_appendix = texfigure.Manager(
    pytex, './',
    number=3,
    python_dir='python',
    fig_dir='figures',
    data_dir='data',
)
\end{pycode}
% spell-checker: enable %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                   Appendix 1                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Non-equilibrium Ion Populations}\label{nei}

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$,
\begin{equation}\label{eq:nei}
    \frac{\partial f_k}{\partial t} = n_e(R_{k+1}f_{k+1} + I_{k-1}f_{k-1} - I_kf_k - R_kf_k)
\end{equation}
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 \citep[see section 6 of][]{young_chianti_2016}. Setting the left hand side of \autoref{eq:nei} 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 \citet{masai_x-ray_1984,hughes_self-consistent_1985,shen_lagrangian_2015}, we can write \autoref{eq:nei} in matrix form,
\begin{equation}\label{eq:nei_mat}
    \frac{\partial}{\partial t}\mathbf{F} = \mathbf{A}\mathbf{F},
\end{equation}
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,
\begin{equation}\label{eq:rate_mat}
    \mathbf{A} = n_e
        \begin{pmatrix}
            -I_1 & R_2 & 0 & \dots & 0 \\
            I_1 & -(I_2 + R_2) & R_3 & \dots & 0 \\
             & \ddots & \ddots & &  \\
            \vdots & I_{k-1} & -(I_k + R_k) & R_{k+1} & \vdots \\
             & & \ddots & \ddots & \\
            0 & \dots & 0 & I_{Z} & -R_{Z+1} 
        \end{pmatrix}.
\end{equation}

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 \citep{macneice_numerical_1984,bradshaw_numerical_2009}. To solve \autoref{eq:nei_mat}, we use the ``deferred correction'' method of \citet{npl_modern_1961}, as pointed out by \citet{macneice_numerical_1984},
\begin{equation*}
    \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),
\end{equation*}
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 \autoref{eq:nei_mat} yields an expression for $\mathbf{F}_{j+1}$,
\begin{equation}\label{eq:nei_iterative}
    \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,
\end{equation}
where $\mathbb{I}$ is the identity matrix. To solve \autoref{eq:nei_mat}, we need only compute $\mathbf{A}_j$ for each $T(t_j)$, set $\mathbf{F}_0$ to the equilibrium ion populations, and iteratively compute \autoref{eq:nei_iterative} for all $j$. We solve \autoref{eq:nei_mat} for all elements in \autoref{tab:elements} and for each strand in the \AR{}.

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. \citet{macneice_numerical_1984} provide two rules for adaptively adjusting the time step to ensure $f_k$ changes sufficiently slowly. Additionally, \citet{shen_lagrangian_2015} provide an alternate scheme for choosing the time step \textit{a priori} based on the input electron density and the eigenvalues of \autoref{eq:rate_mat}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                   Appendix 2                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Effective AIA Response Functions}\label{effective_response_functions}

% spell-checker: disable %
\begin{pycode}[manager_appendix]
from synthesizAR.instruments import InstrumentSDOAIA
from scipy.interpolate import splev
aia = InstrumentSDOAIA([0,1]*u.s, observer_coordinate=None)
em = EmissionModel.restore(os.path.join(manager_appendix.data_dir, 'base_emission_model.json'))
fig,axes = plt.subplots(
    2, 3, sharex=True, sharey=True,
    figsize=texfigure.figsize(
        pytex, 
        scale=1,
        height_ratio=2/3,       
        figure_width_context='columnwidth'
    )
)
T = np.logspace(5,8,100)
p = 10**(15)*u.K/(u.cm**3)
const_p_indices = np.array([(i,np.argmin(np.fabs(em.density.value-d.value))) 
                            for i,d in enumerate(p/em.temperature)])
with h5py.File(os.path.join(manager_appendix.data_dir, 'effective_aia_response.h5'), 'r') as hf:
    colors = qualitative_palette(len([k for k in hf['94']])-1)
    for i, (ax, channel) in enumerate(zip(axes.flatten(), aia.channels)):
        grp = hf[channel['name']]
        real_response = splev(T, channel['temperature_response_spline'])
        ax.plot(T, real_response, ls='-',color='k',)
        ax.plot(em.temperature, 
                np.array(grp['response'])[const_p_indices[:,0],const_p_indices[:,1]],
                color='k',ls='--',)
        elements = sorted([e for e in grp if e != 'response'],
                            key=lambda x: plasmapy.atomic.atomic_number(x))
        for j,element in enumerate(elements):
            ax.plot(em.temperature, 
                    np.array(grp[element])[const_p_indices[:,0],const_p_indices[:,1]],
                    color=colors[j], ls='--', label=plasmapy.atomic.atomic_symbol(element),)
        if i==1:
            ax.legend(ncol=len(elements), loc="lower center", bbox_to_anchor=(0.5,1.02),
                      frameon=False)
        ax.text(6e7,1e-24,'{} $\mathrm{{\AA}}$'.format(channel['name']),
                fontsize=plt.rcParams['legend.fontsize'],
                horizontalalignment='right', verticalalignment='top')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_ylim([2e-30,2e-24])
ax.set_xlim([1e5,8e7])
ax.set_yticks([1e-29, 1e-27, 1e-25])
axes[1,0].set_ylabel(r'$K_c$ [DN cm$^5$ s$^{-1}$ pixel$^{-1}$]')
axes[1,0].set_xlabel(r'$T$ [K]')
plt.subplots_adjust(wspace=0.,hspace=0.)
fig_aia_response = manager_appendix.save_figure('aia-response')
fig_aia_response.caption = r'SSW temperature response functions (solid black) and effective temperature response functions for the elements in \autoref{tab:elements} (dashed black) for all six EUV AIA channels. The colored, dashed curves, as indicated in the legend, denote the contributions of the individual elements to the total response. For this calculation, we have assumed equilibrium ionization and a constant pressure of $10^{15}$ K cm$^{-3}$. We do not account for the time-varying degradation of the instrument.'
fig_aia_response.figure_env_name = 'figure*'
fig_aia_response.figure_width = r'\columnwidth'
fig_aia_response.placement = ''
fig_aia_response.fig_str = fig_str
\end{pycode}
\py[manager_appendix]|fig_aia_response|
% spell-checker: enable %

As discussed in \autoref{atomic}, 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 \autoref{eq:nei_mat} in place of the equilibrium ion population fractions. \autoref{fig:aia-response} shows the effective temperature response functions for the six EUV channels on AIA compared to those calculated from \texttt{aia\_get\_response.pro} in SSW. Even though we include a limited number of transitions from the CHIANTI database (see \autoref{tab:elements}), 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 \citep[see Section 2.1.6 of ][]{boerner_initial_2012}.

% spell-checker: disable %
\begin{pycode}[manager_appendix]
fig = plt.figure(figsize=texfigure.figsize(
    pytex,
    scale=0.5,
    height_ratio=1.0, 
    figure_width_context='columnwidth'
))
ax = fig.gca()
min_T = 1e300*u.K
max_T = 0*u.K
colors = heating_palette()
ar = synthesizAR.Field.restore(os.path.join(manager_appendix.data_dir, 'base_noaa1158'), lazy=True)
i_loop=680
loop = ar.loops[i_loop]
for i,h in enumerate(heating):
    loop.parameters_savefile = os.path.join(manager_appendix.data_dir, f'{h}', 'loop_parameters.h5')
    ax.plot(loop.electron_temperature[:,0], loop.density[:,0], color=colors[i], alpha=0.75,
            label=h.split('_')[0].capitalize())
    min_T = min(min_T, loop.electron_temperature.min())
    max_T = max(max_T, loop.electron_temperature.max())
p = 1e15*u.K/(u.cm**3)
T = np.linspace(min_T,max_T,1000)
ax.plot(T, p/T, color='k')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlim(1.75e5,1e7)
ax.set_ylim(1.1e7,6e9)
ax.legend(ncol=3,loc="lower center", bbox_to_anchor=(0.5, 1.02), frameon=False,)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['bottom'].set_bounds(2e5,1e7)
ax.spines['left'].set_bounds(2e7,6e9)
ax.set_ylabel(r'$n$ [cm$^{-3}$]')
ax.set_xlabel(r'$T$ [K]')
fig_hydro_phase_space = manager_appendix.save_figure('hydro-phase-space')
fig_hydro_phase_space.caption = r'$n-T$ phase-space orbits for a single strand for the three heating scenarios as defined by the legend. The black line indicates a constant pressure of $10^{15}$ K cm$^{-3}$.'
fig_hydro_phase_space.figure_width = r'0.5\columnwidth'
fig_hydro_phase_space.fig_str = fig_str
\end{pycode}
\py[manager_appendix]|fig_hydro_phase_space|
% spell-checker: enable %

Furthermore, during the evolution of a strand, the pressure is not constant for any of our heating scenarios as evidenced by \autoref{fig:hydro-phase-space}. 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 \autoref{eq:ppuv} at the correct temperature and density.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                   Appendix 3                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Computing Time Lags}\label{timelag_details}

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,
\begin{equation}\label{eq:cc_pre}
    \mathcal{C}_{AB}(\tau) = \mathcal{I}_A(t)\star\mathcal{I}_B(t) = \mathcal{I}_A(-t)\ast\mathcal{I}_B(t)
\end{equation}
where $\star$ and $\ast$ represent the correlation and convolution operators, respectively, $\tau$ is the lag and
\begin{equation*}
    \mathcal{I}_c(t)=\frac{I_c(t)-\bar{I}_c}{\sigma_{c}},
\end{equation*}
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 \autoref{eq:cc_pre}, using the convolution theorem, and then taking the inverse Fourier transform, $\mathscr{F}^{-1}$, gives,
\begin{align}\label{eq:cc}
    \fourier{\mathcal{C}_{AB}(\tau)} &= \fourier{\mathcal{I}_A(-t)\ast\mathcal{I}_B(t)}, \nonumber\\
    &= \fourier{\mathcal{I}_A(-t)}\fourier{\mathcal{I}_B(t)}, \nonumber\\
    \mathcal{C}_{AB}(\tau) &= \inversefourier{\fourier{\mathcal{I}_A(-t)}\fourier{\mathcal{I}_B(t)}}.
\end{align}
Scaling \autoref{eq:cc} by the length of the intensity time series $I(t)$ yields the same result as that of the correlation defined in Section 2 of \citet{viall_evidence_2012}. Furthermore, the \textit{time lag} between channels $A$ and $B$ is defined as,
\begin{equation}\label{eq:timelag}
    \tau_{AB} = \argmax_{\tau}\,\mathcal{C}_{AB}(\tau).
\end{equation}

By exploiting the definition of the cross-correlation as given in \autoref{eq:cc}, 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\times10^6$ times. We use the highly-optimized and thoroughly tested Fourier transform algorithms in the NumPy library for array computations \citep{oliphant_guide_2006} combined with the Dask library for parallel and distributed computing \citep{rocklin_dask_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 \AR{} in less than ten minutes. For comparison, doing the same calculation by calling the IDL function \texttt{c\_correlate.pro} on each pixel in serial would take $\approx14$ hours for all 15 channel pairs.
