A Statistical Study of IRIS Observational Signatures of Nanoflares and Non-thermal Particles

Nanoflares, along with magnetohydrodynamic (MHD) Alfv$\acute{\rm e}$n waves, are thought to be an important mechanism to explain the heating of the solar corona, which, despite being of fundamental importance in astrophysics, is still not well understood (e.g., Klimchuk, 2006; Testa et al., 2015; Testa & Reale, 2022). Nanoflares can be produced by magnetic reconnection via braiding of coronal magnetic field lines caused by photospheric motions (Parker, 1988). The effect of these small energy releases (thought to be of the order of $10^{24}$-$10^{25}$ erg) is often difficult to observe in the corona because of the high-conductivity in the corona efficiently spreading out the energy, and the weak initial emission due to the low emission measure as well as non-equilibrium processes. The lower atmosphere of the coronal loops is more sensitive to the heating release in the initial phases and it provides useful diagnostics of coronal heating (e.g., Testa et al., 2013, 2014).

In the standard model of the solar flares, relativistic particles, which are accelerated by magnetic reconnection in the corona, penetrate along the magnetic fields and deposit energy in the chromosphere via collisions (e.g., Holman et al., 2011). Similarly, nanoflares are also expected to generate brightenings in the chromosphere and transition region, and, if indeed they are scaled-down version of larger flares, they might also show presence of accelerated particles. However, these non-thermal electron beams in nanoflares are likely characterized by shorter duration and smaller energies than in larger flares. Indeed, when non-thermal particles are observed in smaller heating events (nanoflare to microflares), their distributions are typically characterized by smaller low-energy cutoff and steeper slopes ($E_{\rm C}\sim 5-15$ keV, $\delta\gtrsim 7$; Hannah et al. 2008; Testa et al. 2014; Wright et al. 2017; Testa et al. 2020; Glesener et al. 2020; Cooper et al. 2021) than larger flares, and therefore also penetrating less deep in the low atmosphere.

Recently, small scale transient brightenings have been reported as the signature of nanoflares in high resolution and high cadence observations. Testa et al. (2013) reported that the footpoints of hot loops, known as moss, occasionally show high variability in Hi-C 193 Å data. These footpoint brightenings are characterized by timescales of the order of $\sim$15 s, which is much shorter than the known time scale of the typically observed moss variability (e.g., De Pontieu et al., 1999; Antiochos et al., 2003). Thus, smaller energies and shorter timescales, compared with larger flares, are required to explain these short term brightenings. The Interface Region Imaging Spectrograph (IRIS; De Pontieu et al. 2014) provided valuable clues on the heating mechanisms of transient hot loops (Testa et al., 2014, 2020). IRIS provides high temporal and spatial resolution UV spectra formed in the chromosphere and transition region. The comparison of the observed temporal evolution of spectral properties of IRIS lines with the prediction of numerical simulations of nanoflare heated loops indicated that the heating by non-thermal electrons is more plausible than conduction to explain some of the observations (Testa et al., 2014, 2020). Polito et al. (2018) and Testa et al. (2020) performed RADYN simulations to show the effect of several physical parameters of the nanoflares (e.g., duration, total energy, low-energy cutoff of non-thermal electron distribution), as well as the initial conditions of the loops, on the observable spectra. They showed that, when non-thermal electrons are present, the energy deposition layer in the lower atmosphere will be determined by the density of coronal loops and the hardness and total energy of the electron beam. Consequently, the properties of the heating also determine the sign of the Doppler velocities, and the intensities of the chromospheric and transition region lines. Bakke et al. (2022) recently used the same RADYN simulations to investigate additional chromospheric diagnostics of the heating properties based on ground-based spectral observations. These numerical models provide a useful framework to interpret the observed spectra.

The recent IRIS observational studies of these footpoint brightenings have analyzed a limited number (about a dozen; Testa et al. 2014, 2020) of events, manually selected, precluding any definitive conclusion on the general properties of nanoflares, and presence and properties of non-thermal particles, in active regions, outside large flares. In order to overcome these limitations we performed statistical studies of properties of nanoflare events, by using an automated selection procedure to identify a large sample of events, and therefore also reducing some selection biases possibly present in previous studies. For this purpose, we exploit the several EUV images observed by Atmospheric Imaging Assembly (AIA; Lemen et al. 2012) on board the Solar Dynamics Observatory (SDO; Pesnell et al. 2012) to detect the moss location, and investigate the IRIS spectra obtained at those locations for four different active regions. We apply a modified version of the automated algorithm of Graham et al. (2019) to automatically detect footpoint variability in co-aligned AIA-IRIS datasets of several active regions, and to extract the IRIS spectral line properties (intensity, Doppler shift, broadening) of Mg II and Si IV for all transient brightenings caught under the IRIS slit. Through the analysis of the statistical distribution of the physical quantities, and of the correlations between different parameters, we determined the general characteristics of nanoflares and their effects on the lower atmosphere. We also compare our finding with previous observational studies and numerical simulation results, and discuss the effect of nanoflares on the lower layer in terms on the plasma dynamics and heating mechanism.

We searched the IRIS database for observations suitable to detect transient loop footpoint brightenings and their chromospheric and transition region emission. We limited the search to sit-and-stare mode, including the Si IV 1402.77 Å and Mg II h&k spectral windows, with raster scan cadence of less than 15 s to ensure high cadence for the spectral observations. We searched for IRIS datasets including several observations of the same active region over several days, so we can also investigate the dependence of the observed heating properties on line-of-sight, and on the active region activity level. Using ssw_hcr_query.pro in IDL solarsoft, we searched the IRIS observing time for every active region from active region number 11800 to 13050 on the solar disk. It is possible that our approach means that data before the assignment of an active region number might not be included. In addition, only IRIS data coincident in time with relevant AIA data without gaps due to missing data or eclipses were selected. Here we present the analysis of IRIS observations for four suitable active regions.

As in Graham et al. (2019), the small scale transient brightening moss regions are identified in AIA data. We used the co-aligned AIA datacubes produced for each IRIS observation, and available on the IRIS search page¹¹1https://iris.lmsal.com/search/. After exposure time normalization, we applied the rapid varying moss detection algorithm (Graham et al., 2019) to the co-aligned AIA data. Here we briefly describe the algorithm. We identify the network region with AIA 1700 Å emission above a threshold ($>$80 DN). Then, we select datasets where hot loops are present because we want to investigate the moss variability related to coronal heating events producing hot loops. The Fe XVIII emission in the AIA 94 Å passband traces plasma hotter than $\sim$4 MK. In particular, to separate the Fe XVIII emission from the other cooler contributions to the 94 Å emission we employed the following simple relation using three different AIA channel data (Del Zanna, 2013):

We chose the pixels where Fe XVIII emission is greater than 5 DN, and the size of the coronal hot loop is larger than 100 pixels which is an empirical value for the minimum size of a coronal loop. To select bright moss emission, we identify the regions with AIA 193 Å intensity greater than 1250 DN. Even after all these criteria have been applied, it is possible that the selected areas might be associated with other phenomena such as filament eruptions or large flares. To eliminate these possibilities, we used the differential emission measure method (DEM, Cheung et al. 2015): we excluded areas where the emission measure in two temperature bands $\log T=5.6$ - $5.8$ and $\log T=6.7$ - $7.0$ is greater than $2\times 10^{26}$ cm^-5 and $2.5\times 10^{27}$ cm^-5, respectively, and their area is less than 15 pixels. We manually confirmed that this method effectively removes flares and filament eruption in most cases. Because we are interested in transient brightening events, we investigated the temporal variation in 171 Å and 193 Å and selected the local maximum in the light curve. From the first derivative of the light curve, we selected the times when its sign changed from positive to negative, known as zero-crossings. (e.g., Testa et al., 2013; Graham et al., 2019).

After identifying the position and time of the moss brightenings in the AIA data, we checked the corresponding IRIS slit position and observing time. If the IRIS dataset includes the 1400 Å slit-jaw images (SJIs), we repeated the alignment process of these with the AIA 1600 Å images to increase the accuracy of the IRIS-AIA alignment. We collected every IRIS raster pixel which is corresponding spatially and temporally, within 1 arcsecond and 6 seconds, to the selected AIA moss position and time. We set an additional criterion using the Si IV spectral line. We produce the Si IV light curves deriving the temporal evolution of the Si IV total intensity obtained by integrating over 1402.77$\pm 2$ Å. The light curves obtained for an interval of $\pm 150$ seconds from the moss detection instant in each pixel. Then, we investigated if an intensity peak exists within $\pm 30$ seconds from the moss detecting instant, and the intensity peak value is greater than 3$\sigma$ of light curve from the base intensity. We also checked that the intensity peak has a lifetime (FWHM) shorter than 60 seconds. If the Si IV light curve satisfies these three conditions the event is marked as a moss brightening and included in our sample. Then, using the chi-square values from the spectral fitting, as described later, we discard the bottom 1%, which effectively eliminates poorly fitted spectra. Finally, we conducted the investigation on the spectral properties at the peak time of the Si IV light curve in each pixel. For the Mg II spectral analysis, we used the peak time of Mg II k light curve in each pixel obtained by integrating within $\pm 0.65$ Å of its central wavelength. As a result, a total of 1082 pixels were obtained from 747 moss brightening events, and only the peak time properties were taken from each pixel. The information and number of selected pixels for each active region is summarized in Table 1.

We collect the spatial information for the IRIS moss brightenings. The heliocentric coordinate of the pixels provides information about the inclination of the local vertical with respect to the line of sight. Under the reasonable assumption that, on average, the magnetic field in the transition region footpoints of hot active region loops is typically vertical, this then provides a constraint on dependence with respect to the magnetic field direction. The spatially averaged Fe XVIII emission for whole data field of view (FOV) provided with IRIS data was collected as an auxiliary value for the environment of the event.

We extract the spectral parameters by fitting the IRIS spectra normalized by their exposure time. For the Si IV spectra, we use a single Gaussian function to fit the 1402.77$\pm 2$ Å wavelength range and obtain the three parameters – amplitude, centroid, and width – from which we derive the intensity amplitude, Doppler velocity, and non-thermal velocity, respectively. The non-thermal velocity was acquired by subtracting the thermal broadening and the instrumental broadening as follows

where $w_{nth}$ is the non-thermal velocity, $w_{obs}$ is the $1/e$ width of Si IV line which is corresponds to $\sqrt{2}$ times the Gaussian width, $w_{th}$ is the thermal velocity which corresponds to about 6.9 km s^-1 for the Si iv 1402 Å, and $w_{I}$ is the instrumental broadening which is order of 3.3 km s^-1 for the IRIS FUV spectral band (De Pontieu et al., 2014).

We also calculated the plasma electron density using the diagnostic based on the O IV 1399 Å and 1401 Å line ratio²²2https://iris.lmsal.com/itn38/diagnostics.html##density-diagnostics. We obtained the O IV total line intensities by integrating between 1399.77$\pm$0.25 Å, and 1401.16$\pm$0.25 Å respectively, and then calculated the line ratio. We exclude the samples of which total line intensity is lower than 3 sigma. The sigma is defined by following equation:

where $\sigma_{line}$ is sigma for the total line intensity, $N_{pix}$ is the number of spectral pixels for each O IV line, and $\sigma_{bg}$ is the standard deviation of background spectra in the wavelength range from 1404.25$\pm$0.25 Å, where no noticeable spectral line exists (Curdt et al., 2004). We collect the density information in 55% (593/1082) of the pixels. The others have the line ratio outside of the theoretically expected range between 0.17 and 0.42 for density diagnostics, or fail the 3 sigma detection criterion. These two groups have a large overlap.

For the Mg II h&k lines, we measured the velocity using two different methods. First, we calculate the centroid of the spectra using the center of gravity method in the range within $\pm 0.65$ Å of their central wavelengths, and convert these values to velocities,

where $v_{center}$ is the centroid velocity, $c$ is speed of light, $\lambda$ is wavelength, $\lambda_{0}$ is central wavelength of the lines, and $I_{\lambda}$ is observed spectrum. This centroid velocity roughly reflects the average motion within the chromosphere because the formation height of the Mg II h&k lines spans a relatively wide height range. Second, we measure the Mg II h3 and k3³³3See Figure 3.2 in https://iris.lmsal.com/itn39/Mg_diagnostics.html##the-mg-ii-h-k-lines Doppler velocities. As in Schmit et al. (2015), we adopted a combination of two symmetric linear functions, and a positive and a negative Gaussian functions to fit each line. This function is coded as iris_mgfit.pro in IDL Solarsoft. We determined k3 and h3 positions from the fitted spectral curve through iris_postfit_get.pro which is also included in the IDL Solarsoft, then these spectral positions were converted to the velocities. If the measured h3 and k3 velocities well represent the dynamics at $\tau=1$ layer like in numerical simulations of the quiet Sun (Leenaarts et al., 2013), they can provide valuable information about the upper chromospheric dynamics. However, the findings by Leenaarts et al. (2013) might not necessarily apply to the very dynamic active region plasma we are investigating here. Futhermore, in some cases the observed Mg II spectra are either single peaked, or with multiple ($>$2) peaks, so it is not always easy to determine the h3 and k3 positions. In fact, about 19% of Mg II k lines (206/1082) and Mg II h lines (203/1082) show single peak profiles, and about 16% cases (107/1082) show single peak profiles in both spectral lines. The single peaks are often accompanied by highly shifted tiny k3 or h3 self absorption features. To avoid misidentifications, we excluded those from the Mg II parameters analysis. On the contrary, the centroid method is free from this problem, so we mainly use centroid velocities to investigate chromospheric dynamics, and the h3 and k3 velocities are used as supplementary parameters. As we will discuss later, another advantage of the centroid analysis is that the comparison with the RADYN simulations (Polito et al., 2018; Testa et al., 2020) is more straighforward as the synthetic Mg II profiles are often quite complex, rendering determination of the k3 and h3 Doppler shifts difficult.

Another important feature is the Mg II triplet line at 2798.823 Å. Numerical simulations in Testa et al. (2020) showed that this line is expected to be observed as an emission line when non-thermal electron of sufficient energy ($\gtrsim 10$ keV) are present, as they will deposit energy and locally heat the lower chromospheric layer. In general, the Mg II triplet spectrum exhibits a complicate shape which is like a miniature of Mg II h&k line: it has dips at both wings and self absorption at the core similar to the Mg II h1&k1 or h3&k3, so it is not easy to determine whether it is an emission line or not. To quantify its behavior, we calculated the equivalent width of the Mg II triplet line by integration of spectrum over 2798.823$\pm$0.2 Å interval

where EW is equivalent width, and $I_{c}$ is the continuum determined by quadrature fitting within 2797.5 Å to 2802.5 Å range. We defined the Mg II triplet to be in emission if its equivalent width is a positive value.

For the comparison with numerical simulations, we adopted the models from Polito et al. (2018) and Testa et al. (2020), which are produced by using the RADYN code (Carlsson & Stein 1992, 1995, 1997a; Allred et al. 2015). That is a 1-dimensional hydrodynamic code including non-local thermodynamic equilibrium radiative transfer. We obtain the optically thin Si IV synthetic spectra as described in Testa et al. (2014), Polito et al. (2018) and Testa et al. (2020). The Mg II spectra are calculated by using the RH 1.5D radiative transfer code (Uitenbroek 2001; Pereira & Uitenbroek 2015) based on the atmospheric model from the RADYN simulations. The detailed information of the simulation models is summarized in Table 2. All models considered here assume the same initial atmosphere characterized by $\sim 1$ MK loop top temperature, coronal density of $\sim 5\times 10^{8}$ cm^-3, and loop length of 15 Mm (see Polito et al. 2018 and Testa et al. 2020 for more details). We analyzed the synthetic IRIS Mg II and Si IV spectral profiles at the peak time of their light curve using the same procedure applied to the observed spectra. This allowed us to derive spectral parameters of synthetic spectra that can be directly compared with ones derived from the observed spectra. We note that different models are characterized by different peak times for the transition region and chromospheric emission.

Figure 1 shows SDO/AIA observations of the four selected target active regions, and in particular, we show images in the photospheric intensity, line of sight magnetic field maps, 94 Å (in active region cores typically dominated by Fe XVIII emission; Testa & Reale 2012; Cheung et al. 2015), and 193 Å images (where active region moss is bright). Even though the active regions evolved as they crossed the solar disk during their life time, their characteristics are reflected in these images. For example, in the case of active region 12524, it shows a relatively simple configuration of the magnetic polarities, and it is well reflected in their Mount Wilson magnetic classification in Table 1. It shows dimmed 94 Å intensity compared to the other active regions. Throughout the observed period this active region was less active than the other active regions. Here we are interested in identifying and characterizing small heating events and exclude general flares, and we filter them out also using the DEM analysis, as described in the previous section. The number of observed flare events in each active region is summarized in Table 1. We found more small-scale transient brightening pixels in active regions where flaring activity is strong and the magnetic field configuration is complex. It is a well known fact that the flare activity of active regions is strongly correlated with their Mount Wilson magnetic classification (Toriumi & Wang, 2019, and reference therein). Although the probability for observing the brightening pixels is affected by the slit location and observing time, it appears that the small scale-transient brightenings have a relation to the magnetic field configuration as well.

The transition region response to coronal heating events provides very powerful diagnostics of the heating properties. Here we have carried out a statistical study of IRIS spectral observations of moss brightenings caused by small-scale heating events (nanoflares). We have further developed and algorithm to detect the moss brightenings automatically from the SDO/AIA and IRIS data (Graham et al., 2019). We have found more than 1000 moss brightenings, and analyzed the IRIS Si IV and Mg II emission, extracting their spectral parameters, and analyzing their distribution and the relations among different parameters.

The automatic detection algorithm allowed us to detect a broader range of events, compared with the manually selected sample we analyzed in previous work (Testa et al., 2014, 2020). Overall, we find that the parameters derived from our much large sample of events are mostly in agreement with the previous smaller studies. In particular, we note that the Si IV Doppler velocity has a similar, roughly symmetric, distribution between approximately $-40$ and $40$ km s^-1, and that the Mg II triplet emission (here identified by positive values for the equivalent width of the line) is found for a non-negligible fraction of the sample ($\gtrsim 30$ %). Some observed properties also show some differences. For instance the non-thermal broadening distribution of our sample covers a much broader range of values than previous work. We note that in Testa et al. (2020), rather than analyzing the absolute values of non-thermal broadening, we had derived the relative increases compared to the quiescent spectra outside of the brightenings. By taking that into account, the values found in the smaller sample of Testa et al. (2020) would correspond to a range of roughly 15 to 60 km s^-1, not too dissimilar to the bulk of our distribution, although their distribution peaked at the lower end, while we find a mean value of $\sim 32$ km s^-1. Also, we detected a large number of weaker events, as indicated by the histogram showing the distribution of Si IV total intensity (Figure 2e) and amplitude (Figure 2f). Testa et al. (2020, 2014) reported about a dozen of transient moss brightenings. For several of their events the maximum amplitude of the Si IV was $\gtrsim$ 40 DN pix${}_{x}^{-1}$ pix${}_{\lambda}^{-1}$ s^-1, slightly larger than our results with mean value of about 30 DN pix${}_{x}^{-1}$ pix${}_{\lambda}^{-1}$ s^-1. This preference for weaker events from the automatic detection algorithm is presumably due to the fact that smaller events occur more frequently. It may indicate that the brightenings we found are scaled-down version of large flares, and their occurrence rate perhaps follow a power law (e.g., Aschwanden et al., 2000). We briefly tested that the Si IV total intensity roughly follows power law distribution with power of -0.71, however, further studies are required to obtain the total energy of nanoflare from the Si IV spectral data.

We examined the dependence of the moss brightenings on the coronal environment. In Figure 2b, the flare productive active regions show relatively higher Fe XVIII intensity which indicates the existence of a larger amount of hot plasma in the corona. The Fe XVIII emission shows a moderate positive correlation with the total intensity of Si IV line, and with the equivalent width of the Mg II triplet line (Figure 7), which can be used as a rough indicator of the cutoff energy of the electron beam, according to the numerical simulations (Figure 2k, and Testa et al. 2020). In addition, we found that the moss brightenings occurred more frequently in the flare productive active region. Thus, we conjecture that in the vigorous active regions with hot plasma, higher energy electron beams can be generated more frequently.

As noted above, compared with our previous work based on small samples of these heating events, here we derived more parameters, and the size of our sample allowed for statistical analysis of the relations between the parameters. For the Si IV line (Figure 3), we found that the total intensity is well correlated with the Gaussian amplitude of the line, as expected. In our sample we found no significant correlation between non-thermal velocity and Si IV Gaussian amplitude, in contrast with results of previous work on an entire active region (De Pontieu et al., 2015). As we discussed above (sec. 3.2.1), this might be explained by the significant presence of multiple emission components or because the atmospheric structure disturbed by accelerated particles differs from that of the typical active regions. Another very interesting correlation we found is between the chromospheric (Mg II h&k) and the transition region Doppler velocities, which show significant positive correlation (Figure 4). This finding points to the fact that the heating event might impact the different layer of the atmosphere similarly. In particular, the Mg II h and k line have a very strong correlation of velocities, either measured by the line centroid or the position of the h3/k3 line. In both cases the velocity of the k line (formed at slightly larger heights than h) is slightly larger than for the h line, for both redshifted and blueshifted cases, suggesting a deceleration in both cases. However, it is interesting to note that the velocities from the Mg II centroid or the h3/k3 position only show a modest correlation (r=0.27) indicating that the Mg II h&k lines are mostly not symmetric, and, in turn, this suggests that these heating events are characterized by strong flows and gradients in the chromosphere.

To investigate the nature of the non-thermal velocity, we analyzed its relation with the distance from the solar disk center, i.e., with the inclination with respect to the line-of-sight. We found that the non-thermal velocity does not have a significant correlation with the distance from the solar disk center (Figure 5). Previous studies reported that center-to-limb variation of non-thermal velocity or line width is negligible (e.g., Chae et al., 1998; Ghosh et al., 2021). Our work confirms this lack of correlation is also valid in moss brightening regions. This could occur because the turbulent motions are isotropic, or there may be different processes that act along the field and perpendicular to the field and that are of equal magnitude (Tian et al., 2014; De Pontieu et al., 2014; De Pontieu et al., 2015).

Another finding about the Si IV non-thermal velocity is that the observed values are higher than numerical simulations, and than previous reports. The median value of our measurements is 32 km s^-1, and the distribution has a long tail to values that are higher than 60 km s^-1. In contrast, the simulation results show less than 25 km s^-1 (Figure 2g). This may be attributed to the fact that the simulations do not contain turbulent motions and assumed a single loop, while several thin adjacent magnetic strands, heated around the same time, are likely observed within a single IRIS pixel. Clearly such superposition within one pixel can lead to an increase in non-thermal broadening, given the different velocities along independently heated loops. In De Pontieu et al. (2015), they measured the non-thermal velocities of Si IV line in an active region, and reported that their distribution is close to a normal distribution with peak value of about 18 km s^-1. Our selected events are for moss regions, heated impulsively by conduction or accelerated electron beams, which might contribute to the high non-thermal broadening tail observed in active regions (De Pontieu et al., 2015; Testa et al., 2016). In some of the Si IV spectra in our sample, the very high non-thermal velocities are due to the presence of multiple components.

Thus, it is clear that the moss brightenings are usually characterized by violent dynamics. This is evidenced not only by the non-thermal velocity of the Si IV line, but also by the asymmetric profiles of Mg II spectral lines. An asymmetric profile implies that the h3 or k3 position is shifted with respect to the line centroid, so the upper chromospheric motion is not consistent with the average chromospheric motion. Moreover, the difference of h2v and h2r intensities indicates strong upward or downward motion (Leenaarts et al., 2013; Carlsson & Stein, 1997b).

By comparison with the numerical simulations, our results help to understand the nanoflare mechanism and provide constraint on the physical parameters. Our numerical simulations with various model parameters reproduce fairly well the Si IV and Mg II spectral properties when nanoflares occur. Most of the parameters have similar values to the observations. In agreement with our previous work (Testa et al., 2014; Polito et al., 2018; Testa et al., 2020), the observed Si IV blueshifted profiles, as well as the Mg II triplet emission are only reproduced by simulations including heating by beams of accelerated electrons. However, several parameters deviate from the observational data, pointing to interesting issues that need to be investigated in more details, and to the need for additional models that can more consistently explain the observations. Some examples include the relatively small non-thermal velocity, and the lack of significantly large negative value of the Mg II centroid velocities in the simulations (Figure 2g, 2i and 2j). As we speculated earlier in the paper, the former might be due to the fact that the single loop models might not be an adequate comparison for the observations in each IRIS pixel where many strands, possibly with significantly different initial conditions and heating properties, might overlap. The latter might be due to shortcomings of the model in reproducing realistic chromospheric conditions as discussed at length in previous literature (e.g., see discussions in Polito et al., 2018). Another puzzling discrepancy between simulations and observations is the relation between the Si IV Doppler velocity and the equivalent width of the Mg II triplet. The numerical simulations characterized by blueshifts also have stronger Mg II triplet emission, whereas the two parameters appear to have a very weak but positive correlation in the observational data (Figure 6c). It is possible that our definition of the equivalent width does not adequately capture Mg II emission in the observations, but this discrepancy will be thoroughly investigated in future work.

In summary this work has presented a thorough analysis of the observational signatures of the lower atmospheric (chromosphere and transition region) response to coronal small heating events. The comparison with RADYN simulations of nanoflare-heated loops has highlighted the overall agreement between the predictions and the observations for individual, single observables (e.g., intensity). However, we have also found some discrepancies, and there are no specific models that explain all observed parameters simultaneously. These limitations shed light on the shortcoming of the models and emphasize the need for more diverse and realistic numerical simulations for a more comprehensive reproduction of the observed signatures and interpretation of the data.