The morphology of average solar flare time profiles from observations of the Sun's lower atmosphere

We study the decay phase of solar flares in several spectral bands using a method based on that successfully applied to white light flares observed on an M4 dwarf. We selected and processed 102 events detected in the Sun-as-a-star flux obtained with SDO/AIA images in the 1600~{\AA} and 304~{\AA} channels and 54 events detected in the 1700~{\AA} channel. The main criterion for the selection of time profiles was a slow, continuous flux decay without significant new bursts. The obtained averaged time profiles were fitted with analytical templates, using different time intervals, that consisted of a combination of two independent exponents or a broken power law. The average flare profile observed in the 1700~{\AA} channel decayed more slowly than the average flare profile observed on the M4 dwarf. As the 1700~{\AA} emission is associated with a similar temperature to that usually ascribed to M dwarf flares, this implies that the M dwarf flare emission comes from a more dense layer than solar flare emission in the 1700~{\AA} band. The cooling processes in solar flares were best described by the two exponents model, fitted over the intervals t1=[0, 0.5]$t_{1/2}$ and t2=[3, 10]$t_{1/2}$ where $t_{1/2}$ is time taken for the profile to decay to half the maximum value. The broken power law model provided a good fit to the first decay phase, as it was able to account for the impact of chromospheric plasma evaporation, but it did not successfully fit the second decay phase.


INTRODUCTION
Flares are explosive events that occur in solar and stellar atmospheres. Flares are observed across a wide range of wavelengths, such as radio, visible, X-rays and gamma rays and the emission responsible for these observations originates from many different regions within solar and stellar atmospheres, from photospheres to coronae. It is generally believed that both solar and stellar flares are produced by the same mechanism, namely magnetic reconnection (Shibata & Magara 2011). However, discrepancies are readily observed, most notably the energy of the flares, which, for stellar flares, can often be several orders of magnitudes larger than even the largest solar flares.
While stellar flares have been studied for decades, it is only recently, through observations made by NASA's Kepler mission (Borucki et al. 2010), that a statistically large sample of flares for a single star has been obtained. Davenport et al. (2014) detected over ★ E-mail: A-M.Broomhall@warwick.ac.uk 6,000 flares on an M dwarf star in just 11 months of Kepler data. Using a subset, containing 885 flares, which were classified as "classical" because they contained just a single peak, Davenport et al. created an empirical flare template, that was well-represented by a polynomial in the fast rising phase, and two exponential decays, one each for the impulsive and gradual decay phases respectively. The first decay phase corresponds to radiative cooling losses and the second one is thought to be related to thermal conduction losses (Aschwanden 2004). Considering the large number of flares used, the scatter around this flare template is remarkably low.
Kepler obtained broadband white light observations (between 4200 and 9000 Å) of stellar intensity and has vastly increased the number of solar-like stars on which flares have been observed, prompting a number of statistical studies (e.g. Maehara et al. 2015;Notsu et al. 2019), many of which consider flaring rates, with a particular interest in determining how likely the Sun is to produce an energetic superflare. While flares are a common feature of Kepler light curves, the same cannot be said for white light flare observations of the Sun. Solar white light flares are rare because they tend to be relatively short in duration (a few minutes) and have a low contrast, making "Sun-as-a-star" observations of white light flares particularly challenging. Nevertheless, Kretzschmar (2011) performed a statistical study of Sun-as-a-star flares observed in a number of data sets, spanning a range of wavelengths, including Total Solar Irradiance (TSI). Kretzschmar demonstrated that white light emission is ubiquitous in solar flares, regardless of the energy of that flare, and that, while there is some dependence on flare strength, the blackbody temperature of solar flares is between 8, 000-10, 000 K, which is similar to estimates for M dwarf flare temperatures (Hawley & Fisher 1992). In a similar analysis to that of Davenport et al., Kretzschmar produced average flare profiles, using 2,100 flares, and found that the TSI profile only contained the impulsive phase, with no evidence of the gradual decay phase. Kretzschmar (2011) speculated that this could be because the gradual decay was below the noise level.
However, using resolved observations of the Sun, it has been shown that white light flares are common even for relatively weak flares (e.g. Matthews et al. 2003;Hudson et al. 2006;Jess et al. 2008). Kawate et al. (2016) use highly resolved white light observations of a single flare to produce a large number of light curves within the flaring region. The authors find that 58% of the light curves can be well represented by a single-component decay phase, while the other 42% are better represented by two components, akin to the morphology found by Davenport et al. (2014). Kawate et al. also found that, where two phases are favoured, the cooling times correspond to those expected for the chromosphere and corona respectively (Xu et al. 2006). Namekata et al. (2017) used the Helioseismic and Magnetic Imager (HMI) onboard the Solar Dynamics Observatory (SDO) to produce white light flare light curves, finding that, for a given energy, solar white light flares are an order of magnitude longer than stellar superflares observed with Kepler. Although, it is worth noting that the energies of the solar and stellar flares do not overlap, and so this result is based on extrapolations. Furthermore, the solar light curves were obtained by selecting spatial regions corresponding to the flare, rather than using Sun-as-a-star data. Although potential explanations for the discrepancy are proposed in terms of differences in cooling or magnetic field strength, the authors remain understandably cautious, stating that a better understanding of the mechanisms responsible for white light emission would help clarify the situation.
The exact emission mechanisms responsible for white light flares are not yet totally clear (see discussions in Kleint et al. 2016;Kawate et al. 2016, and references therein). It is generally assumed that continuum white light emission originates from the region encompassed by the mid-photosphere and lower chromosphere. White light flares can show good temporal correspondence with hard Xray emission (e.g. Fang & Ding 1995), implying that energetic electrons play an important role. While these electrons are capable of reaching the chromosphere, the mechanism by which the energy is transported to the photosphere is still debated. To determine how analogous solar flare morphology is to the white light stellar flares observed by Davenport et al. (2014) this article considers 304 Å, 1600 Å and 1700 Å data from the Atmospheric Image Assembly (AIA; Lemen et al. 2012) instrument onboard the SDO. Although these wavelengths are not traditionally considered as white light, and are outside the waveband observed by Kepler, the lines are used here because they originate from various locations within the photosphere and chromosphere (as described in Section 2).
Although TSI and Sun-as-a-star SDO/HMI data may be more akin to Kepler data, solar flares are infrequently observed in these data, even in the case of large eruptive events. As a result, analyses often rely on assumed theoretical models. For example, Emslie et al. (2012) had to complement the direct measurements of bolometric irradiance by estimations based on modelling. Moreover, the temporal resolution of SDO/AIA data is less than 1 minute, in contrast to TSI. This condition is important for the analysis of time profiles of solar flares which are more dynamic than stellar ones.
A better characterisation of the underlying shape of a flare permits studies of more transient flare light curve features, such as quasi-periodic pulsations (QPPs; see Van Doorsselaere et al. 2016;Kupriyanova et al. 2020, for recent reviews), which are notoriously difficult to detect and characterise robustly. Detection mechanisms that rely on detrending or model fitting (e.g. Dominique et al. 2018;Pascoe et al. 2017;Broomhall et al. 2019) would benefit from information concerning the underlying flare shape.
This article aims to investigate the solar-stellar flare connection, and emission mechanisms associated with white light flares by comparing the flare profile associated with various layers of the Sun's lower atmosphere with the flare profile associated with white light observations of flares on an M dwarf obtained by Davenport et al. (2014). The target of the study was to provide an instrument for the analysis of cooling during the decay phase and revealing cases related to additional sources of energy release. As previously mentioned, we use SDO/AIA data from the 1600 Å, 1700 Å and 304 Å channels. In Section 2, we describe the data in more detail, including how they were combined to determine median flare profiles for each channel. These median flare profiles are then fitted with both a combination of two exponential functions and a broken power law, as described in Section 3. The fits are discussed in detailed Section 4, and the main conclusions are summarised in Section 5.

OBSERVATIONS AND DATA PROCESSING
We obtained the total flux of the Sun-as-a-star using the images obtained by AIA (Lemen et al. 2012) on-board SDO. The channels used here demonstrate the chromospheric and photospheric emission (O'Dwyer et al. 2010): 1600 Å (transition region and upper photosphere), 1700 Å (temperature minimum, photosphere) and 304 Å (chromosphere, transition region). The image cadence is 24 s for the 1600 Å and 1700 Å channels and 12 s for the 304 Å channel.
The initial flare selection was performed using the GOES flare catalogue 1 . Each event was visually inspected and was required to have a "classical" flare time profile, consisting of a fast rise followed by a slow decay without any flattening or additional peaks in soft X-rays (SXR). When determining which flares to include in our sample the time interval onset was taken to be about 10 minutes before the flare maximum in GOES X-ray flux and the duration was defined as the length of time taken for the flux to decrease to the pre-flare level. More rigorous time scales were defined in the subsequent analysis, as described below. The initial list consisted of 359 flares from B5 to X9.3 GOES class.
After processing the AIA images with the standard package of SolarSoftware (Bentely & Freeland 1998;Freeland & Handy 1998, SSW), we obtained the total flux of the whole image for each time moment and produced a preliminary time profile of the flare. Thus, we processed the AIA images as if they were Sun-as-a-star or without an extracting the flare area and the resultant flux was equivalent to data of instruments observing without spatial resolution. To enable a combined analysis of the flares of different strength and duration, the profiles were normalised in both flux and time. First, each time profile was normalised using the flux maximum observed in each flare, meaning the maximum in normalised flux was unity for all flares in the sample. Second, we defined the time moment of the flux maximum as zero, so the time of the rise phase had negative values and the decay phase time had positive values. Finally, to get the same time scale for events of different duration, we presented the time series in normalised units defined by the time taken for the intensity to decrease to half the maximum for each time profile ( 1/2 ), following the methodology of Davenport et al. (2014). Then the duration of each time profile was limited to a range of -5 to 10 time bins ( 1/2 ). The time profiles in each channel were again checked to ensure there were no additional peaks during the decay phase that were missed when selecting the initial sample. The criteria for inclusion in the final sample was that any fluctuations in the time profile should be less than 30% of the flux value. We used the time profile derivative for control of this criterion. Application of this criterion also minimises uncertainties in the determination of the 1/2 value. As we processed solar flares from X to B GOES class and saturation in emission could not be excluded, this current criteria allowed us to escape significant contribution from saturation when forming time profiles.
The final sample contained 105 events in total from the 1600 Å and 304 Å channels (104 from 1600 Å and 102 from 304 Å). All profiles are presented on the same plot, with data shown as red crosses in Figures 1 and 3 for the 1600 Å and 304 Å channels, correspondingly. The X class flare ratio is 6%, the M class flare ratio is 58%, the C class flare ratio is 35% and B class flare ratio is 1%. Only 53 time profiles of emission in the 1700 Å channel showed both significant response and satisfied the all criteria outlined above (see Figure 2).
To construct the average time profile we re-sampled time profiles to a time resolution of 0.001 1/2 using linear interpolation, as  was done by Davenport et al. (2014). The average time profiles were defined as the median flux value computed at each time bin and are shown in Figures 1-3 as blue lines. One can see that the dispersion of initial data is significantly higher in the 1700 Å band plot, despite containing a smaller number of flare profiles.
The standard errors, as shown in Figures 1-3, are given by the median of values above and below the average value for high ( ℎ ) and low values ( ), correspondingly. This is also referred to as the interquartile range. As these errors are mainly asymmetric during the decay phase we used half of the sum of these two values ( = ( ℎ + )/2) as the error for the function fitting. The timeaveraged time profiles obtained for the solar flares demonstrate the similar behaviour to the M4 dwarf flare time profile obtained by Davenport et al. (2014) (see in Figure 4). During the decay phase, the evolution is fully agreed up to 1/2 equal about 1 . After this moment, the solar intensity decays slower in all three wavelength bands relative to the M4 dwarf flare emission. Gryciuk et al. (2017) determined an analytic template of the decay phase of a solar flare composed of a single exponent with various coefficients taking into account peculiarities of the flare profile shape. The template was derived and applied to soft X-ray observations that describe coronal plasmas. However, the decay phase of a solar flare is a composition of the cooling process and various heating processes, especially in lower layers of the solar atmosphere. As our time profiles clearly demonstrated a steep initial decay followed by a slower decline in flux, we decided to apply the two-phase model. This model was successfully used by Davenport et al. (2014) for the flares observed on an M4 dwarf. Despite the difference in the atmospheric structure (temperature and density stratification) between the Sun and M dwarfs, the flare evolution is similar in both cases (Allred et al. 2006).

DECAY PHASE FITTING
The average time profiles, obtained using the natural logarithm of flux, are presented in Figures 5-7. The flux data were fitted using two straight lines, 1 ( ) = 1 exp 1 1/2 and 2 ( ) = 2 exp 2 1/2 , corresponding to the two different decay phases. Initially, the same time regions as Davenport et al. (2014) were used for the fitting. These ranges are 0 < 1/2 < 0.5 for 1 ( ) and 3 < 1/2 < 6 for 2 ( ). The results are shown in red in Figures 5-7. In addition, alternative time regions were chosen to improve the fitting of each phase: 0 < 1/2 < 1.5 for 1 ( ) (shown in blue) and 3 < 1/2 < 10  for 2 ( ) (shown in green). The results of the fitting are shown in Table 1. We also fitted the time profiles with a broken power-law model to reveal the time of the phase change and with the aim of finding a better description of the average time profiles. The function was: where break , index 1 and index 2 are all parameters determined by the fitting process. Table 2 shows the results of fitting equation 1 to the median flare profiles, with 2 values. A comparison of fitting a broken power-law model with the previous fitting, of two independent exponents, is presented in Figure 8. One can see that the broken power-law model describing the decay phase of the average time profiles provides a good fit to the data until = 6 1/2 . We note that the broken power-law function is a better fit to the flare profile at around 2 1/2 .

DISCUSSION
In this study we obtained averaged time profiles and fitted templates describing the decay phase of a solar flare for different spectral bands. As previously mentioned, there are two main mechanisms determining the cooling processes of a flare decay: radiative cooling and thermal conduction (Aschwanden 2004). According to the model of plasma cooling, thermal conduction losses initially dominates during the first part of the decay phase. Later we observe the domination of radiation losses during the second phase (see, Cargill et al. 1995;Aschwanden & Tsiklauri 2009). Cargill et al. (1995) also noted that if radiative losses dominate during the initial phase of decay, they still dominate during the entire decaying phase. However, here we can discriminate between the steep (called impulsive by Davenport et al. 2014) and gradual phases during the decay phase of the analysed time profiles. For this reason, we focused our study on the fitting of the averaged time profiles with models consisting of two components only, as done in previous studies (see e.g. Davenport et al. 2014). Domination of thermal conduction or radiative losses relates to the cooling times due to these processes, and depends on temperature, density and loop length in the case of thermal conduction. Thus, the flux behaviour during the decay phase depends on two factors -the ratio between cooling by radiation and cooling by thermal conduction and the ratio between the temperature and density of the generating emission region. The temperature of spectral lines that dominate in the formation of emission of the chosen spectral bands depends on the formation height of the emission. According to O'Dwyer et al. (2010), emission in the 304 Å band is mainly formed by the doublet H II line (303.78Å) that has log( ) = 4.7 and corresponds to the chromosphere and transition regions. We note that the SDO/AIA 304Å band also has a contribution from two spectral lines with log( ) > 6 (see O'Dwyer et al. 2010, and references in it). The first line is the Si XI 303.33 Å line with log( ) = 6.2 whose contribution is significant in quiescent emission of off-limb active regions only. The Ca XVIII 302.19 A line with log( ) = 6.85 contributes to flare emission. However, its fraction of total emission was estimated as 0.05.
The quiet Sun emission in the 1600 Å band is generated predominantly by the C IV line and continuum emission, both with log( ) = 5. This temperature value implies that the level where the emission of this band is formed is higher than the formation height of the 304 Å band. The continuum dominates the emission of the 1700 Å band for quiet Sun regions. Thus, the emission is formed at log( ) = 3.7, and it typically originates from the photosphere. The results obtained by Simões et al. (2019) refined the contribution of the spectral lines emitted within the 1600 Å and 1700 Å bands. The authors confirmed the share of the C IV doublet emission to the 1600 Å band, but they also revealed the domination of the C I 1656 Å multiplet contribution over the continuum for the 1700 Å band. This means that emission in the 1700 Å band is mostly formed at log( ) = 4-4.2. Therefore, the emission from this band should be the closest in temperature to M4 dwarf flares, where log( ) ≈ 4 (Hawley & Fisher 1992).
However, we should take into account that the structure of the atmosphere of an M dwarf star differs from the solar atmosphere. The temperature in the M dwarf atmosphere rises from cooler than the solar photosphere to hotter than the solar corona within a more narrow height range (see, Allred et al. 2006). Moreover, the initial energy release should occur lower in an M dwarf atmosphere than in the solar one (see Figure 8a, Allred et al. 2006).
If we fit the first phase of cooling using the same time interval as Davenport et al. (2014), the coefficient 1 of decay for solar flux is lower than that obtained by Davenport et al. for the M4 dwarf even in the case of the 1700 Å band (0.800 vs 0.965). Increasing the fitting interval for the first cooling phase up to 1.5 1/2 decreased the coefficient 1 to 0.668 for the 1700 Å band. We note that the difference between the coefficients obtained for different time intervals exceeds the error bars. Therefore, cooling during the first steep phase was found to be even slower for solar flares than for flares on the M4 dwarf. As the 1700 Å channel and M dwarf flares are associated with similar temperatures this implies that the plasma responsible for M4 dwarf flare emission is denser than that associated with the 1700 Å solar flares, and so potentially originates from a deeper layer. The fact that the coefficient decreases when the fitting range is extended could be related to the impact of chromospheric evaporation, which would cause an increase in flux at around 1-3 1/2 . Fitting a broken power law to the 1700 Å band gave the decay coefficient index 1 equal to 0.714 and the time of the phase change as about 1.3 1/2 . This implies that a new phase with a different rate of cooling onsets after 1.3 1/2 . The coefficient index 2 for this phase, according to the fitting by a broken power law, is about 0.2. This value is close to the value obtained for M4 dwarf flare (0.29) but it does not fit time profile values above 5 1/2 . This could be due to the impact of the evaporation of the heated chromospheric plasma to the flare time profile (see, for example, Fletcher et al. 2011) . As predicted by the modelling of Allred et al. (2006) for the 304 Å band, the contribution of the evaporation for solar flare flux is more significant than for M dwarfs. Figure 9 contains a comparison of the 304 Å band and the appropriate model from Allred et al. (2006), where we can see the model flux drop below the observed flux at around = 1/2 . Then the impact of chromospheric evaporation is seen in the form of an increase in the modelled flux. Although we note the model flux again drops below the observed flux for > 2 1/2 .
In the case of fitting for t1=[0, 0.5] 1/2 , the relation between the decay coefficients of the first phase for the 1600 Å, 304 Å and 1700 Å bands confirms that the temperature of the emission formation is as discussed above. The value of break in the broken power-law fitting characterises the transfer from the radiative cooling to the contribution of chromospheric evaporation. We obtained the lowest value of break for the 1700 Å band (1.282) and the highest value for the 1600 Å band. This is the same ordering as suggested by the temperature-height model of the solar flare atmosphere. An alternative way to specify a temperature-height model is through an analysis of the delay between maximums of flare emission in different spectral ranges. The delay between signals obtained in the UV and EUV bands are actively used for analysis of heighttemperature dependence of solar atmosphere over the sunspot (see, Reznikova et al. 2012). As one can see in Table A1, the delay between the maximum of flares in 304Å band and 1600Å band can be both positive and negative. The most common time delay was 17 seconds but a group of events demonstrated a negative delay of about 7 seconds. The uncertainties related to an image cadence of 24 seconds are 14 seconds (using the same approach as Reznikova et al. (2012)). Thus we can conclude that negative delays are within error bars, and most of the positive delays are around this value. Such small delays, close to the level of uncertainties, agree with classical F1 and F2 flare models by Machado et al. (1980) suggesting a height difference between the layers with temperatures corresponding to 304Å and 1600Å bands of about 3 km.
When fitting the two exponents to the median flare profiles, we obtained two sets of parameters for each of the first and second phases by fitting over different time intervals (see Table 1). The difference between the parameters obtained for the same phase, but for the different time intervals, exceeds the parameter uncertainty. Thus, it is necessary to determine which time intervals should be fitted to best represent the cooling processes. We believe that the decay during the second phase should be fitted using parameters obtained for t2 = [3, 10] 1/2 because the flux takes longer to decay in the solar atmosphere compared to the M dwarf flare profile. To reveal the decaying related to radiative and thermal conduction losses during the first phase, we should avoid the chromospheric evaporation contribution to emission. As mentioned above, the breakpoint, break in the broken power-law fitting parameters characterises the transfer from the radiative cooling to the contribution of the chromospheric evaporation. For both the 1700 Å and 304 Å bands, the obtained values of break were below 1.5. While, for these bands, the parameters obtained by fitting the different time intervals showed a significant difference, the parameters obtained by fitting the different time intervals for the 1600 Å band, where break was above 1.5, did not demonstrate a significant difference. All these facts indicate the impact of chromospheric evaporation to time profiles is not negligible (depending on observational wavelength), meaning that fitting over the interval t1 = [0, 1.5] 1/2 is not appropriate and the range t1 = [0, 0.5] 1/2 is favourable.
Based on a comparative analysis of the fitting results for solar and M dwarf flare templates, we can conclude that, for the Sun, the optimal template describing cooling processes consists of two exponents, fitted for t1 = [0, 0.5] 1/2 and t2 = [3, 10] 1/2 respectively.
We note that the second cooling phase of solar flares turned out to be more complicated than for M4 dwarf flares. The cooling of the 304 Å band during the second phase is marginally faster than for the 1600 Å band, which is hotter and originates from less dense plasma. It is possible that, during cooling the relative contribution of spectral lines with different temperatures to the emission of the 1600 Å band changes with time, which results in a faster decrease of emission flux. We also would like to note that the 304 Å band time profile demonstrated unusual behaviour. As emission of this band mostly forms by the emission of a single spectral line (O'Dwyer et al. 2010), its time profile should be more akin to density evolution (Aschwanden & Tsiklauri 2009). This fact is confirmed by modelling (see Allred et al. 2006). However, the observed time profile shows an exponential decay, which is more characteristic of temperature evolution.
As mentioned above, there are two processes, namely thermal conduction or radiation losses, and domination of one of the process over the other would result in different behaviours of temperature. The temperature evolution during radiative cooling can be described as ( , ) = 0 ( )[1 − (1 − ) / 0 ] 1/(1− ) , where 0 is the radiative cooling time at the start of the radiative phase, is the coefficient of radiative losses function and is the coordinate along the magnetic field (here and after, Cargill et al. 1995). The coefficient is assumed to equal to -0.5 for log(T)>5. As the temperature range of the data analysed here is less this value, we cannot use this assumption. The dependence ( ) = 0 (1 + / 0 ) −2/5 describes the temperature evolution for the case of static conductive cooling (where 0 is the conductive cooling at the beginning). We performed fitting using these two functions for the analysed spectral bands, using 0 , 0 and (1 − ) as variables. The results can be seen in Figure 10. The averaged time profiles for all bands are better fitted by the function corresponding to radiative cooling. Only the 304Å band time profile demonstrates good agreement with the fitting based on conductive cooling for a short time interval when < 0.1 1/2 .

CONCLUSIONS
Based on the results of reconstruction and analysis of averaged time profiles of solar flare decay phases obtained within 304 Å, 1600 Å and 1700 Å spectral bands, we can conclude the following: • The processes of cooling during the decay phase of solar flares are the most closely described by fitting of a combination of two independent exponents for t1 = [0, 0.5] 1/2 and t2 = [3, 10] 1/2 .
• The parameters obtained for the 1700 Å solar flare time profile, which is theoretically closest in temperature to the M dwarf flares, are consistent with models that imply that the white light M dwarf flare emission is formed in a higher density of layer.
• Fitting a broken power-law model allows the contribution of chromospheric evaporation to be taken into account. However, it is not sufficient for a full description of the second part of the decay phase. Table A1 contains details of the flares in our sample.

DATA AVAILABILITY
The datasets were derived from sources in the public domain: SDO/AIA data was obtained from Joint Science Operations Center (JSOC) (http://jsoc.stanford.edu); GOES data were obtained from the GOES Flare Catalogue (https://hesperia.gsfc.nasa.gov/goes/goes_event_listings/). This paper has been typeset from a T E X/L A T E X file prepared by the author.