Statistical Study of the Kinetic Features of Supra-arcade Downflows Detected from Multiple Solar Flares

We have developed a tracking algorithm to determine the speeds of supra-arcade downflows (SADs) and set up a system to automatically track SADs and measure some interesting parameters. By conducting an analysis of six flares observed by the Atmospheric Imaging Assembly on the Solar Dynamics Observatory, we detect more smaller and slower SADs than prior work, due to the higher spatial resolution of our observational data. The inclusion of these events with smaller and slower SADs directly results in lower median velocities and widths than in prior work, but the fitted distributions and evolutions of the parameters still show good consistency with prior work. The observed distributions of the widths, speeds, and lifetimes of SADs are consistent with log-normal distributions, indicating that random and unstable processes are responsible for generating SADs during solar eruptions. Also, we find that the fastest SADs occur at approximately the middle of the height ranges. The number of SADs in each image versus time shows that there are “rest phases” of SADs, when few SADs are seen. These findings support the idea that SADs originate from a fluid instability. We compare our results with a numerical simulation that generates SADs using a mixture of the Rayleigh–Taylor instability and the Richtmyer–Meshkov instability, and find that the simulation generates quantities that are consistent with our observational results.


Introduction
The current sheet, which connects coronal mass ejections (CMEs) to post-eruption flare loops and is the place where magnetic reconnection takes place, plays a significant role during CME evolution Forbes et al. 2006;Lin et al. 2007;Ciaravella & Raymond 2008;Mei et al. 2012;Liu 2013;Reeves et al. 2019;Chen et al. 2020;Yu et al. 2020). The supra-arcade fan is the visible plasma sheet with hot (10 MK) plasma that appears above the observable flare arcades (Hanneman & Reeves 2014;Reeves et al. 2017;Cai et al. 2019). Turbulent motions in the fans have been reported by McKenzie (2013) and Freed & McKenzie (2018), and the plasma β (the ratio of gas pressure to magnetic pressure) in the fans can exceed unity in some cases (McKenzie 2013), suggesting that the turbulent flows may have an important impact on the magnetic fields in the fans.
In the fans, there are downward-moving dark voids, called supra-arcade downflows (SADs; Švestka et al. 1998;McKenzie & Hudson 1999;Innes et al. 2003). They were initially observed by the Yohkoh Soft X-ray Telescope (SXT; McKenzie 2000; Khan et al. 2007), then detected with a variety of instruments, such as the Transition Region and Coronal Explorer (extreme ultraviolet, or EUV; Asai et al. 2004 (Reeves & Golub 2011;Warren et al. 2011;Savage et al. 2012;Liu et al. 2013;Hanneman & Reeves 2014;Innes et al. 2014;Gou et al. 2015;Guidoni et al. 2015;Chen et al. 2017), which is sensitive to plasma at about 10 MK (Boerner et al. 2012). SADs are density-depleted structures, and most SADs are cooler than the surrounding fan plasma (Hanneman & Reeves 2014;Li et al. 2016;Reeves et al. 2017;Xue et al. 2020). SADs can contribute locally to the heating of the surrounding plasma by adiabatic heating (Reeves et al. 2017;Xue et al. 2020;Li et al. 2021) and viscous heating (Reeves et al. 2017), which are caused by the compression of the plasma and viscous motions, respectively. In addition, several studies have revealed that SADs correspond to hard X-ray bursts (Asai et al. 2004;Khan et al. 2007) and cause X-ray emission enhancement when they collide with the flare arcades (Samanta et al. 2021), suggesting the significant role of SADs in flare energy release.
In order to explore the process of SAD creation,  analyzed the frequency distributions of the sizes and fluxes of SADs by examining 16 flares observed by SXT. They found that the areas follow a log-normal distribution and the fluxes follow both a log-normal and an exponential distribution. They devised a simple SAD growth scenario with minimal assumptions, which explains the log-normal distribution.  examined a large sample of SADs and supraarcade downflowing loops (SADLs; retracting flux tubes viewed from a perpendicular angle) in eruptive solar flares from four instruments, and provided statistical samples of parameters such as the height, velocity, acceleration, area, elapsed time, magnetic flux, and shrinkage energy. They found that the speeds of the SADs are an order of magnitude slower than the assumed Alfvén speed, and concluded that SADs are post-reconnection loops. The higher-resolution data from AIA led to a refinement of this idea, where the SADs are the wakes behind the shrinking flare loops (Savage et al. 2012).
The origins of SADs are not fully understood, but studying the statistics of the kinetic features of SADs can help us to better understand SAD behaviors from a big-picture perspective, which will provide insights into possible SAD models. In this article, we investigate the statistics of the kinetic features of the SADs by examining six flares observed by SDO/AIA. The plasma velocities are calculated by using Gunnar Farnebäck's tracking algorithm (Farnebäck 2003), and the SADs are tracked and measured automatically by a system that we set up. Section 2 describes the AIA observations, and Section 3 describes the methods used in this work. Section 4 presents the results of the probability distributions, the evolution of the parameters, and the correlations between the parameters. Section 5 shows a comparison of our observational results with the results of a numerical simulation. The discussion and conclusions are presented in Section 6. Detailed descriptions of the six cases that we analyzed can be found in the Appendix.

Observations
The data for the six flares used in this study were obtained using the 131 Å channel of SDO/AIA. The AIA contains four telescopes, with four UV and seven EUV wavelengths, covering a full-disk field of view (FOV) up to ∼1.3 R e . These EUV channels (94, 131, 171, 193, 211, 304, and 335 Å), which are sensitive to a variety of coronal temperatures (from 10 5.5 to 10 7.5 K; O'Dwyer et al. 2010), take images at a temporal cadence of 12 s, with a spatial resolution of ∼0 6 per pixel. To process the AIA data, we first deconvolve the point-spread function from the original AIA Level 1.0 images by using the aia_deconvolve_richardsonlucy routine in the SolarSoftWare (SSW; Freeland & Handy 1998) package. Then we derotate, rescale, and align the AIA Level 1.0 data to create AIA Level 1.5 data, by using the aia_prep routine in the SSW package.
The six flare events that we analyze here occurred on the solar limb (from the SDO/AIA viewpoint), which is beneficial for us to examine the SADs. An overview of the SADs in the six cases observed using SDO/AIA 131 Å is presented in Figure 1. The information about the six flares is summarized in Table 1, in which we see the variety of eruption classes (GOES classifications), time spans, and total numbers of SADs in our samples. The N2Ts (the ratio of the number of SADs to the time span) of four cases (cases A, B, C, and E) with comparable flare ribbon lengths are very similar. Case D, with relatively shorter flare ribbons, has a smaller N2T, and case F, with relatively longer flare ribbons, has a larger N2T (the details are available in the Appendix).

Methods
In this section, we use case A as an example to illustrate the methods used in this work.

Velocity Measurements
The velocity results here are derived from an optical flow calculation. Optical flow is the pattern of object motion in an image caused by the relative movement between the object and the observer. 5 We utilize Gunnar Farnebäck's algorithm (Farnebäck 2003), which is a two-frame motion estimation algorithm based on polynomial expansion. The algorithm calculates the dense optical flow (the optical flow for all the points in the frames), which is beneficial for us in understanding the SAD properties and the surrounding fan environment. The open source code for the Farnebäck algorithm is available in the Open Source Computer Vision Library 6 (OpenCV), and we made some modifications to the original OpenCV code for better fitting in our work. These include: (1) changing the input interface from an image file to a fits file for better resolution and calculation accuracy; and (2) using the most suitable parameters, such as the size of the pixel neighborhood used to find the polynormal expansion and standard deviation of the Gaussian, for the better tracking of the SADs in the frames. The results of the modified algorithm show good performance in calculating the flow fields of the supra-arcade fans (see the movies in the .tar.gz package)-the Farnbäck algorithm reproduces the flow fields with lower noise than the previous method used in Reeves et al. (2017). An example of the calculation of the velocity results is shown inside the red box in Figure 2(a).

Method for Finding SADs
The method for finding SADs is based on the basic features of SADs, i.e., their moving downward and their darker appearance than the surrounding fan plasma in the AIA 131 Å channel. In addition to using the AIA Level 1.5 images, we implemented the Multi-Scale Gaussian Normalization method (MGN; Morgan & Druckmüller 2014) to enhance our image contrast. The MGN method, which is based on the localized normalizing of the data at many different spatial scales, is an applicable method for solar image processing. The MGN method that we used here reveals structures in the fans that are not easily noticed in the original AIA images, enhances the SAD edges, and therefore further improves the accuracy and efficiency of the automatic detection of SADs. As an example, the background of Figure 2(b) shows the MGN result from the AIA image inside the red box of Figure 2(a).
The thresholds of the velocity, intensity, and MGN value for finding the SADs in the six cases are slightly different, because the fan environments are different. For each case, we first test a couple of reasonable thresholds and review the results. We then choose the most appropriate thresholds for the best detection results. We do not extract the pixels of very faint SADs (usually they are SADs that just appear in the FOV), whose intensities are near the background noise level, because the velocity results of very faint SADs are not trustworthy. The process of extracting the pixels of the SADs and marking the SAD numbers involves the following four steps. First, we obtain the preliminary pixels that meet the conditions of the velocity exceeding 18 km s −1 and the intensity and MGN values being (typically) at least 30% darker than the surrounding fan. Next, we exclude sparse pixels by only reserving the pixels within the neighborhoods that have at least (typically) a 30% coverage rate of the preliminary pixels. We then label the separate blocks according to the connectivity among the selected pixels, and we only reserve the blocks that have at least 10 pixels as SAD candidates in an image. Finally, we identify the blocks in sequential images as the same SAD, according to the block's evolution. We mark the SAD numbers according to their appearance time for further statistical analysis. An example of the extracted pixels and marked SAD numbers can be seen in Figure 2(b). Movies of the velocities and extracted pixels of the SADs in the six events are available in the .tar.gz package.

SAD Measurements and Term Description
We set a "baseline," which is parallel to the looptop surface in the later phase of the flare, when the loop arcade is clearly visible, in order to calculate the height. We do not choose the solar surface as the baseline for calculating the height because the angle  between the looptops and the solar surface in most cases (see Figures 19-23) is not negligible. An example of a baseline is shown as the white dotted line in Figure 2(a). For an SAD, there is a specific point that has the closest distance between the baseline and the SAD-we call this point the SAD "head." We connect the heads of the corresponding SADs from the sequential frames to illustrate the SAD motion tracks as the solid colored lines shown in Figure 2(c). The tracks of the SADs show that the SADs move toward the baseline and further confirm that the way we chose the baseline was reasonable. In each flare, not all SADs move along the same direction (a good example can be found in Figure 24(a)), and even for an individual SAD, the direction of motion usually changes slightly over time, as indicated by the SAD motion tracks shown in in Figure 2(c). Generally, the SADs elongate along the direction in which they are moving during their evolution. These characterizations are implied by the velocity vectors in Figure 2(a) and the shapes of the SADs shown in Figure 2(b) (or in the movies in the . tar.gz package). Thus, in order to measure the parameters reasonably and effectively, instead of defining one fixed direction, we define the line starting from the head and going through the SAD's direction of motion (or the elongated direction, after the SAD has evolved for a while) as the "trunk" of the SAD, shown by the dashed lines presented in Figures 3(b) and 3(d). We do not take noisy and faint SAD tails into account, in order to avoid possible disturbances occurring at the top of the fan that could skew our subsequent width analyses. The trunk extends as far as there are distinct boundaries with the fan plasma, as seen in the SDO/AIA 131 Å or MGN images; an example is displayed in Figure 3. The height that we discuss herein is the distance between the baseline and the point in question, and the width is measured along the direction perpendicular to the trunk. Considering that the height is a relative value, depending on the baseline we choose, we do not directly compare the heights from different cases, to avoid introducing errors. As an example, the solid lines in Figures 3(b) and 3(d) display the lines that indicate the largest widths of the SADs. We see that as the SADs evolve, the trunk lines and width lines change accordingly, as shown at two times for SAD 4, SAD 8, SAD 13, and SAD 15 in Figures 3(b) and 3(d).

Results
In this section, we present the statistical results of the kinetic features of SADs. We present fitted distributions to the observed distributions of the widths, velocities, and lifetimes of the SADs in each image, to reveal information about SAD creation and evolution. We first examine the observed distributions for each flare, to avoid influences brought from other flares. We then combine the SADs from the six events to check whether the combination preserves or destroys the main characterizations of the observed distributions for an individual event. We also examine the heights, velocities, widths, and numbers of SADs as a function of time to see the overall evolution in the cases. In addition, we explore the possible correlations between widths and velocities, heights and velocities, and heights and widths, to attempt to provide some hints about SAD models. Figure 4 shows the fitted distributions for the SAD widths for the six cases. We measure the maximum widths for the SADs in each image, and the widths in the distributions are found by averaging these widths during the SAD lifetimes. The fitting of a log-normal distribution, shown as the curve in each panel, is performed by a maximum-likelihood estimator (MLE; Redner & Walker 1984;Venables & Ripley 2002;Myung 2003) to the unbinned data. The p-value of the Anderson-Darling test (AD; Anderson & Darling 1952, 1954Marsaglia & Marsaglia 2004), which rejects the hypothesis of a log-normal distribution if the value is less than 0.05, is denoted in each panel as well.

Probability Distributions
A characteristic of a log-normal distribution is the downturn of the frequency distribution for widths smaller than the median. The downturn in the distribution of Figure 4 appears clearly. Although there are different features among these cases (the medians and the width ranges at least are different, as shown in Figure 4), the lognormal is a good fit for the observed distribution of the widths in each case, which means that a log-normal distribution of the widths is an inherent attribute caused by SAD creation and evolution. A log-normal distribution is relevant when an observation is the product of fluctuations and random variables. In the solar physics context, log-normal distributions have been found in the widths of plasmoids during a CME observed by coronagraphs (Figures 2 and 4 in Guo et al. 2013; Figure 6 in Patel et al. 2020). In fluid mechanics, the late-time dissipation rate during the Richtmyer-Meshkov instability (RMI) also shows a log-normal distribution (Figure 14(b) in Weber et al. 2012), and log-normal distributions occur frequently across other branches of sciences, appearing in size distributions of aerosols, clouds, the pores in a cocoa press cake, the distribution of species abundance, and so on (Limpert et al. 2001). Figure 5 displays the fitted distributions for the maximum velocities of the SADs during their lifetimes for the six cases. Although the medians among the cases are not close, the lognormal distribution is a good fit to the observed velocity distributions, showing an inherent feature of the velocities during SAD creation. Similarly, a log-normal distribution of the plasmoid velocities during a CME is shown in Figure 8   The Alfvén speed in the lower corona is estimated to be ∼1000 km s −1 , and the outflows in the two-dimensional models of Sweet-Parker and Petschek reconnection have Alfvénic velocities. The velocities shown in Figure 5 are much slower than the Alfvén speed in the lower corona, which is consistent with the previous observational results in  and Warren et al. (2011), raising questions about whether SADs are really reconnection outflows. Figure 6 shows the maximum velocity versus maximum width, in which one dot represents one SAD, and the ρ and p-v in each panel denote Spearman's rank correlation coefficient (Best & Roberts 1975) and the corresponding p-value indicating sufficient  evidence of a significant correlation, if it is less than 0.05, respectively. We note that as the width increases, the velocity generally increases. Figure 6 shows that case D has a lower amount of small and slow SADs than the other cases, which is also demonstrated in the distributions of the widths (the width with the highest frequency in case D is higher than in the other cases) in Figure 4 and the velocities (the velocity median in case D is higher than in the other cases) in Figure 5.
We consider the velocities with displacements larger than one pixel (0.436 Mm) in two frames to be trustworthy, and since the time interval for calculating the velocities in two frames is 24 s, we estimate the velocity limit to be 18 km s −1 . From Figure 6, we know that as the width decreases, the velocity generally decreases, so there is another type of velocity limit as the SADs get smaller and approach the width limit. This limit on the velocity measurement occurs because the SADs below the limit would on average be too small to be measured, though not necessarily too slow to be measured. This type of velocity limit is not the same value for each case; we do not know the exact value of this type of velocity limit, but we can see the vertical dashed lines of the width limit in Figure 6 for reference (for simplicity, we call the former fixed velocity limit the "measured velocity limit," and the latter the "correlated velocity limit"). These two types of velocity limit combine together to decide the final minimum velocity that we can detect, and which type of velocity limit dominates depends on the specific cases. For example, the velocities of the SADs in case A are not close to the measured velocity limit, but as we can see, some dots in case A of Figure 6 are very close to the width limit indicated by the vertical dashed line. In this case, it is possible that there are slower but smaller SADs that we are unable to measure because of their size. On the other hand, the slowest SADs in case F are not very close to the width limit, but they are very close to the measured velocity limit, indicating that the limitation on the detection of the size of the SADs does not restrict the velocity measurement. The dashed lines in Figure 6, which indicate the width limit and the measured velocity limit, help us mark the "unreachable regions" in the cases, given the current observational conditions, as indicated by the line-filled regions in Figure 6. We note that the slowest SADs in five of the six cases (except case D) are already close to the boundaries of the unreachable regions, so the slowest SADs that we detect in these five cases are near the minimum velocity that we can detect, given the current observational conditions. Although the downturn of the frequency distribution for the smallest velocities in each case ( Figure 5) already manifests the feature of a log-normal distribution, future instruments with higher cadence and spatial resolution will definitely perfect the observed log-normal distributions of the velocities. Figure 7 shows the distributions of the lifetimes of the SADs. As we can see in Figure 7, relatively short-lifetime SADs dominate the whole process during each flare, and we fit both a log-normal distribution and exponential distribution to the observed lifetime distributions (further details can be found in the discussion about the fitted model selection around Table 4 in the Appendix). The lifetime median in case D is two to three times the medians in the other cases, because there are quite a number (8/39) of SADs whose lifetimes are longer than 30 minutes in case D, while we have a few SADs whose lifetimes are longer than 30 minutes in the other cases. Also, the bins above 30 minutes in the histogram are sparsely distributed, due to the small sample size of the SADs (39) in case D. The inconsistency of case D with the other cases might have been caused by the tilted fan from the observational viewpoint in case D (see Figure 22 and the discussion in the Appendix), Figure 6. The maximum width versus maximum velocity, where one dot represents one SAD, and the line-filled regions indicate unreachable regions in the observations. ρ denotes the Spearman rank correlation coefficient, and p-v denotes the corresponding p-value.
which resulted in the observation of some SADs that were obscured by other SADs, and therefore led to a smaller number of observable SADs ( Table 1) and the longer lifetimes of the SADs in case D than in the other cases.
The SADs that we can detect are limited by the temporal cadence of the instrument. We are not able to observe SADs with lifetimes of less than 12 s, and the size of the bins is limited by multiples of 12 s. In the future, if there are EUV imaging instruments with a temporal cadence faster than 12 s, then SADs with lifetimes shorter than 12 s might be detected, giving more evidence for SAD lifetimes being distributed as log-normal or exponential.
We further examine the distributions of the widths, velocities, and lifetimes of the SADs by combining the SADs from the six cases ( Figure 8). Figure 8(a) shows that the distribution of the widths has a median of 5.0 Mm, after combining the SADs from the six cases, and that the lognormal is still a good fit to the observed distribution of widths, which further confirms that a log-normal distribution of SAD widths is an inherent attribute of the SADs. McKenzie & Savage (2011) examined 120 SADs from different flares observed by SXT, and found that the SAD areas follow a lognormal distribution. We have measured width, not area, but the two results are nevertheless consistent in producing log-normal distributions of SAD sizes, even though the spatial resolution in our observations is four to eight times higher than in most of the observations in . This result suggests that the creation of the SADs is not a fractal process in which the measurement of the widths would yield a power-law distribution (Shibata & Tanuma 2001;Nishizuka et al. 2009) and there would be more detections of relatively small SADs with higher spatial resolution. Furthermore, there are similar distributions of SAD widths whether we combine the cases or not, indicating that the major processes governing the widths of the SADs from different flares are the same, and that it is reasonable to study the general features of SADs by combining SADs from different flares, as McKenzie & Savage (2011) and A log-normal distribution of the SAD velocities still exists, with the median of 66.8 km s −1 (Figure 8(b)), after we combine the SADs from the six cases. The result of their Figure 5(b)) also showed a log-normal distribution of the SAD velocities measured from dozens of flares, and we note that the velocity median in  is 146.13 km s −1 , which is higher than our velocity median. One reason for the difference in the velocity medians is that  have a higher fraction of energetic flares in their sample (see Table 1 in . As a result, they had quite a few SAD contributions whose velocities exceeded 350 km s −1 , with a highest velocity of about 1000 km s −1 , while we do not observe SADs with velocities exceeding 350 km s −1 in our six cases at all. Also, Figure 6 shows that the velocities of the SADs are correlated with the widths of the SADs, and that the SADs whose widths are generally smaller than 3.6 Mm have velocities that are generally slower than 120 km s −1 . The spatial resolution of the work in  is about 3.6 Mm, and, as we can see from Figure 6, we have a significant amount of SADs (27%) whose widths are smaller than 3.6 Mm with velocities slower than 120 km s −1 . Therefore, the contribution of the smaller and slower SADs due to the higher spatial resolution of SDO/AIA is another reason that the velocity median in Figure 8(b) is lower than the velocity median of . This explanation can be further confirmed by noting that the velocity median in case D is much higher than in the other cases ( Figure 5), and is closer to the result of , and that there is no SAD with a maximum width smaller than 3.6 Mm in case D ( Figure 6).
It is interesting to note that the lifetimes of the SADs follow a log-normal distribution with a median of 4.0 minutes, after we combine the SADs from different cases, as shown in Figure 8(c), while it is not possible to distinguish between lognormal or exponential distributions for SAD lifetimes when examining an individual event (further details can be found in the discussion about the fitted model selection around Tables 4 and 5 in the Appendix). Since SADs are found in the context of a fan, we note that the fan environment varies from region to region. For example, as shown in Figure 2, the detected height of the upper fan boundary is lower in the part of the fan near arrow B than it is in the part of the fan near arrow A. Given that fans are usually brighter at lower heights, the SADs initially appearing around region A are at a great height with a low surrounding intensity, and those appearing around region B are at a lower height with a greater surrounding intensity. It is unclear why the fan around location B does not extend higher, as it does near location A, but it must have something to do with the local environment above the fan. The various local environments of each zone of the flare fan lead to the general property that SADs that are initially detected at lower heights appear surrounded by brighter plasma, as displayed in Figure 9.

Evolution of the Parameters
Compared to the overall trend shown in Figure 9, the increase in the initial detected heights as a function of time for a certain color (with a comparable intensity in the surrounding plasma) is more pronounced. If we assume that SADs with a comparable initial detected intensity in the surrounding fan are formed under the same conditions, then the increase over time in the initial detected heights of the SADs with same color (with a comparable intensity) of surrounding fan suggests that the regions with the same local conditions or the sites forming SADs above the AIA fan also increase over time.
The initial detected heights of the SADs as a function of time in case E, however, do not show an increase, which is mainly attributed to the fact that the heights of some of the fan regions in case E decrease monotonically over time. It is possible that some energy transfer processes in the fan in case E lead to a decrease of the fan height during the period of the SAD's appearance, which is different from the evolution of the fan height in other cases. It is also possible that an increase in the initial detected heights of the SADs over time still exists in case E, but that there are undetected SADs that disappear before they arrive at the observed fan during the later phase, because part of the fan is out of the temperature response range of the 131 Å imager, especially since there are almost no SADs in the faint part of the fan after 80 minutes, while the initial detected heights of bright SADs (red) show an increasing trend. However, we do not see any "missing parts" of the fan in SDO/AIA 94 Å, which mainly has a lower response temperature than SDO/AIA 131 Å. It is very hard to tell if any "missing parts" of the fan appear in SDO/AIA 193 Å, which mainly has a higher response temperature than SDO/AIA 131 Å in flares, since SDO/AIA 193 Å is also very sensitive to the background corona at 1.5 MK, and, in this event, the background remains very bright, since the time before the beginning of the flare. Therefore, it is not clear here whether the decrease in the fan height over time observed by SDO/AIA 131 Å in case E is due to density depletion or heating. Future studies on temperature-related features in fans will help us to understand the evolution of detected fan heights, and further reveal the relation between the initial heights of the SADs and the fan heights.
The  maximum velocities and the time of appearance of the SADs in our results (Figure 10). In general, our results show that the fastest SADs do not appear at the earliest times, but that it instead takes several tens of minutes (28-50 minutes for cases A-E) after the SADs are initially observed for the SADs to reach their highest velocities. These observations provide insights into the turbulent features of fans compared with the hypothesis that the SADs are outflows of magnetic reconnection, based on the classic eruption model. As the flare progresses, the region that outflows the encounter at the top of the arcade rises. Given that the reconnection outflow speed is generally inversely proportional to the height within a current sheet (Lin 2004), if the SADs were the outflows of magnetic reconnection, then the maximum velocities of the SADs would monotonically decrease over time. We also note, in case F, that a group of the fastest SADs appear at the earliest time. Nevertheless, the maximum velocities of the SADs in case F do not show a pronounced decrease over time. The difference between our results and those of  is mainly because we have detected more slower SADs, and we notice that a few very fast SADs appear in the middle evolution time in the results (Figure 8(d)) of  as well. Figure 11 shows the maximum width evolution of the SADs, and it does not support a correlation between the maximum width and the corresponding time, which is consistent with the results of their Figure 9(b)). But we note in our results that the widest SADs usually appear in the early evolution phase of the fans, if we exclude case D. The inconsistency of case D with the other cases might be caused by the tilted fan from the observational viewpoint in case D (see Figure 22 and the discussion in the Appendix), which resulted in the observations of some SADs being obscured by other structures.
It is interesting to check the number of SADs versus time, to see the overall evolution of the system. As we can see from Figure 12, the evolution of the number of SADs is made up of several active phases, interspersed with rest phases (when few SADs are seen) in between. In an active phase, the number of SADs first increases, then it arrives at the maximum number in the active phase, then the number decreases, and finally the system arrives at a rest phase, which is often followed by another active phase. The evolution of the number of SADs reveals information about the productivity of SAD creation. The most productive active phase is not necessarily the first one, as we can see in Figure 12. The most productive active phase occurs during the first active phase in cases B, E, and F, but it occurs in the later active phase in the other cases. The existence of several active phases in each case indicates that the creation process of the SADs is intermittent. Similar features corresponding to rest phases have similarly been reported in the evolution of the number of plasmoids in a current sheet (see Figure 2 of Lynch et al. 2016).
In an extreme case, if a rest phase is durable enough, then no new SADs appear around the time period, and the previous SADs in the fan have already disappeared, meaning that we can have some frames with no SADs. The frames with no SADs are obvious at some times in case C, such as at about 8, 28, 74, and 93 minutes. Figure 13 shows the maximum velocity versus corresponding height in each case, and it does not support a correlation between those parameters. We normalize the distance between the minimum height and the maximum height as 1, and denote the normalized height with the maximum velocity in each case (in case F, we also denote the height with the second maximum velocity, since the dot with the maximum velocity may be an Figure 9. The initial detected height evolution of the SADs, where we set the appearance time of the first SAD as time 0, and one dot represents one SAD, with the reddish color representing a greater initial detected intensity of the surrounding fan. ρ denotes the Spearman rank correlation coefficient, and p-v denotes the corresponding p-value.

Correlations between the Parameters
outlier, so the height of the second maximum velocity is more reliable). The denoted heights are close among these cases, in the sense that the fastest SADs usually appear in the middle of the range of heights. A similar hill-shaped distribution of the velocity versus height is shown in Figure 8(c) of  as well. These hill-shaped distributions of the velocity versus height, which show the various velocities for a certain height, do not support the model of SADs being the outflows of magnetic reconnection, in which the velocities would decrease as a function of height. The normalized height of the fastest SAD in case F is smaller than in the other cases, partly because the fastest SADs in cases A-E appear at 28 ∼ 50 minutes after the initial observed SADs, the fastest SAD in case F is at the time of the initial observed SAD (Figure 10), and the initial detected heights of the SADs generally increase over time (even in case E, the decay of the initial detected height over time does not occur before the appearance of the fastest SAD), as shown in Figure 9. Figure 14 shows the maximum SAD width plotted against the height. The correlations between the height and the width vary considerably among the cases. There are correlations in cases C, D, and E, where we can see that the widths generally increase as the heights increase, and there are no clear correlations in the other cases. The plasma properties of the fan may contribute to how the widths of the SADs change as a function of height, which is a question that we will explore in a later paper.

Comparison with a Numerical Simulation
The numerical simulation we compare with here is from Shen et al. (2022), in which they perform a threedimensional magnetohydrodynamic simulation with a realistic flare-reconnection geometry to explore the origin of SADs. They postulate that SADs are self-organized structures formed in the turbulent interface region below the lower tip of the current sheet due to the Rayleigh-Taylor instability (RTI) and the RMI. We use the system that we set up in the observations to detect the SADs and measure the parameters in the synthetic SDO/AIA 131 Å images in the numerical simulation (but we can obtain the velocities from the simulation results directly).
Similar to the observational results, we plot the distribution of the SAD widths by averaging the maximum widths during the SAD lifetimes in Figure 15, where L 0 is the changeable characteristic length in the system (in Shen et al. 2022, they use L 0 = 150 Mm). It turns out that the width distribution in the numerical simulation also follows a log-normal distribution, with a median width of 0.0063 L 0 , which shows good consistency with the observational results shown in Figures 4 and 8(a). The fact that the width reasonably follows a log-normal distribution in the numerical simulation is attributed to the SADs being formed in the turbulent interface region, due to a mixture of the RTI and the RMI, a process that involves random fluctuations. Figure 16 shows the evolution of the initial height of the SADs in the numerical simulation, where one dot represents one SAD, with the reddish color representing a greater initial detected intensity and t 0 being the changeable characteristic time of the system (in Shen et al. 2022, they use t 0 = 109.8 s). We can see from Figure 16 that the initial heights of the SADs generally increase over time, and that the increase for a certain color (with comparable intensity) is somehow more pronounced, which is consistent with the observational results shown in Figures 9 and  8(e) of . The initial heights of the SADs as a function of time show an increase in the simulation, because the height of the turbulent interface where the SADs can be found generally rises over time.    fastest SADs, normalizing the distance between the minimum height and the maximum height as 1. Figure 17 does not suggest a correlation between the height and the velocity of the SADs, since, for a particular height, there are SADs with various velocities. The hill-shaped distribution, with the fastest SADs being situated at approximately the middle of the height Figure 13. The maximum velocity versus corresponding height, where one dot represents one SAD. The dashed lines and texts mark the relative heights of the fastest and second-fastest (case F) SADs, with the labels indicating the normalized distance between the minimum height and the maximum height. Figure 14. The maximum width versus corresponding height, in which one dot represents one SAD. ρ denotes the Spearman rank correlation coefficient, and p-v denotes the corresponding p-value. range, resembles the distributions obtained from the observations (Figures 13 and 8(c) in . The fastest SADs that appear at approximately the middle height in the simulation are mainly related to advantageous conditions for the instability growing at that height (interested readers are referred to Shen et al. 2022). In addition, the characterization that the speeds of the SADs are much slower than the assumed Alfvén speed in the observations ( Figure 5, Warren et al. 2011) is well demonstrated in Figure 17. The reconnection outflows and SADs in Shen et al. (2022) have distinctively different velocities: the velocities of reconnection outflows are near the local Alfvén speeds, and the velocities of the SADs are 10%-40% of the velocities of the reconnection outflows. Figure 18 shows the evolution of the number of SADs in each image in the numerical simulation, where some minima that may correspond to rest phases in the observations are pointed out by blue arrows. The SADs in the simulation are self-organized structures formed in the turbulent interface region due to the RTI and RMI, and the turbulent interface groups appear intermittently, so the SADs in the simulation appear in different groups, resulting in the rest phases shown in Figure 18.

Discussion and Conclusions
In this paper, we present a statistical study of the kinetic features of SADs detected from six flares. We determine the speeds of the SADs by using a modified Gunnar Farnebäck algorithm, and we identify the SADs and measure the parameters (height, width) automatically using the system that we set up. The distributions of the widths, speeds, lifetimes, and numbers of SADs in each image yield good fits of lognormal distributions, both when we examine the SADs in each event separately and when we combine the SADs from different events together, indicating that the log-normal distributions in these parameters are consequences of the inherent processes of the creation and evolution of SADs. A log-normal distribution is relevant when an observable is a product of fluctuating and random variables, indicating that random and unstable processes are involved in the creation and evolution of SADs. Unlike a normal distribution, which has a Figure 15. The distribution of the SAD widths by averaging the maximum widths during the SAD lifetimes in the numerical simulation. The curved line shows the fit result of a log-normal function performed by an MLE to the unbinned data, M denotes the median value of the widths, and the p-value is obtained from the AD test. L 0 is the changeable characteristic length in the system. Figure 16. The initial height evolution of the SADs in the numerical simulation, where one dot represents one SAD, with the reddish color representing a greater initial detected intensity, and the curve being the maximum height of the fan as a function of time. ρ denotes the Spearman rank correlation coefficient, and p-v denotes the corresponding p-value. t 0 is the changeable characteristic time of the system, and L 0 is the changeable characteristic length in the system. Figure 17. The maximum velocity versus corresponding height in the numerical simulation, where one dot represents one SAD. The dashed lines and text labels mark the relative heights of the fastest and second-fastest SADs, normalizing the distance between the minimum height and the maximum height as 1. L 0 is the changeable characteristic length in the system, and v A is the changeable characteristic Alfvén speed of the system. symmetric distribution, the log-normal distributions of these parameters show skewed distributions, suggesting that the systems prefer relatively narrow, slow SADs, with short lifetimes, and there are a relatively low number of SADs in each frame (Figure 8).
The width distributions in our work show good consistency with the area distribution of Figure 1(a) in , even though our data have generally higher spatial resolution, indicating that the creation process of the SADs is not a fractal process, where the distribution would be a powerlaw distribution and we would detect more instances of relatively small SADs with the higher spatial resolution.
The log-normal distributions of the velocities are shown in both Figure 5(b) of  and Figure 5 and Figure 8(b) of the present work, but we note that the velocity medians ( Figure 5 and Figure 8(b)) in our results are slower than those in . By examining the plots of the velocity versus width (Figure 6), we find that the velocities of the SADs generally increase as the widths increase. Additionally, there are a remarkable number of smaller and slower SADs (>20%) in our observations, exceeding the spatial resolution limit of most of the observations in , so we have more instances of smaller and slower SADs, giving rise to lower medians of velocities than in prior work.
Generally, the initial detected heights of the SADs increase over time, and the widths of the SADs do not monotonically change over time, which is consistent with the prior work of . However, our results show weaker correlations between the velocities and time than the result shown in Figure 8(d) of . On the other hand, the fastest SAD in  occurred in the middle of the evolution phase. The fastest SADs appearing in the middle of the evolution phase might be related to advantageous conditions for the instability growing around a certain time period, but this idea still needs further investigation in the future.
By checking the evolution of the number of SADs in each image, we find that the evolution of the number of SADs is composed of several active phases with rest phases, with a few SADs appearing in between, indicating that the creation process of SADs is intermittent.
The plots of the velocity versus height do not support strong correlations, but they instead show hill-shaped distributions, with the fastest SADs appearing at approximately the middle of the height range. Moreover, the correlation of the width versus height is only shown in some of the cases that we examine. The widths of the SADs may correlate with the local pressure in the fan, rather than their heights, which we will investigate in a future work.
Finally, we implement the system that we set up for the observations to measure the parameters of the SADs in a numerical simulation, and compare the results with the observations. As in the observations, we find that the width of the SADs in the simulation follows a log-normal distribution, because the SADs in the numerical simulation are formed in the turbulent interface region, due to a mixture of RT and RM instabilities, a process that involves random fluctuations. The initial heights of the SADs in the simulation generally increase over time, and the detection of the SAD heights is constrained by the maximum height of the fan. The velocity of the SADs versus their height shows a hill-shaped distribution, with the fastest SAD appearing at approximately the middle of the height range in the simulation, which might be related to advantageous conditions for the instability growing at the middle height (interested readers are referred to Shen et al. 2022). These results generally reproduce the features shown in the observations. Moreover, the evolution of the number of SADs shows that there are rest phases, as in the observations, and these rest phases can be well explained by the mechanism of SAD creation in the simulation, as SADs are self-organized structures formed in the turbulent interface region, and the turbulent interface group appears intermittently. We conclude here that the model proposed by the numerical simulation is consistent with our results.
We note here that the flares examined in this work show several retracting loops at the heads of the SADs, features that were elaborately illustrated in Savage et al. (2012), who suggested that the SADs are the wakes of the retracting loops. On the other hand, the morphology and spatial positions of the retracting loops demonstrated in Savage et al. (2012)-or in this work-are also quite compatible with the three-dimensional model of Shen et al. (2022;see Figures 1(A) and 4(B) and the related discussion therein.) Reeves et al. (2017) found that in the region of the fan without SADs observed by SDO/AIA, the plasma cools at a rate that is slower than the estimated conductive cooling, which may be the consequence of unobserved flows. If the SADs are self-organized structures formed in the turbulent interface region, as in the scenario provided by the numerical simulation of Shen et al. (2022), that means that we could have undetected flows in the fans with our current instruments. Though turbulent motions have been reported by McKenzie (2013) and Freed & McKenzie (2018), imaging instrumentation in hot bandpasses with higher spatial resolution will in the future be useful for further examining turbulent motions, testing models, and determining the role that SADs (or flows) play in energy transfer in fans. We find that the statistical kinetic features that we present in this paper can provide some clues from a different perspective for future SAD models to consider.
The authors thank the anonymous referee for comments that improved the paper. We are grateful to Dr. Vinay Kashyap for valuable comments and suggestions. The work of X.X. is supported by the scholarship of the China Scholarship Council, the Chinese Academy of Sciences (CAS) grants XDA17040507 and QYZDJ-SSWSLH012, NSFC grant 11933009, and the grants associated with the project of the Group for Innovation of Yunnan Province 2018HC023 and the Applied Basic Research of Yunnan Province 2019FB005. K.R. is supported by NSF grant AGS-1923365. C.S. is supported by NSF grant AST-2108438. J.I. is supported by the NSF-REU Solar Physics Program at SAO grant AGS-1850750 and the Chandra grants AR9-20004X and AR6-17003X. The AIA data are provided courtesy of NASA/SDO and the AIA science team. The SolarSoftWare (SSW) system is built from Yohkoh, SOHO, SDAC, and Astronomy libraries, and draws upon contributions from many members of those projects.

Appendix A Fitted Model Selection
The fitted model selection between log-normal and exponential models in the distributions of the SAD widths, velocities, and lifetimes is performed by the Bayesian Information Criterion (BIC; Schwarz 1978). If the ΔBIC between two models is greater than 10, then the model with lower BIC gives a better fit; if not, they are indistinguishable (Kass & Raftery 1995).
From Table 2, we note that the BIC of a log-normal model is lower than the BIC of an exponential model, and ΔBIC > 10 in all six cases of the SAD widths that we examined. Therefore, we confirm that the SAD widths follow a log-normal distribution rather than an exponential distribution, in each case, as shown in Figure 4.
From Table 3, it is easy to notice that the BIC of a lognormal model is lower than the BIC of an exponential model, and ΔBIC > 10 in all six cases of the SAD velocities that we examined. Therefore, we confirm that the SAD velocities follow a log-normal distribution rather than an exponential distribution, in each case, as shown in Figure 5.
From Table 4, we note that ΔBIC > 10 only occurs in one of the six cases, meaning that there is not enough evidence to select either the log-normal or exponential distribution for the SAD lifetimes in most cases. Thus, we plot both fitted lognormal and exponential distributions in Figure 7.
From Table 5, we note that the BIC of a log-normal model is lower than the BIC of an exponential model, and ΔBIC > 10 in the SAD widths, velocities, and lifetimes after we combine the SADs from the six cases. Therefore, we confirm that the SAD widths, velocities, and lifetimes follow a log-normal distribution rather than an exponential distribution after we combine the SADs from the six cases, as shown in Figure 8.

B.1. Case A
On 2011 October 22, a flare of magnitude M1.3 erupted on the northwest of the Sun (case A). It started at about 10:00 UT, and the X-ray intensity achieved its peak at about 11:10 UT (Figure 19(c)). The decay phase lasted for several hours, during which dozens of SADs appeared. The orientation of the arcade loops can be seen in the observation from STEREO-Ahead (STEREO-A; Figure 19(b)), which shows two flare ribbons almost parallel to each other. This observation, combined with the angular separation between the Earth (SDO is located in a nearly geosynchronous orbit, so the Earth's viewpoint can represent SDO's viewpoint) and STEREO-A (Figure 19(b)) during the time period of the flare eruption, indicates that there are relatively few projection effects from the SDO/AIA viewpoint. Clear SADs are seen in the SDO/AIA 131 Å image, as shown in Figure 19(a). The time period and the area of the fan in which we examine the SADs are indicated by the shaded region in Figure 19(c) and the red box in Figure 19(a), respectively.   . It started at about 00:15 UT, and the X-ray intensity achieved its peak at about 01:20 UT (Figure 20(b)). The decay phase lasted for several hours, during which dozens of SADs appeared, as indicated in the snapshot observed from SDO/AIA 131 Å shown in Figure 20(a).
Compared with case A, the boundaries of the SADs in case B are less clear and the fan is more turbulent overall, as Figure 19(a) and Figure 20(a) show. The time period and the area of the fan in which we examine the SADs are indicated by the shaded region in Figure 20(b) and the red box in Figure 20(a), respectively. Unfortunately, STEREO did not have observational data for this flare.

B.3. Case C
On 2012 July 17, a flare of magnitude M1.7 erupted on the southwest of the Sun (case C). It started at about 12:00 UT, and the X-ray intensity achieved its peak at about 17:00 UT (Figure 21(c)). Unlike case A and case B, the increase of the soft X-rays in case C was long in duration and moderate, and dozens of SADs mainly appeared during the rise phase of the flare, as Figure 21(a) and Figure 21(c) indicate. The orientation of the arcade loops observed from STEREO-A 195 Å and the angular separation between the Earth and STEREO-A are shown in Figure 21(b). The time period and the area of the fan in which we examine the SADs are indicated by the shaded region in Figure 21(c) and the red box in Figure 21(a), respectively.

B.4. Case D
On 2012 January 16, a flare of magnitude C6.5 erupted on the northeast of the Sun (case D). It started at about 02:30 UT, and the X-ray intensity achieved its peak at about 04:00 UT (Figure 22(c)). The increase and decay of the soft X-rays were long in duration and moderate, and dozens of SADs appeared before and after the X-ray flux peak. The arcade loops observed from STEREO-Behind (STEREO-B) 195 Å (Figure 22(b)) show that the flare ribbons are relatively short and not that parallel to each other, compared to other cases in this work.
Combined with the angular separation between the Earth and STEREO-B shown in Figure 22(b), we expect there to be projection effects, due to the tilted fan in the SDO/AIA observations. The appearance of slender SADs (because some structures were obscured by others due to the observational viewpoint) in this case further confirms the expectation. All in all, the moderate eruption, short ribbons, and tiled fan from the observed viewpoint result in small number of SADs in case D. The time period and the area of the fan in which we examine the SADs are indicated by the shaded region in Figure 22(c) and the red box in Figure 22(a), respectively.

B.5. Case E
On 2012 January 27, a flare of magnitude X1.7 erupted on the northwest of the Sun (case E). It started at about 17:00 UT, and the X-ray intensity achieved its peak at about 18:20 UT (Figure 23(c)). The decay phase lasted for several hours, during which dozens of SADs appeared. The STEREO-A observations show three main orientations of the arcade loops in this event. Group A (indicated by arrow "A" in Figure 23(b)) and group B (indicated by arrow "B" in Figure 23(b)) of the loops are almost perpendicular to each other, and the orientation of group C is in between the orientations of group A and group B. The orientation differences of the arcade loops lead to differences in the observation of the SAD features above the arcade loops by SDO/AIA 131 Å. As shown in Figure 23(a), the plasma above group A manifests in the images as a cusp, and no SADs were observed. There are a few SADs above group B and group C during the observation time, but the directions of motion of the SADs are different in the different groups. The SADs above group B appear to move almost parallel to the x-axis in Figure 23(a), and the direction of motion of the SADs above group C has a large angle of separation to the x-axis. The time period and the area of the fan in which we examine the SADs are indicated by the shaded region in Figure 23(c) and the red box in Figure 23(a), respectively. . It started at about 12:00 UT, and the X-ray intensity achieved its peak at about 13:25 UT (Figure 24(c)). The decay phase lasted for several hours, during which dozens of SADs appeared. The observations from STEREO-A show relatively long and single orientation ribbons in this flare, and the SADs appear in a broad spatial range along the arcades collectively observed from SDO/AIA 131 Å (Figure 24(a)). Some turbulent features were also seen from the SDO/AIA 131 Å observations (Figure 24(a)). The angular separation between the Earth and STEREO-A during the time period of the flare eruption is shown in Figure 24(b), and the time period and the area of the fan in which we examine the SADs are indicated by the shaded region in Figure 24(c) and the red box in Figure 24(a), respectively.