1 Introduction

Understanding the basic physical properties of coronal loops is essential if we are to determine the dominant energy-transport mechanisms operating from the solar interior through the chromosphere and into the corona. Observed across all wavelengths, loops are the intrinsic building blocks of the magnetically closed solar atmosphere. These plasma loops range in temperature from a few 10,000 K up to 10 MK, can channel plasma flow up to 100 km s-1 (Fredvik et al., 2002), and show significant transverse oscillations following large impulsive events (Nakariakov et al., 2019). But do current observational capabilities indicate that we can resolve and definitively determine the inherent nature of these ubiquitous structures?

The launch of the Atmospheric Imaging Assembly (AIA: Lemen et al., 2012, 0.6′′ pixel) onboard NASA’s Solar Dynamics Observatory (SDO) has provided 24/7 observations of the whole Sun over a wide EUV wavelength range. Coupled to the Extreme Ultraviolet Imaging Spectrometer (EIS: Culhane et al., 2007) on Hinode, the resulting work on solar atmospheric loops is extensive (e.g. McIntosh and De Pontieu, 2009; Brooks, Warren, and Ugarte-Urra, 2012; Reale, 2014; Xie et al., 2017). Subsequently, the superior resolution images from NASA’s sounding rocket High-resolution Coronal imager (Hi-C: Kobayashi et al., 2014, 0.1′′ pixel) has led to observational assertions of the detection of magnetic braiding (Cirtain et al., 2013), nanoflares in inter-moss flares (Winebarger et al., 2013), and the determination of plasma loop parameters (Brooks et al., 2013; Peter et al., 2013; Scullion et al., 2014; Aschwanden and Peter, 2017; Williams et al., 2020b). It is still not fully settled whether the loop structures we observe with our current capabilities are actually spatially resolved or if they may be comprised of several individually isolated sub-resolution strands (Peter et al., 2013; Xie et al., 2017; Williams et al., 2020a).

When examining the make-up of coronal loops, one important factor for consideration is whether the observed loop structure is isothermal or multi-thermal along the line of sight. Differential emission measure (DEM) analysis of loop cross-sections have demonstrated both (Warren et al., 2008; Brooks, Warren, and Ugarte-Urra, 2012). Also consider the work by Schmelz, Christian, and Chastain (2016) where the authors create a coronal-loop inventory of 99 loop observations, finding 61 required multi-thermal loop models to reproduce the observations, another 28 being nearly or effectively isothermal, with another 10 firmly isothermal. Thus, when determining the temperature at a specific location in a loop (at the apex, say), if this yields a multi-thermal cross-section, then there is the argument that there could be either other plasma along the line of sight or sub-resolution elements existing in the loop that contribute a range of temperatures at that loop location. On the other hand, if there is a single dominant temperature recorded, then either the structure at that point is resolved (is monolithic) and hence at that temperature or if there are strands, their temperatures at the observed location are all very similar – they appear coherently. Determining a resolved fundamental spatial scale or the presence of sub-elements within a loop is an important step in addressing, for example, how coronal plasma is possibly being heated.

One of the favoured approaches attempting to explain the heating of coronal loop structures is the nanoflare model (Parker, 1983, 1988). Here coronal plasma is heated to million degree temperatures due to the occurrence of numerous, small-scale, localised bursts of energy somewhere within the overall loop envelope. The energy release could be between \(10^{23}\) – \(10^{27}\) erg, with possible differing mechanisms of the underlying physics being magnetic reconnection (Parker, 1988) or magnetohydrodynamic wave dissipation (Moriyasu et al., 2004). Antolin et al. (2008) demonstrate a link between the power-law index and the operating heating mechanism in a loop using two 1.5D models; one of Alfvén wave dissipation and one of random energy bursts to simulate nanoflares. From this they obtain different coronae for the two models, which can be seen in simulated Hinode/XRT flow patterns. It is expected the dissipation of Alfvén waves would lead to frequent, short-lived heating events along individual magnetic-field lines (Asgari-Targhi and van Ballegooijen, 2012). However, Antolin et al. (2021) have recently observed the presence of “nanojets” in coronal loops. These “nanojets” are postulated to be a signature of coronal heating caused by localised magnetic-reconnection activity and the results strongly suggests nanoflare storms arise due to curved magnetic-field lines reconnecting at small angles. Investigating the heating frequency of individual strands in the corona could further reveal the nature of nanoflare events, that is whether they are caused exclusively by small-scale magnetic reconnection events or whether they may also be the result of Alfvén waves.

Additionally, observations from different wavelength ranges provide strong constraints on the heating mechanism of loop structures (Cargill, 1993, 1994; Cargill and Klimchuk, 1997, 2004). Following this, Mendoza-Briceño, Erdélyi, and Di G. Sigalotti (2002) present numerical calculations detailing the response of coronal plasma to periodic, random micro-scale heating pulses within a magnetic loop. The results show successive energy bursts can maintain the coronal plasma to typical coronal temperatures, as well as providing good qualitative agreement with Transition Region and Coronal Explorer (TRACE) loops observations (Aschwanden, Schrijver, and Alexander, 2001). More recently Cargill, Warren, and Bradshaw (2015) demonstrated that the time between nanoflare events on individual magnetic strands should be between 500 – 2000 seconds, where the energy release may be as small as a few \(10^{23}\) erg.

One numerical approach that provides a quick and accurate answer to the coronal response of a loop to heating is the use of zero-dimensional (0D) field-aligned hydrodynamical models (EBTEL: Klimchuk, Patsourakos, and Cargill, 2008; Cargill, Bradshaw, and Klimchuk, 2012a,b; Cargill, Warren, and Bradshaw, 2015). As each loop is modelled as a single grid point, 0D models allow the coronal-loop evolution to be simulated very efficiently with minimal computational cost, whilst still being capable of obtaining time-dependent, spatially averaged loop quantities that are comparable to fully resolved 1D simulations. Consequently, this method for simulating the evolution of coronal loops has proved popular (Viall and Klimchuk, 2011; Cargill, Warren, and Bradshaw, 2015; Barnes, Cargill, and Bradshaw, 2016).

However, 1D simulations are still greatly useful, particularly those which model a multi-stranded structure as opposed to monolithic loops. For example, Price and Taroyan (2015) forward model the synthetic profiles of a well-documented Hinode/EIS structure (McIntosh and De Pontieu, 2009) using various numbers of sub-element strands to try and deduce the nature of the outflows observed. Their results support a scenario whereby long loops formed of multiple strands undergo periodic heating and cooling. Similarly, Susino et al. (2013) demonstrate that whilst 1D multi-stranded models can predict observed values derived from TRACE and AIA filter ratios, they cannot explain all of the characteristics of warm over-dense loops as such models assume spatial interactions between individual strands can be neglected, which may not always be the case (Hood et al., 2016; Reid et al., 2018, 2020).

Building upon and improving the computational approach of the multi-stranded loop code from Sarkar and Walsh (2008, 2009), this article uses synthesised emission as would be observed by AIA to investigate whether over a long time series of data (six hours) it is possible to recover an imposed heating frequency as well as any information about potential coronal loop sub-elements that are beyond the spatial resolving power of the instrument – such as the recent high-resolution strand observations by Williams et al. (2020a,b). Section 2 details the updated numerical approach that now incorporates the unresolved transition region (UTR) method developed by Johnston et al. (2017a,b). In Section 3, an analysis of the loop-apex temperature evolution from 18 different loop configurations is detailed. The data are employed to recreate synthetic AIA emission in Section 4 where the resulting light curves are examined using methods such as emission-line ratio comparisons and time-lag analysis (Viall and Klimchuk, 2011). The findings of this study and their implications on the determination of heating rate and strand number are discussed in Section 5.

2 Improved Multi-Stranded Loop Model

The numerical code employed in this article is written in FORTRAN and is based on the Lagrangian remap method (LareXd) developed by Arber et al. (2001), which has been modified for multi-stranded coronal loops (Sarkar and Walsh, 2008, 2009). A unique feature of this code is that each individual sub-element strand of a loop is an individual, independent simulation, all of which are then amalgamated in post-processing to form a single coherent coronal loop. The code is managed in such a way that the loop/strand parameters (such as nanoflare energy and loop/strand length and radius) can be edited without the need for recompiling the model, allowing for multiple instances/scenarios to run concurrently. In this article, we continue the development of the multi-stranded hydrodynamic loop (MSHDL) simulation by treating the transition region as an unresolved discontinuity across which energy is conserved by imposing a jump condition: the Unresolved Transition Region (UTR) method from Johnston et al. (2017a,b).

As will be discussed in this section, adopting the UTR method means the resolution of the MSHDL simulations can be reduced with minimal impact on the “observables” compared to other loop models (Johnston et al., 2017a,b), whilst also decreasing the computational time of each strand; an important factor when each loop is comprised of numerous strands. Similar approaches have been also been adopted in previous studies (for example, see Mikić et al., 2013; Johnston et al., 2019; Johnston and Bradshaw, 2019; Van Damme et al., 2020).

2.1 Initial Model Setup

The plasma-\(\beta \) is the ratio of thermal pressure to magnetic pressure within a given plasma. In the corona, and structures such as coronal loops, the \(\beta < 1\). This means the plasma is “tied-in” to the magnetic field, and subsequently, the shape, size, and orientation of the plasma structures are determined by the field lines. In addition to a low plasma-\(\beta \), the corona is also a highly conducting medium. These properties allow for the assumption that the plasma dynamics occur along magnetic-field lines, with negligible feedback between the magnetic-field and thermodynamic parameters. As is discussed by Sarkar and Walsh (2008, 2009), these physical attributes support the use of a 1D hydrodynamic code, which neglects thermal conduction across field lines when modelling entire loop-like structures.

With this in mind, here a loop is modelled as an amalgamation of a collection of individual, narrower, sub-loop element, which here are called “strands”. Each strand is assumed to have the same cross-sectional radius and total length. The MSHDL code is a 1D hydrodynamic adaptation of LareXd (Arber et al., 2001) that simulates each of these strands individually, which are then combined in post-processing to form a single loop structure. These loops can then be compared with observed loop structures that may appear as monolithic in nature but may actually comprise several to hundreds of structures whose spatial scales are below current instrument capabilities.

As aforementioned, the modelled loop is a compilation of single strands [\(i\)], each obeying the standard equations for mass, momentum, and energy conservation in curvilinear abscissa along a semi-circular loop:

$$\begin{aligned} \frac{\mathrm{D}\rho _{i}}{\mathrm{D}t} &=-\rho _{i} \frac{\partial v_{i}}{\partial s}, \end{aligned}$$
(1)
$$\begin{aligned} \rho _{i} \frac{\mathrm{D}v_{i}}{\mathrm{D}t}& = - \frac{\partial p_{i}}{\partial s} + \rho _{i} g + \rho _{i} v_{i} \frac{\partial ^{2}v_{i}}{\partial s^{2}}, \end{aligned}$$
(2)
$$\begin{aligned} \frac{\rho _{i}^{\gamma }}{\gamma -1}\frac{\mathrm{D}}{\mathrm{D}t} \left (\frac{p_{i}}{\rho _{i}^{\gamma }}\right )& = \frac{\partial }{\partial s}\left (\kappa \frac{\partial T_{i}}{\partial s}\right ) - n_{i}^{2}Q\left (T_{i} \right ) + H_{i}\left (s,t\right ), \end{aligned}$$
(3)

where

$$ p_{i} = \frac{\mathrm{R}}{\tilde{\mu }}\rho _{i} T_{i} $$
(4)

and

$$ \frac{\mathrm{D}}{\mathrm{D}t} \equiv \frac{\partial }{\partial t} + v_{i} \frac{\partial }{\partial s}. $$
(5)

Following standard notation, \(\rho _{i}\), \(p_{i}\), and \(T_{i}\), are the density, thermal pressure, and temperature of strand \(i\). The velocity [\(v_{i}\)] is along the curvilinear abscissa [\(s\)] of the strand. The adiabatic index, and mean molecular mass are given by and \(\tilde{\mu }=0.6~\textrm{mol}^{-1}\), respectively. Gravity is assumed to be constant and is taken to be the surface gravity: \(g = 2.74 \times 10^{4}~\textrm{cm}\,\text{s}^{-2}\). The number density [\(n_{i}\)] is defined as \(10^{21}\) [c\(\textrm{m}^{-3}\)] \(\rho _{i}\) [kg \(\textrm{m}^{-3}\)]. The conductivity of the plasma in the direction of \(s\) is denoted by \(\kappa = 9.2\times 10^{-7}T^{5/2}~\textrm{erg}\,\text{s}^{-1}\,\textrm{cm}^{-1}\,\textrm{K}^{-1}\). \(\mathrm{R} = 8.3 \times 10^{7}~\textrm{erg}\,\text{mol}^{-1}\,\textrm{K}^{-1}\) is the molecular gas constant, and \(Q(T_{i})\) is the optically thin radiative loss function (Rosner, Tucker, and Vaiana, 1978). The remaining term [\(H_{i}\left (s,t \right )\)] is our heating input for the individual strands, which is discussed briefly in Section 2.1.1 and presented in more detail by Sarkar and Walsh (2008, 2009).

For the simulated loop, a coronal length of \(L=100\) Mm is adopted. The loop foot-points are at and , with the loop apex situated at \(s = 0\). The loop radius is 1 Mm giving the loop an aspect ratio of 50:1, which is consistent with high-resolution observations for loops longer than 50 Mm (Peter et al., 2013). The area of each strand is approximated as \(A_{\mathrm{strand}} = A_{\mathrm{loop}}/N_{\mathrm{strands}}\), where \(A\) is cross-sectional area, and \(N_{\mathrm{strands}}\) is the total number of strands forming the loop.

The boundary conditions are prescribed as

$$ T\left (-L/2,t\right ) = T\left (+L/2,t\right ) = T_{\mathrm{ch}} = 10^{4}~\textrm{K} $$
(6)

and

$$ p\left (-L/2,t\right ) = p\left (+L/2,t\right ) = p_{\mathrm{ch}} = 0.315~\textrm{Pa}, $$
(7)

at the loop footpoints at the base of the chromosphere, which is 5 Mm deep. Note that \(T_{\mathrm{ch}}\) and \(p_{\mathrm{ch}}\) are the chromospheric temperature and pressure.

To define the temperature of the loop from the amalgamated strands, the emission-measure temperature is calculated, as is described by Sarkar and Walsh (2009), by

$$ T = \frac{\sum \limits _{i} \rho _{i}^{2}\left (s,t\right ) T_{i}\left (s,t\right )}{\sum \limits _{i} \rho _{i}^{2}\left (s,t\right )} . $$
(8)

2.1.1 Mimicking Nanoflare Heating

The heating is prescribed on each strand as a series of small, localised, energy-release episodes akin to the nature of nanoflares. The nature of these small energy bursts are governed by a random-number generator, such that their spatial and temporal distributions are unique and distinct for each strand of the loop. This heating approach is described fully by Sarkar and Walsh (2009), and so they are just restated briefly here.

Each heating episode has a minimum energy threshold, which may be set by the user – although a typical energy would be \(10^{24}~\textrm{erg}\). During the simulation, each strand is subject to a user-specified number of energy episodes whose spatial deposition is random but is weighted such that the energy release predominantly occurs near the base of the strands between loop positions ±35 – 39 Mm (Figure 1). The length scale of each heating event is fixed at 2 Mm whilst the duration is randomly distributed between 50 – 150 seconds.

Figure 1
figure 1

The top (bottom) panel shows the spatio-temporal plot indicating the duration, position, and magnitude of each energy episode within a loop consisting of 64 (8) strands, which is subjected to localised heating events every 250 seconds (1500 seconds) with minimum flare energy \(2.5\times 10^{23}\) erg (\(1.2\times 10^{25}\) erg). The magnitudes of the energy episodes are shown on the same colour table, whilst the blue lines indicate the nanoflare number densities as a function of loop position for the two scenarios. The total energy deposit of each loop is equivalent and hence together, the two panels provide a comparison of the heating distribution between a loop that is heated by many small-magnitude energy episodes (top) and a loop that is subjected to a few energy episodes that have a large magnitude (bottom).

Full size image

2.2 Modelling the Transition Region

In this article we consider the transition region to be a thin boundary layer between the chromosphere and corona where there is a rapid increase (decrease) in temperature (density). In coronal-loop simulations, the ability to accurately resolve the behaviour of this boundary layer is important. For example, Bradshaw and Cargill (2013) demonstrate that inadequately reproducing the transition region leads to artificially low coronal densities. In this article, the approach outlined by Johnston et al. (2017a,b) is adopted in the MSHDL code by treating the lower transition region as a discontinuity that responds to changing conditions – such as mass flow and plasma heating – through the imposition of a jump condition derived from an integrated form of energy conservation. This method has been shown to capture accurately the coronal-density evolution of fully resolved 1D models (e.g. Bradshaw and Cargill, 2013), whilst also being computationally significantly more efficient.

Here, the key points for modelling the transition region using the UTR method are briefly outlined, as per Sections 2 and 3 of Johnston et al. (2017a,b). The base of the UTR is defined as the location where the temperature first reaches, or drops below the chromospheric temperature (\(10^{4}\) K) when moving from the strand apex towards the chromosphere along the curvilinear abscissa [\(s\)]. Below the base of the UTR, the temperature is held fixed at the chromospheric temperature – \(10^{4}\) K.

The top of the UTR is defined as the final location at which the following criterion is satisfied when travelling away from the apex along \(s\):

$$ \frac{L_{R}}{L_{T}} \leq \delta < 1, $$
(9)

where \(L_{R}\) and \(L_{T}\) are the spatial and temperature length scales of the simulation as defined by Johnston et al. (2017a). In order to ensure that multiple gridpoints are maintained across the temperature length scale at this region, \(\delta =0.15\) is selected. These two definitions allow for easy identification of the UTR across all time steps of the simulation.

Having identified the UTR, the second part of the UTR method is the imposition of the jump condition at the top of the UTR, which takes the form of a local velocity correction which conserves the total energy in the UTR and is imposed at each time step. The jump condition is reproduced here from Equation 12 of Johnston et al. (2017a)

$$ \frac{\gamma }{\gamma -1}p_{0}v_{0} + \frac{1}{2}\rho _{0}v_{0}^{3} + \rho _{0}\Phi _{0}v_{0} = l_{\mathrm{UTR}}\bar{H} - (R_{\mathrm{UTR}}+F_{c,0}), $$
(10)

where the left-hand terms denote the flux for enthalpy, kinetic energy, and gravitational potential energy, respectively. The right-hand terms describe the averaged volumetric heating rate per unit cross-sectional area, the integrated radiative losses in the UTR [\(R_{\mathrm{UTR}}\)] and heat flux [\(F_{c}\)]. \(\Phi \) is the gravitational potential, \(\bar{H}\) is the spatially averaged volumetric heating rate, and \(l_{\mathrm{UTR}}\) is the length of the UTR. The subscript 0 in Equation 10 denotes quantities measured at the top of the UTR. The full details on the UTR method are given by Johnston et al. (2017a).

2.2.1 Comparison of Two Numerical Approaches

Here, a comparison is made between the two versions of the multi-stranded loops code: MSHDL (Sarkar and Walsh, 2008, 2009) and the updated code MSLUTR (Multi-Stranded Loops with Unresolved Transition Regions) employing the established UTR jump condition (Johnston et al., 2017a,b). For this comparison, a single localised release of energy is modelled at position \(s=-23\) Mm depositing a total energy of \(2.3\times 10^{24}\) erg over a length of 2 Mm for a duration of 1218 seconds. The temporal evolution of this energy burst near the \(-L/2\) footpoint is shown for the two codes at 20 second intervals in Figure 2.

Figure 2
figure 2

Single-strand velocity profile comparison of the previous MSHDL code (black) vs. new MSLUTR method (blue) for a single burst event near the \(-L/2\) footpoint in 20 second intervals. The base of the transition region is denoted by the vertical dashed-red line. Only one footpoint is shown so that the differences between the models at the transition region can be seen more clearly.

Full size image

The deposited energy drives a plasma flow that propagates down the strand leg until the conduction front reaches the transition region. Initially, the discrepancies between the MSHDL and MSLUTR are negligible. This can be attributed to the dynamics being set by the direct in-situ heating for the first 60 seconds. However, once the enhanced downward heat flux reaches the transition region at 60 seconds notable differences start to arise. This is the start of the evaporation phase, and by 80 seconds the leading front of the flow in MSLUTR is \(\approx 1\) Mm ahead of the corresponding flow in MSHDL.

By 100 seconds it can be seen that MSHDL is underestimating the evaporative upflow compared to MSLUTR. The largest difference in the velocity profiles between the two regimes occurs at \(\approx -47\) Mm; where MSLUTR velocity \(\approx 85\) km s-1 vs. \(\approx 50\) km s-1 for MSHDL.

In addition to MSHDL underestimating the coronal velocity, under-resolving the lower transition region also results in spurious oscillations propagating upwards from the chromosphere into the corona. These oscillations can be seen clearly at times t = 80, 100, 120, and 140 seconds. In contrast, these oscillations are significantly less prominent for MSLUTR but a small perturbation can still be seen to the right of the base of the transition region between 80 seconds – 140 seconds (red line, Figure 2). However, unlike MSHDL, the perturbation is contained within the UTR and does not propagate upwards into the coronal part of the structure. Furthermore, as is evidenced in previous work, adopting the UTR method provides comparative results with other fully resolved 1D simulations of impulsive heating (Johnston et al., 2017a,b) and the development of thermal non-equilibrium (Johnston et al., 2019).

Whilst Figure 2 highlights the improvements to modelling a single strand, one of the main benefits to using this code over other numerical schemes is the fact that it investigates the multi-stranded nature of loop structures. Thus, it is not enough for improved accuracy on a single strand during a single impulsive event if the run-time and/or spatial resolution required to achieve this greater accuracy makes multi-strand analysis impractical. To this end, a maximum velocity convergence test is performed using the aforementioned scenario where the number of gridpoints are varied using the MSLUTR code (Table 1). From this convergence test, it is deduced that a simulation resolution of nx = 2000 is sufficiently accurate with maximum coronal velocities within \(\approx 4.5\)% of an over-resolved simulation (nx = 15,000). This choice allows for improved accuracy compared to the MSHDL code whilst allowing the run-time to remain sufficiently manageable for potentially hundreds of strands to be simulated in a timely manner.

Table 1 MSLUTR velocity convergence test relative to \(nx = 15{,}000\) simulation.
Full size table

3 Multi-Stranded Analysis

This section provides results from the improved version of the multi-stranded loop code MSLUTR. The model loop is configured to be 100 Mm long with a radius of 1 Mm and initial chromospheric and coronal temperatures of 104 K and 106 K, respectively. The number of strands and energy-release episodes (i.e. nanoflares) the loop is subjected to varies in such a way that each loop configuration has approximately the same total amount of energy deposited (mean: \(1.10\times 10^{28}\) erg, \(\sigma \): \(0.04\times 10^{28}\) erg) during the simulation run, i.e. the greater the number of nanoflares and/or strands within a loop, the smaller the magnitude of energy release per nanoflare that occurs. Since the total energy deposition is approximately the same across all the simulations, direct comparisons can be made between the different heating frequencies investigated. The loop is simulated for a total time of 7.5 hours, but the first 1.5 hours are ignored. The reasoning behind this is that each strand requires several acoustic timescales to elapse due to i) the sudden initial injection of energy generating large-scale perturbations to damp throughout the system and ii) several heating events needing to occur to allow a quasi-steady equilibrium to be reached. Ignoring the first 1.5 hours is an arbitrary value that is sufficiently long to account for these hence allowing 6 hours to be analysed from each generated data-set.

In the literature regarding nanoflare heating frequencies, it is common to use the repetition time when discussing frequency due to the units being more meaningful, and as such we adopt this nomenclature for the remainder of this study (Table 2 also lists the nanoflare frequency in Hz for completeness). For example, Cargill, Warren, and Bradshaw (2015) argue that the repetition time of nanoflares along an individual strand to be in the range of 500 – 2000 seconds. With this in mind, nanoflare repetition times are explored that correspond to low-frequency heating (\(\approx 1500\) seconds: termed LFH), intermediate–low-frequency heating (\(\approx 1250\) seconds: ILFH) intermediate frequency heating (\(\approx 1000\) seconds: IFH), intermediate-high frequency heating (\(\approx 1500\) seconds: IHFH), high frequency heating (\(\approx 500\) seconds: HFH), and ultra-high frequency heating (\(\approx 250\) seconds: UHFH). Three sets of loops with increasing numbers of strands are explored – 8, 16, and 64 elements. The details of the resulting 18 loop configurations examined are shown in Table 2. For the configurations tested in this study, the mean apex temperatures of the loop (as calculated using Equation 8) is in the narrow range of between 2.71 – 2.96 MK, a consequence of the injected energy being approximately equal for each loop configuration (mean: \(1.10\times 10^{28}\) erg, \(\sigma \): \(0.04\times 10^{28}\) erg) despite the varying number of strands and nanoflares taking place.

Table 2 Loop-apex information for the 18 different setups examined.
Full size table

In Figure 3 the standard deviation of loop-apex temperature compared to their mean temperatures (Table 2) is shown for the six hour window displayed as a function of nanoflare repetition time. To understand these plots, consider a loop with a fixed number of strands; as the time between events increases (the repetition time increases), the magnitude of each individual nanoflare must also increase to maintain a fixed total energy deposition budget across the entire simulation. Consequently, more sparse events of larger size in the strands will create larger variations in the responding physical parameters of the loop. In contrast, many strands with a shorter delay between nanoflares have smaller energy-release magnitudes and subsequently alter the physical parameters of the loop significantly less – hence a significantly reduced standard deviation.

Figure 3
figure 3

The average standard deviation in temperature at the loop apex for all simulations plotted against nanoflare repetition rate.

Full size image

4 Synthesised AIA EUV Emission

Synthetic AIA emission is calculated for the 18 loop configurations shown in Table 2 by the following equation outlined by Sarkar and Walsh (2009):

$$ EM_{\mathrm{AIA}} = \sum _{i}^{N_{\mathrm{strands}}}\rho _{i}^{2} C \left (T_{i}\left (s\right )\right )\; \mathrm{d}s, $$
(11)

where \(i=1,2,3,\ldots\;N_{\mathrm{strands}}\) is the number of strands. \(EM_{\mathrm{AIA}}\) is the synthetic AIA emission, \(\rho _{i}\) is the strand density, and \(C\left (T_{i}\left (s\right )\right )\) is the AIA temperature response function, which is dependent on temperature [\(T_{i}\left (s \right )\)] and is measured for all positions \(s\) within the loops. The AIA temperature-response functions obtained from CHIANTI v9.0.1 are shown in Figure 4.

Figure 4
figure 4

AIA temperature-response functions used to generate synthetic channel temporal profiles for the 18 loop configurations examined.

Full size image

4.1 Examination of Time Lags Between EUV Line Pairs

Following the approach of Viall and Klimchuk (2011), the synthesised AIA emission is employed to examine any possible time lags between pairs of lines in the resulting emission-line time series. Here, eight gridpoints at the loop apex from each simulation are convolved to match the resolution of an AIA pixel, for which the AIA emission is then obtained from Equation 11 to produce the synthetic light curves; this is demonstrated in Figures 5 and 6. The cross-correlation of two light curves that are closest in temperature (e.g. AIA 211 and 193) are computed using the c_correlate function in IDL (Fuller, 1995). As with Viall and Klimchuk (2011), a positive (negative) time lag means a peak in emission occurs for the hotter AIA channel before (after) the cooler AIA channel emission peaks.

Figure 5
figure 5

Synthetic AIA channel time profiles over an example three-hour period for the lowest (left: 1500-second nanoflare repetition time on 8 strands) and highest (right: 250-second nanoflare repetition time on 64 strands) frequency heated loop scenarios in this study. The figure follows the same colour scheme as Figure 4 with AIA 94 (red), 335 (magenta), 211 (orange), 193 (blue), 171 (green), and 131 (black) shown in descending order for the two examples in normalised units with an appropriate offset so that the line profiles can be stacked.

Full size image
Figure 6
figure 6

Synthetic AIA channel time profiles over an example three-hour period for the lowest (left: 1500-second nanoflare repetition time on 8 strands) and highest (right: 250-second nanoflare repetition time on 64 strands) frequency heated loop scenarios in this study. The figure follows the same colour scheme as Figure 4 with AIA 94 (red), 335 (magenta), 211 (orange), 193 (blue), 171 (green), and 131 (black), which have been re-normalised with respect to the AIA 193 channel to show the relative emission of each channel.

Full size image

Mimicking Figure 7 of Viall and Klimchuk (2011), Figure 5 shows an example three-hour time window from the two most extreme cases of heating frequency simulated, as viewed across the AIA channels. The left (right) panel represents the case of a loop with the lowest (highest) heating frequency for 8 (64) strands that are subjected to nanoflare events with a minimum energy release of \(1.2\times 10^{25}\) erg (\(2.5\times 10^{23}\) erg) every \(\approx 1500\) seconds (\(\approx 250\) seconds). The emission in Figure 5 is normalised relative to their maximum intensity of each channel and then offset to produce a stacked plot. In contrast, Figure 6 displays the intensity of each channel normalised relative to the maximum intensity of the exceptionally bright AIA 193 channel from the three-hour time period under investigation.

In both of the extreme cases presented here, it is possible to follow the plasma evolution across the AIA channels throughout the three-hour period as it heats in response to the nanoflare events and subsequently cools. By examining the resulting time lags of the calculated emission between lines adjacent in peak emission temperature, any possible relationship of this to the input-heating frequency can then be explored.

Focusing on the LFH scenario presented in Figure 5 (left) and between \(t=0-1\) hour, a large rise and fall in emission can be seen in AIA 94, that is followed later by a rise and fall in AIA 335, and so on, which, following the nomenclature outlined by Viall and Klimchuk (2016), results in a positive time lag of \(\approx 480\) seconds between neighbouring AIA channel pairs based on the peak-temperature ordering shown in Figures 4 and 5. This time lag is due to the infrequent, large-magnitude energy release from each nanoflare (\(\gtrsim 1.2\times 10^{25}\) erg) generating large variability in the light curves. Conversely, the UHFH example (Figure 5 right) shows minimal variability in the light curves and subsequently no easily discernible time lag is present between the AIA channels due to the large number of smaller energy nanoflares taking place (\(\gtrsim 2.5\times 10^{23}\) erg).

In Figure 6, the light curves of Figure 5 are re-normalised to the maximum intensity of the AIA 193, which given the peak average loop temperature of \(\approx 2.7\) MK (Table 2) is exceptionally bright across all 18 loop configurations analysed in this study. In doing so, it can be seen that AIA 171, 193, and 211 have vastly superior counts compared to the hotter channels (335, 94, and 131). This is largely due to i) the loop parameters being tuned such that the resulting strand temperatures strongly correspond to the peak emission temperatures of the cooler AIA channels and ii) the hotter channels’ temperature response functions having lower relative peaks and broader distributions compared to the relatively cooler channels. Additionally, multiple peaks with contributions from cooler temperatures (\(\approx 1\) MK) complicate analysis of the hotter emission line intensities (O’Dwyer et al., 2010). Considering these factors, the following analysis will purely focus on the emission from the AIA 171, 193, and 211 channels.

Figure 7 provides a means for quantifying the light curves’ variability over the six-hour duration for the emission lines analysed. Here, the cross-correlation of the AIA channel pairs 211–193 and 193–171 are computed and summed together to generate a cumulative time lag to test whether any correlation between time lag and nanoflare repetition time can be identified.Footnote 1 As would be expected, a general trend can be seen whereby the time lag is shorter for loops whose strands are heated more frequently (smaller repetition time) compared to stands that are heated more infrequently (larger repetition time). For the loops comprised of the most frequently heated strands, the time lag between emission is \(\approx 250\) seconds, which is on the order of the input-heating repetition time. As the time between heating events increases, so too does the time lag, although for repetition times of \(\gtrsim 1000\) seconds the time lag significantly differs from the nanoflare repetition time. As the heating frequency further decreases towards our lowest frequency examples (1500-second repetition time), the deviation between repetition time and time lag increases with the time lag plateauing at around the 800 – 1000 second level in all of the multi-stranded scenarios in the case study.

Figure 7
figure 7

The cumulative time lag for emission detected in AIA 211 to cool and be detectable in the subsequent AIA channels 193 and 171 (as demonstrated for Figure 5) versus the nanoflare repetition time. All the strand-heating frequency scenarios tested are shown for loops consisting of 8 (blue), 16 (black), and 64 (red) strands. The green dashed line displays where the one-to-one correspondence between the heating frequency and time lag would occur.

Full size image

These time lags are a direct consequence of the competition between the fundamental cooling processes operating on the strand plasma (both radiative and thermal-conductive losses) and subsequent plasma draining versus the re-energisation arising from another nanoflare event. The fundamental cooling time of a strand (and hence the overall amalgamated loop) is directly proportional to the loop length (Barnes, Bradshaw, and Viall, 2019), which in this study remains fixed. Hence, similar time lags are observed for loops whose nanoflare repetition times are longer than this cooling time.

Therefore, given the above behaviour, it is likely very difficult to conclusively infer the heating rate of sub-element strands from any single specific time-lag observation unless the heating repetition time is much shorter than the overall characteristic cooling timescales. However, there is an alternative approach that requires utilising the possible detection of several similar observed time lags across the EUV channel pairs but over an extended time period, at least say several times the estimated cooling timescale.

Consider a simple, monolithic loop structure with a nanoflare repetition time greater than the loop cooling time; each nanoflare energy-release event in the loop would result in a rise and fall in the emission lines with the resulting calculated time lag of one event being close to the cooling time. However, if these time lags are observed repeatedly throughout the observational window, then the time between, say, the onset of the observed time lag in each line from one determined time-lag period to the next would provide an estimation of the nanoflare repetition time. However, the situation can become more involved when a multi-stranded approach is introduced, with the possibility of several nanoflares taking place nearly simultaneously and contributing to a confusing cumulative loop-emission profile.

With this in mind, consider the three-hour observational window shown in Figure 8 for the least-frequently heated loop in the case study (eight strands subjected to LFH). Here, the light curves of AIA channels 211 (orange), 193 (blue), and 171 (green) are plotted and there are five distinct periods (numbered 1 – 5) of rise and fall in emission that are associated with a positive time lag that are observable in all three channels. Measuring the time lag between emission-channel pairs from peak-to-peak yields a mean time lag of \(792 \pm 146\) seconds with a 95% confidence level. As is seen from Figure 7, this mean time lag corresponds well to the fundamental cooling timescale of the loop. Subsequently, calculating the repetition time between each of these time lags yields a value of \(\approx 1500\) seconds, which is equivalent to the mean period of nanoflare repetition. Thus, measuring the period between a series of successive time lags that are greater than the fundamental cooling timescale of a particular loop may then be used as a tool to indicate how frequently the loop is being heated.

Figure 8
figure 8

Synthetic AIA channel time profiles over an example three-hour period for the lowest (1500-second nanoflare repetition time on 8 strands) frequency heated loop scenario in this study for AIA 211 (orange), 193 (blue), and 171 (green). The emission has been time-averaged for 120 seconds and the five distinct rise and fall in emission from nanoflare heating are indexed in all three AIA channels.

Full size image

However, with that said it is worth noting that additional peaks can be seen, such as the peak between 4 and 5 in AIA 211 in Figure 8, which does not appear to correspond to any peaks in 193 and/or 171. This could be because this is an additional heating event that prolongs the duration for which the loop plasma remains at temperatures of ≈1.5 MK and, as such, a single peak is seen rather than two. Additionally, the time-series plots shown in Figure 8 only analyse a segment of the loop at the apex. Hence, the plasma at a particular point within the loop subjected to a heating event seen in one channel could be swept to another portion of the loop by mass flow, and subsequently the plasma may not appear in the neighbouring AIA channels. For this reason, the visual inspection of Figure 8 only considers the peaks in intensity that can be traced through the three AIA channels shown. This highlights the complexity of analysing the light curves of individual events within multi-stranded loops, and thus any analysis doing so must proceed with caution.

4.2 AIA Emission Channel Ratios

The mean AIA emission-channel ratios over the six-hour period analysed are shown in Figure 9 for 171/193 and 193/211 for all 18 loop configurations in the case study. Here, strong correlations are seen between the heating frequencies of each loop and their emission-channel ratios. As one would expect, strands that are subjected to fewer, but larger, magnitude nanoflares result in a loop that has greater 171/193 and 193/211 emission-channel ratios compared to a loop whose strands are subjected more frequently to lower magnitude nanoflares. This is because a loop that is heated more frequently has less opportunity to cool to lower temperatures before being re-energised compared to a less-frequently heated loop, as is evidenced by the standard deviation in apex temperature (Figure 3).

Figure 9
figure 9

Emission-line ratios for AIA 171/193 (left) and AIA 193/211 (right) for 8- (blue), 16- (black), and 64-stranded loops (red) shown as a function of nanoflare repetition time.

Full size image

4.2.1 Determining Heating Frequency and Strand Topology

In this subsection, a new technique is described that could prove a useful method for determining within a real, observed loop an estimate of the possible number of strands contained in the structure along with the subsequent heating frequency. This method combines AIA emission-channel ratio (Figure 9), standard deviation in loop-apex temperature over a period of time (Figure 3), and a numerical model such as the multi-stranded loop code.

Forward modelling a loop with AIA observations alone, it is only possible to measure the emission-channel ratios and thus it is only possible to constrain the nanoflare repetition time, as the number of strands appears to have little to no effect on the emission-channel ratios (Figure 9) for this study. In order to additionally constrain the number of strands that may be contained within a real AIA loop, co-aligned spectrometer data are needed such that a mean temperature and an associated standard deviation can be calculated and compared to the numerical models.

For example, consider a supposed observed loop that has a standard deviation in temperature of \(\approx 0.4\) MK over a six-hour time period with AIA 171/193 and 193/211 emission-channel ratios of 0.92 and 2.12, respectively. From Figure 3 there are three potential candidates that replicate the standard deviation in apex temperature, all of which have different heating frequencies and strand numbers (8 strands HFH, 16 strands ILFH, and 16 strands LFH). However, their emission-channel ratios for 171/193 and 193/211 are all distinctly different (see Table 3), and thus the model whose emission-channel ratios most closely match the observed loop could potentially reveal the heating frequency and number of strands contained within a real loop. In this example, the real loop would likely be comprised of 16 sub-element strands that are subjected to LFH.

Table 3 Emission-channel ratios for forward modelling.
Full size table

In order for these comparisons between observations and simulations to be made, there are several factors that must first be considered. Firstly, estimates of the loop length and radius must be undertaken. The former could be determined by performing non-linear force free-field extrapolation of EUV imager data yielding the magnetic-field line geometry; the length of the magnetic-field lines can then be determined and used as a proxy for loop length. The latter can be determined by measuring the cross-sectional widths of emission profiles, such as is demonstrated by Williams et al. (2020a,b).

Using co-aligned spectrometer data, estimates of the loop density and temperature can be made over a period of time. From this, the mean loop-apex temperature and the standard deviation could be calculated, and when combined with the loop geometry this can be used to generate several models that replicate these properties with different numbers of sub-element strands and heating frequencies. Then the AIA-channel ratios can be determined for the real and forward-modelled loops, and as previously outlined they can be compared to each other to estimate the number of strands and their most likely nanoflare repetition time(s). Finally, the inclusion of spectroscopic data would provide a further constraint on the forward-modelled loops by checking that the average (Doppler) velocities and densities of the observed loop are consistent with those generated by the model over the period analysed.

5 Summary and Concluding Remarks

The multi-stranded loop code employed by Sarkar and Walsh (2008, 2009) has been updated to model the numerically challenging transition region as a discontinuity (Johnston et al., 2017a,b). This allows the multi-stranded loop code to better resolve the interaction of the transition region with down-flowing plasma as a result of a nanoflare event. It is demonstrated that the new code leads to an increase in coronal velocity and density compared to previous versions (Sarkar and Walsh, 2008, 2009) – something that may often be underestimated in coronal-loop models in response to impulsive heating (Bradshaw and Cargill, 2013).

With an under-resolved transition region, artificially low coronal densities can often be obtained because the downward heat flux “jumps” across the unresolved region to the chromosphere, underestimating the upflows (see, e.g., Bradshaw and Cargill, 2013; Johnston and Bradshaw, 2019; Johnston et al., 2020). However, it is demonstrated that the updated code leads to an increased evaporative-upflow velocity when compared with the previous method. This results in a higher coronal density that shows good agreement with a fully resolved model, when subjected to impulsive heating (Bradshaw and Cargill, 2013). The parameter space is explored for loops comprised of a low (8), intermediate (16), and high (64) number of strands that are subjected to a range of nanoflare repetition times, which are consistent with those specified by Cargill, Warren, and Bradshaw (2015).

Synthetic AIA loop-top data are generated and the light-curve time series are explored using the time-lag analysis method (Viall and Klimchuk, 2011) for the cooler AIA channels (211, 193, and 171).Footnote 2 As noted by Viall and Klimchuk (2016), the time-lag technique will identify any cooling plasma present within a loop following an impulsive event by cross-correlating the rise and fall in emission of light-curve pairs. The light curves of an impulsively heated loop can result in a cyclical rise and fall of emission in light curves that might be detected and quantifiable using the time-lag method, although it has previously been demonstrated that the majority of time lags are consistent with cooling (Viall and Klimchuk, 2013). The results presented in Figure 7 reveal that the time lag typically increases as the nanoflare repetition time (and the magnitude of energy release per impulsive event) increases. However, loops that consist of fewer sub-element strands display variation as large as 100 – 200 seconds from one nanoflare repetition time to the next (compared to the 64-stranded loops). This variation makes it difficult to relate a particular time-lag value to a given repetition time. For example, the loop consisting of eight strands that is subjected to a nanoflare every 1500 seconds has a time lag of \(\approx 700\) seconds, which makes it not possible to differentiate from several other scenarios in the case study (16 strands subjected to ILFH; 16 and 32 strands subjected to IHFH). Similarly, the gradient of the time lag vs. repetition-time plots flatten for repetition times \(\gtrsim 1000\) seconds, and as such accurately determining the repetition time associated with a given time lag for even the 64-stranded loops is not viable.

With that said, an alternative approach to quantifying the nanoflare repetition times is outlined when the observed time lags are greater than the fundamental cooling timescales of the loop in question. It is demonstrated in this article that measuring the duration between successive time lags in this scenario provides an approximate estimation of the input-heating frequency. However, it is worth noting that for this method to provide a reliable estimate it is required that the input-heating frequency is sufficiently low such that the plasma in the strands making up the loop is able to cool through the channels being analysed.

Comparing the AIA channel ratios of 171/193 and 193/211 for the 18 loops (Figure 9) provides some promising initial results on obtaining the heating frequency of coronal loops. A strong correlation is found between the ratio values and nanoflare repetition times regardless of the multi-strandedness of the loop(s) in consideration. Comparing the AIA channel ratios with observations may prove to be a method by which the strand-heating frequency can be inferred by forward modelling. Furthermore, combining the model loops and AIA observations with spectrometer data will allow for mean and standard-deviation measurements of loop temperature, and, as shown in Table 3, this standard deviation, and the AIA channel ratios can constrain the number of strands and their heating frequency.

For example, if an observed loop has a standard deviation in temperature of \(\approx 0.4\) MK, then from Figure 3 there are three potential candidates that replicate this. However, their AIA channel ratios for 171/193 and 193/211 are all distinctly different (see Table 3), and thus the model whose AIA channel ratios most closely match the observed loop could potentially reveal the heating frequency and number of strands contained within a real loop. Without co-aligned spectrometer data to complement AIA, it would not be possible to accurately determine the mean loop temperature and the associated standard deviation, and thus an estimate of the number of sub-element strands could not occur. Additionally, given the large number of parameters at play (loop aspect ratio, loop length, temperature, geometry, Doppler/flow velocities, etc.) it is likely this analysis would need to be performed on a loop-by-loop basis. This will be explored further in a future study analysing Hinode/EIS loops (see Xie et al., 2017).

As is evidenced in this article and discussed in detail by O’Dwyer et al. (2010), there is currently poor spectroscopic coverage of plasma at temperatures \(\gtrsim 2\) MK. The upcoming NASA sounding rocket Marshall Grazing Incidence X-ray Spectrometer (MaGIXS: Kobayashi et al., 2017; Vigil et al., 2021) aims to fill this void by resolving soft X-ray spectra (Fe xvii – Fe xx) above a 4 MK active region. Future work utilising the data from MaGIXS, as well as the X-ray and EUV spectrometer data from ESA’s Solar Orbiter, is planned to further examine whether it is possible to determine the number of strands and/or their heating frequencies using the techniques presented in this article.