The Lyman-alpha Emission in Solar Flares. I. a Statistical Study on Its Relationship with the 1--8 \AA\ Soft X-ray Emission

We statistically study the relationship between the Lyman-alpha (\lya) and 1--8 \AA\ soft X-ray (SXR) emissions from 658 M- and X-class solar flares observed by the {\em Geostationary Operational Environmental Satellite} during 2006--2016. Based on the peak times of the two waveband emissions, we divide the flares into three types. Type I (III) has an earlier (a later) peak time in the \lya\ emission than that in the SXR emission, while type II has nearly a same peak time (within the time resolution of 10 s) between the \lya\ and SXR emissions. In these 658 flares, we find that there are 505 (76.8\%) type I flares, 10 (1.5\%) type II flares, and 143 (21.7\%) type III flares, and that the three types appear to have no dependence on the flare duration, flare location, or solar cycle. Besides the main peak, the \lya\ emission of the three type flares also shows sub-peaks which can appear in the impulsive or gradual phase of the flare. It is found that the main-peak (for type I) and sub-peak (for type III) emissions of \lya\ that appear in the impulsive phase follow the Neupert effect in general. This indicates that such \lya\ emissions are related to the nonthermal electron beam heating. While the main-peak (for type III) and sub-peak (for type I) emissions of \lya\ that appear in the gradual phase are supposed to be primarily contributed by the thermal plasma that cools down.


INTRODUCTION
Solar flares are sudden brightenings and energetic events in the solar atmosphere (e.g., Fletcher et al. 2011;Shibata & Magara 2011). In the standard flare model (also called CSHKP model ;Carmichael 1964;Sturrock 1966;Hirayama 1974;Kopp & Pneuman 1976), the magnetic energy is released in the corona through magnetic reconnection. The energy is then transported downward to the chromosphere via nonthermal electron beams and/or thermal conduction. Consequently, the chromospheric plasma is heated and goes up into the corona due to a thermal pressure (i.e., chromospheric evaporation). The evaporated hot plasma fills the flare loops that can be well seen in the soft X-ray (SXR) waveband. Based on the observed SXR light curve, the evolution of a flare can be classified into three phases, pre-flare phase, impulsive (or rise) phase, and gradual (or decay) phase (e.g., Fletcher et al. 2011;Hudson 2011). The SXR emission has begun to increase in the pre-flare phase; then it has a rapid rise, reaching its maximum in the impulsive phase because of a significant heating; in the gradual phase, it decays gradually due to the plasma cooling. The heating and cooling processes can also be manifested in some other wavebands. For example, the nonthermal electron beam heating produces a hard X-ray (HXR) emission that resembles the time derivative of the SXR emission (i.e., Neupert effect; Neupert 1968). While the hot plasma that is heated earlier will cools down gradually and be observed in the SXR and extreme-ultraviolet (EUV) wavebands successively.
In the solar UV spectrum, there is a strongest line of hydrogen Lyman-alpha (H I Lyα) at 1216 A (Curdt et al. 2001). This line is formed in the mid-to-upper chromosphere and low transition region (Vernazza et al. 1981), which is optically thick and suffers an opacity effect (e.g., Vial 1982;Woods et al. 1995). It has been reported that the Lyα line shows a significant enhancement during solar flares (Woods et al. 2004;Rubio da Costa et al. 2009;Milligan et al. 2012Milligan et al. , 2014Kretzschmar 2015;Milligan et al. 2020). However, up to now, we still have a poor understanding on this line and its physical origin in solar flares owing to relatively rare observations as well as a complex formation process of this line.
In the past decades, there have been a few studies on Lyα in solar flares by using spectral/imaging/photometric observations. Canfield & van Hoosier (1980) presented the Lyα line profiles and their temporal variations for two flares observed by a slit spectrograph on Skylab. With the L.P.S.P. experiment on OSO-8, Lemaire et al. (1984) preformed a simultaneous observation of a flare in six chromospheric lines including Lyα and found an indication of a downward energy propagation via the temporal behaviour of the different lines. The Lyα emission was found to increase by about 0.6-30% during flares (Brekke et al. 1996;Woods et al. 2004;Kretzschmar et al. 2013;Milligan et al. 2020) with its radiated energy rate estimated to be ∼10 25−27 erg s −1 (Johnson et al. 2011;Milligan et al. 2012Milligan et al. , 2014. In particular, Milligan et al. (2020) carried out a statistical study based on 477 M-and X-class flares observed by the Geostationary Operational Environmental Satellite (GOES) and showed that 95% of these major flares have an enhancement of 10% or less in Lyα, with a maximum of about 30% in all of the flares. The authors also reported that there is a center-to-limb variation in the Lyα emission due to an opacity effect. Besides Lyα enhancements, quasi-periodic behaviours at one or three minutes were detected in the Lyα emission during flares (Milligan et al. 2017;Li et al. 2020). In addition, there were some studies focusing on the relationship of the Lyα emission with soft X-ray (SXR) and hard X-ray (HXR) emissions. Nusinov et al. (2006) reported that the Lyα maximum is reached well before the SXR maximum and the variation of the Lyα emission is synchronous with that of the HXR emission above 50 keV during the impulsive phase of a flare. These "Neupert-effect" features were also observed and confirmed by some other authors via studying different flare events (Rubio da Costa et al. 2009;Milligan & Chamberlin 2016;Milligan et al. 2017;Dominique et al. 2018;Chamberlin et al. 2018). In particular, with the Lyα imaging observations from the Transition Region and Coronal Explorer (TRACE), Rubio da Costa et al. (2009) demonstrated that most of the Lyα emission is co-spatial with the HXR sources which originate from the flare footpoints. These observational results indicate that the Lyα emission in flares has a nonthermal origin.
Some modelling works on flaring Lyα have also been emerging in quite recent years. Using the radiative hydrodynamics code RADYN, Brown et al. (2018) simulated the response of the Lyα line with hard and soft electron beam heatings and compared the synthetic Lyα line profiles with the observed ones from the Extreme-ultraviolet Variability Experiment (EVE) onboard the Solar Dynamics Observatory (SDO). Via RADYN, Hong et al. (2019) also calculated the Lyα line profiles in nonthermal and thermal heating models in which the Lyα line exhibits different features in line asymmetry and light curves. Moreover, Yang et al. (2020) modelled the Lyα and Hα lines as well as the continua at 3600Å and 4250Å in white-light flares and found that the Lyα line has different responses to the nonthermal beam heating with the Hα line and the continua. Besides RADYN, the HYDRO2GEN code was used to calculate the Lyman lines and continuum with a beam injection model (Druett & Zharkova 2019). These simulations show that the Lyα line has a notable response to the nonthermal and also thermal heatings.
In order to study the physical origin of the Lyα emission in solar flares, here we statistically study the relationship between the Lyα and SXR emissions from 658 M-and X-class flares observed by GOES during 2006-2016. We find that in most of the flares, the Lyα emission peaks earlier than the SXR emission and holds the Neupert effect in general, indicative of a nonthermal origin, just as previous studies demonstrated. However, we also find that there are about one fifth flares with their Lyα peak later than the SXR peak, which may primarily have a thermal origin. To the best of our knowledge, the latter result has never been reported in observations before, as least from a statistical aspect. The statistical study in the present work can improve our understanding on the physical origin of the Lyα emission in solar flares. In the following, we describe the instruments and data in Section 2 and our flare dataset in Section 3. The statistical results are shown in Section 4, followed by a summary and interpretation in Section 5. In the last Section 6, we give the conclusion and discussions.

INSTRUMENTS AND DATA
The GOES series spacecraft (GOES-1-17) have been providing the solar X-ray irradiance continuously since 1975. The X-ray Sensor (XRS; Hanser & Sellers 1996) on GOES observes the full-disk Sun in two soft X-ray channels, i.e., 0.5-4Å (short channel) from XRS-A and 1-8Å (long channel) from XRS-B, the latter of which is widely used to define the flare magnitude from A-to X-class. Since GOES-13 (launched in 2006), there also includes an Extreme Ultraviolet Sensor (EUVS; Viereck et al. 2007) that observes the whole Sun in the EUV wavebands. EUVS contains five channels, called A, B, C, D, E, from 50-1270Å with the E channel targeting the Lyα emission. The E channel has a width of ∼90Å covering a wavelength range of 1180-1270Å, whose emission primarily comes from the Lyα line.
In this work, we mainly use the Lyα and 1-8Å SXR data. These data have been converted into irradiances in units of W m −2 , which have cadences of ∼10 s and 2 s, respectively. The Lyα data that have been publicly released are for the period of 2006-2016 (the GOES-13-15 era), namely covering an entire solar cycle. Considering that some Lyα data for a certain time were simultaneously obtained by two or even three spacecraft, here we choose the data on the following basis: from GOES-15 first, if no then from GOES-14, and finally from GOES-13. In order to obtain the flare-induced emission in both of the Lyα and SXR wavebands, we have made a background subtraction by averaging the pre-flare flux from one hour before the flare onset. Note that a correction (divided by 0.7) for the "true" SXR flux is not applied in this study, which is followed by the work in Milligan et al. (2020). In addition to the Lyα and SXR data from GOES, we use the imaging data from the Atmospheric Imaging Assembly (AIA; Lemen et al. 2012) on board SDO and the sunspot data 1 from the World Data Center SILSO at the Royal Observatory of Belgium in Brussels. Duration of the flares. The distribution histogram of the flare duration is given in the top left panel of Figure 1, which is calculated from the start and end times recorded in the GOES flare list. Note that there are 13 (∼2%) flares that last for >120 minutes and are randomly distributed in the time range of 120-720 minutes and we do not show them in this histogram just for a better display. It is seen that the flare dataset contains both short-and long-duration events that last for a few minutes to hours, though most of the flare events end within 30 minutes. Note that GOES only records part of the decay phase of the flare, which actually underestimates the flare duration. It is also seen that the flare duration has no dependence on the flare magnitude, namely both short-and long-duration flares can be M-or X-class flares.
Location of the flares. We plot the distribution histogram of the flare location in the top right panel of Figure 1. The location is represented by the distance of the flare from disk center that is calculated from the X and Y coordinates 2 . Note that there are 19 (∼3%) M-class flares that lack the location information and thus are not included here. One can see that the flare examples in our dataset are located through the disk center (i.e., the distance is 0 ′′ ) to the limb (∼960 ′′ ). Here it should be clarified that the increasing trend of the distribution is mainly caused by a foreshortening effect when the solar hemisphere is projected to the X-Y plane. This trend does not reflect the real distribution of the flare on the solar surface. In fact, according to the flare distribution as a function of heliocentric angle shown in Milligan et al. (2020), there exists a nearly uniform distribution of flares from disk center to the limb. The approximately uniform distribution can also be seen in the top right panel of Figure 4 to some extent. Within the circle of a radius of 500 ′′ (or the heliocentric angle being <30 degrees more or less), where the foreshortening effect is somewhat trivial, we can see that the flare is in general distributed uniformly on the plane. From the histogram of flare location, one can also see that the location distribution has no dependence on the flare magnitude either, i.e., both X-and M-class flares can take place on the disk or at the limb.
Peak fluxes of the flares. In the bottom panels of Figure 1, we show the scatter plots of the peak Lyα flux versus peak SXR flux for the flare examples, with the left one without subtracting the background and the right one with the background subtracted. It can be seen that there seems to be a correlation between the Lyα and SXR peak fluxes whether the background is subtracted or not. We calculate the correlation coefficients for the two cases, which are 0.63 and 0.45, respectively. It should be mentioned that the pattern with subtracting the background is actually similar to the one as reported in Milligan et al. (2020), though we use a different method to find the background level. In the remaining paper, we just adopt the peak fluxes of both Lyα and SXR that have been subtracted the background emission.

RESULTS
After viewing the Lyα and SXR emission curves of all the flares in the dataset, we see some features as follows. (1) There is usually a rapid rise followed by a relatively slow decay in both of the Lyα and SXR emission curves. (2) The Lyα emission can reach its maximum (main peak) earlier or later than the SXR emission. (3) There appear two or even more evident peaks (main peak plus sub-peaks) in numerous flares particularly in the Lyα emission curve. These features can be seen from the light curves as plotted in Figures 2, 5, and 6.

Three Types of Flares based on the Main Peak of the Lyα emission 4.1.1. Example Light Curves
According to the time sequence of the main peaks of Lyα and SXR emissions, we divide the flares into three types: type I (III) has an earlier (a later) peak time in the Lyα emission than that in the SXR emission, while type II has nearly a same peak time (within the time resolution, or ±10 s) between the Lyα and SXR emissions. In other words, the Lyα emission in type I/III flares peaks in the impulsive/gradual phase of the flare, while the one in type II flares peaks around the same time as the SXR flux. These Lyα emissions in different types are expected to be related to different physical processes according to the standard flare model. Figure 2 gives three examples for the three types of flares. The top panel shows an example for a type I flare. One can see that both of the Lyα and SXR emissions exhibit a good single peak and the Lyα emission peaks more than one minute earlier than the SXR emission. Here we also plot the time derivative of the SXR flux that actually corresponds to the Lyα emission very well. This indicates that the Lyα emission in this flare follows a Neupert effect, which has already been reported before though for different flares. The middle panel displays an example for a type II flare. It is seen that the Lyα emission reaches its maximum nearly at the same time as the SXR emission with a time difference of 8 s. In this case, the Lyα emission obviously peaks later than the time derivative of the SXR emission. An example for a type III flare is shown in the bottom panel. We can see that the Lyα emission rises slowly and reaches its maximum about three minutes later than the SXR emission. It should be highlighted that this case has rarely been reported before.

Time Difference between the Lyα and SXR Emission Peaks
It is found that among the 658 flares in our dataset, there are 505 (76.8%) type I flares that consist of 35 X-class and 470 M-class flares, 10 (1.5%) type II flares all of which are M-class flares, and 143 (21.7%) type III flares that consist of 8 X-class and 135 M-class flares. Here we compute the time difference (denoted by t p ) between the Lyα and SXR emission peaks for all the three type flares and show the results in Figure 3. It should be mentioned that we only consider the Lyα peak that falls into the range between the start and end times of each flare. This would actually underestimate the t p value for some type III flares that have a Lyα peak later than the flare end time (see the two flares as shown in the bottom panels of Figure 6). After carefully checking the light curves, we find that this kind of type III flares are quite few and do not substantially affect the statistical results here. According to our definition, type I flares have negative values of t p and type III flares have positive values of t p . Considering that the majority of flares (more than two thirds) have a value of t p within five minutes, here we only show the t p histogram in the range of ±300 s in the top panel of Figure  3. It is seen that both types I and III flares include M-and X-class events, as already revealed from above. In addition, most of the type I flares have a t p within three minutes, while the t p of type III flares is distributed equally within five minutes in general. The middle and bottom panels of Figure  3 show the scatter plots of t p versus peak SXR and Lyα fluxes, respectively. It can be seen that there seems to be no obvious relationship between t p and the peak fluxes of Lyα or SXR.

Dependence on Flare Duration, Flare Location, and Solar Cycle
We also check the dependence of flare types, more accurately speaking, of types I and III owing to quite a small number of type II flares, on the flare duration, location, and solar cycle, which is shown in Figure 4. From the top left panel we see that the flares with different types can be short-and long-duration flares, namely the flare type has no dependence on the flare duration. We can also see that the flare type is independent on the flare magnitude, either. From the spatial distribution of the flares with different types on the solar disk (the top right panel of Figure 4), one can see that there is no correlation between the flare type and the location. In addition, we plot the count variation of different type flares with years or solar cycle in the bottom panel of Figure 4, also overplotted the curve for sunspot counts. It is seen that the counts of flares with different types generally match the sunspot counts over the solar cycle.

Representative Light Curves
Apart from main peak, the Lyα emission curve shows some evident sub-peaks in a large number of flares (roughly two thirds checked by eye), either for type I or type III. Such kind of light curves are plotted in Figures 5 and 6. From Figure 5 one can see that the Lyα emission shows up one or two evident peaks, i.e., main peak and/or sub-peak, in the impulsive phase in the four type I flares, which resemble the time derivative of the SXR flux. Moreover, additional sub-peak appears in the gradual phase, with a time delay relative to the SXR peak. For the type III flares as shown in the top panels of Figure 6, the Lyα emission exhibits a sub-peak in the impulsive phase which also matches the time derivative of the SXR flux, in addition to the main peak that appears in the gradual phase. Note that for the two type III flares in the bottom panels of Figure 6, the Lyα emission does not show any significant enhancements in the impulsive phase but displays a few peaks in the gradual phase. In particular, a stronger peak appears in the late gradual phase. We notice that these two flares are located at or close to the solar limb (see the coordinates in the right corner of the panel), whose loop footpoints may be occulted (or partly occulted) by the solar disk. This would cause a weak Lyα emission increase during the impulsive phase of the flare (Milligan et al. 2020). While the significant enhancement of Lyα in the late phase could be contributed by the flare loops that are still visible above the solar limb. Generally speaking, there usually exhibit an impulsive-phase peak (i.e., main peak for type I flares but sub-peak for type III flares) and a gradual-phase peak (main peak for type III flares but sub-peak for type I flares) in the Lyα emission. Here it should be mentioned that there are also about one third flares which only show an isolated peak (or main peak) in Lyα during the flare period. This case will be discussed in Section 5.

Neupert Effect Check for the Impulsive-phase Peak of Lyα
For the impulsive-phase peak emission of Lyα, it is worthwhile to check its validity of the Neupert effect. First, we take a look at the main peak of Lyα for all the 505 type I flares. The top panel of Figure 7 shows the histogram of the time difference (denoted as t d ) between the Lyα main peak in the impulsive phase and the peak of the time derivative of the SXR flux. We can see that t d , either for X-or M-class flares, has a notable peak around zero and the majority of flares (∼80%) have a t d value within ±2 minutes. This indicates that the main peak emission of Lyα in type I flares holds the Neupert effect generally. The middle and bottom panels of Figure 7 display the t d distribution versus the peak SXR and Lyα fluxes, respectively. It is seen that larger flares have a smaller t d value, which may suggest that larger flares hold the Neupert effect better. However, there should be a caution that this might be owing to a data bias, since we have much fewer X-class flare examples than the M-class ones in our dataset. Regarding t d versus the peak Lyα flux, it seems that there is no correlation between them.
Furthermore, we identify an evident impulsive-phase peak (sub-peak) of Lyα in 62 type III flares by eye and add their t d distribution on the result of type I flares (see Figure 8). It is seen that similar to the main peak of type I flares, the impulsive-phase peak of Lyα in type III flares also holds the Neupert effect mostly. Note that for the remaining 81 type III flares, it is somewhat hard to identify a well isolated impulsive-phase peak of Lyα by eye owing to its very weak emission or too many small bumps in the Lyα light curve.

Delayed Time Check for the Gradual-phase Peak of Lyα
We also study the delayed time relative to the SXR peak for the gradual-phase peak of Lyα and its relationship with the impulsive-phase peak flux of Lyα and the flare loop length which are expected to be related to the initial heating and loop cooling (e.g., Cargill et al. 1995;Yoshida & Tsuneta 1996). Here we firstly identify a good gradual-phase peak (sub-peak) of Lyα in 173 type I flares by eye and plot their delayed time (denoted as t ′ d ) versus the impulsive-phase peak (main peak) flux of Lyα in the top panel of Figure 9. Furthermore, we over-plot the result for the 62 type III flares (delayed time of the main peak, or their t p , versus the flux of the sub-peak of Lyα in the impulsive phase) that are mentioned above. It is seen that there shows a cone shape or a weak correlation with a coefficient of 0.1 between the two parameters. Or we could say that as the impulsive-phase peak flux of Lyα increases, the delayed time of the gradual-phase peak will more likely be larger.
In the bottom panel of Figure 9, we give the relationship of the delayed time for the gradual-phase peak of Lyα with the flare loop length. Note that here we only select 35 flares (24 out of 173 type I flares plus 11 out of 62 type III flares) that are located within the circle of a radius of 500 ′′ (see the top right panel of Figure 4) to show this relationship. In this circle, the flare loop length could be measured more accurately by using AIA images due to the projection effect. Here we derive the loop length mainly from the AIA 171Å images combining with the AIA 1600Å images for the loop footpoint identification. Note that the horizontal error bar for the loop length comes from different sets of flare loops with different lengths in a single flare. From the selected 35 flares we can see that there appears to be a correlation, with a coefficient of 0.5, between the delayed time of the gradual-phase peak of Lyα and the loop length.
To further check the relationship between the delayed time of the gradual-phase peak of Lyα and the impulsive-phase peak flux of Lyα or the flare loop length, we carry out very simplified simulations for the selected 35 flares through the "enthalpy-based thermal evolution of loops" (EBTEL) model (Klimchuk et al. 2008;Cargill et al. 2012a,b). Here we treat each of the flares with a single-loop flare and the heating pulse assumed to be a Gaussian shape is constrained from the Lyα light curve. More specifically, the peak heating rate (denoted as Q p ) and the heating duration (or the Gaussian width, σ) are proportional to the impulsive-phase peak flux of Lyα and the flare rise time (i.e., from the flare onset time to the SXR peak time), respectively. We give an example for the heating pulse in the top left panel of Figure 10. The peak heating rate and the heating duration for all the 35 flares are shown in the bottom left panel, as a function of the loop length. It is seen that in our modeling, Q p (diamond symbols) is in the range of 0-2.5 erg cm −3 s −1 and σ (plus symbols) is from 0-180 s. Note that here we set a maximum of σ = 180 s in the single loop for a few flares that have a relatively longer rise time, according to some flare heating modelings via EBTEL (Raftery et al. 2009;Li et al. 2012Li et al. , 2014a. The top right panel of Figure 10 shows an example for the evolution of average temperature (T ) and electron density (n) of a single loop. Using the average temperature, electron density, and loop length, we compute the cooling time (t cool ) and show its relationship with the loop length in the bottom right panel of Figure 10. Here t cool is defined as 1/t cool = 1/t c + 1/t r , where t c and t r are conductive and radiative cooling times, respectively (Cargill et al. 1995;Raftery et al. 2009;Cargill et al. 2012b;Sun et al. 2013). One can see that as the loop length increases, the cooling time becomes longer, which is consistent with the relationship between the delayed time of the gradualphase peak of Lyα and the flare loop length in observations as shown in the bottom panel of Figure  9. This suggests that the delayed time is a good indicator of plasma cooing, which is related to the flare heating magnitude and especially the flare loop length.

SUMMARY AND INTERPRETATION
The statistical results as well as light curve features for Lyα from our flare dataset are summarized in Tables 1 and 2. According to the fact that the flare types have no dependence on the flare magnitude, duration, location, or peak fluxes of Lyα or SXR, we conjecture that the Lyα emission is likely contributed by some common processes in each of the flare, say flare heating and plasma cooling. In the following, we attempt to more specifically interpret the Lyα emission mainly for type I and III flares based on its relationship with the SXR emission in the framework of the standard flare model (e.g., Raftery et al. 2009;Holman et al. 2011;Hudson 2011).
The impulsive-phase peak (i.e., main peak for type I flares and sub-peak for type III flares) of Lyα basically holds the Neupert effect, which suggests that this kind of emission is closely related to the nonthermal electron beams (Nusinov et al. 2006;Rubio da Costa et al. 2009;Milligan & Chamberlin 2016). It is known that in the standard flare model, the electron beams deposit their energy mostly in the chromosphere via Coulomb collisions, where the Lyα line is mainly formed. The local plasma is heated and then we can see evident emission in the Lyα waveband together with some other bands. In particular, the Lyα emission in this case is supposed to come from the chromospheric footpoints of flare loops, which has been supported by TRACE imaging observations (Rubio da Costa et al. 2009). This nonthermal origin of Lyα is also confirmed in radiative hydrodynamic simulations. With an electron beam heating model, Yang et al. (2020) calculated the time evolution of the Lyα line intensity and found that the Lyα intensity peaks nearly at the same time as the nonthermal heating rate. t p has no dependence on peak fluxes of Lyα or SXR. The flare types are independent on flare magnitude, duration, or location. The flare counts in each type match the sunspot counts over solar cycle. Interpretation The Lyα emission is likely contributed by some common processes in flares. By contrast, the gradual-phase peak (i.e., main peak for type III flares and sub-peak for type I flares) of Lyα that has a time delay relative to the SXR peak is supposed to be mainly caused by the thermal plasma cooling from flare loops. During a flare, the chromospheric plasma is heated and fills in the flare loops. Thus we see prominent hot emission in the SXR and also EUV wavebands. As the flare evolves, the hot plasma in flare loops suffers from a conductive or radiative cooling that usually takes place in the gradual phase. When the thermal plasma cools down to the formation temperature of the Lyα line, we could see evident emission enhancement in the Lyα waveband. This scenario is supported by our observational and modeling results, namely the delayed time of the gradualphase peak of Lyα (t ′ p ) or the plasma cooling time (t cool ) is related to the flare heating magnitude (represented by the impulsive-phase peak flux of Lyα) as well as the flare loop length. In fact, cool flare loops with a chromospheric temperature have been observed in recent years (e.g., Heinzel et al. 2018;Koza et al. 2019;Heinzel et al. 2020). In addition, Milligan et al. (2020) reported that the Lyα energy is comparable to, or about an order of magnitude smaller than the total thermal energy. This implies that the energy from thermal plasma is enough to provide the emission that is radiated by Lyα.
Note that in some of the flares, the Lyα emission can show multiple sub-peaks in the impulsive phase or in the gradual phase in addition to a main peak. The former case suggests that there probably exist multiple heatings during the flare, which likely take place in multiple loop strands (e.g., Hock 2012; Li et al. 2014a;Reep et al. 2018Reep et al. , 2019Reep et al. , 2020. While the latter case indicates that such kind of flares are probably EUV late-phase flares characteristic of different sets of loops with different lengths (e.g., Woods et al. 2011;Hock 2012;Sun et al. 2013;Li et al. 2014b;Woods 2014). Considering that a real flare consists of multiple loop strands (also implied from the error of the loop length measured from AIA images in the present study), multiple heatings/coolings are supposed to exist in a series of loop strands and cause a relatively long rise/decay or multiple peaks in the Lyα light curve.
It should also be noted that there are about one third flares whose Lyα emission does not show a well isolated impulsive-phase peak (i.e., sub-peak for type III flares) or a well identified gradual-phase peak (i.e., sub-peak for type I flares) but only an evident main peak (like the one shown in the top panel of Figure 2). One possibility is that these Lyα sub-peaks are too weak to be identified by eye or to be submerged in the background emission. Another possibility is that these Lyα sub-peaks overlap with the main peak and can not be well separated from the main peak, probably due to a long time heating (say, proceeding into the decay phase) or a very short time cooling (say, with a very short loop length).
Overall, the Lyα emission can be originated from both of a nonthermal electron beam heating and a thermal plasma cooling. These two processes are supposed to play a different role in contributing the Lyα emission in different types of flares. We conjecture that the nonthermal electron beam heating may play a major role in contributing the Lyα emission in type I flares that have a relatively larger impulsive-phase peak of Lyα. While the thermal plasma cooling may be more important in type III flares whose gradual-phase peak of Lyα is greater.

CONCLUSION AND DISCUSSIONS
In this paper, we have performed a statistical study on the relationship of the Lyα emission with the 1-8Å SXR emission in solar flares. It is worth highlighting that there are about one fifth (type III) flares whose Lyα emission peaks later than the SXR emission. This result is complementary to the prior studies on flaring Lyα that usually showed an earlier maximum in Lyα than in SXR (i.e., type I flares). Based on the Neupert effect check as well as the delayed time analysis and modeling, we conclude that the Lyα emission in different types of flares could be of nonthermal origin as well as thermal origin.
The nonthermal origin of Lyα has been illustrated in previous studies (e.g., Nusinov et al. 2006;Rubio da Costa et al. 2009;Milligan & Chamberlin 2016), though only with case studies for a few large flares. Our statistical results in the present work provide a further confirmation for that. Note that here we just use the time derivative of the SXR flux rather than HXR emissions to serve as the evidence for nonthermal electron beams due to unavailable HXR observations for numerous flare examples. Regarding the thermal origin of Lyα, we owe it to a delay of the Lyα peak with a few to tens of minutes relative to the SXR peak, which is further confirmed from the relationship between the delayed (or cooling) time and the flare loop length. In fact, it is natural and reasonable that the Lyα emission can come from the coronal loop cooling (e.g., Milligan et al. 2020). This can be predicted in the well-established standard flare model based upon decades of flare analyses and modelings.
Here we would like to draw attention that the thermal origin of Lyα can be of importance in some cases (say, in type III flares) and that the role of nonthemal and thermal processes may change from flare to flare, which could yield different types of Lyα emission curves. Note that we cannot rule out a thermal (or direct) heating that could also be contributed to the Lyα emission, as the radiative hydrodynamic simulations showed by Hong et al. (2019). Considering that a statistical study cannot touch upon all of the features of the Lyα emission, we propose that some more detailed case studies are needed to further explore the physical property of the flaring Lyα particularly with comprehensive observations from multiple instruments. For instance, the HXR observations from RHESSI can provide more clear evidence for the flare heating, i.e., nonthermal electron beams or not (e.g., Holman et al. 2011). The spectral observations at multiple temperatures from SDO/EVE can be used to better trace the cooling pattern of thermal plasma (e.g., Thiemann et al. 2017). Moreover, the imaging observations from SDO/AIA can give us the spatial information on both the flare loops and loop footpoints.
The present study can improve our understanding on the physical nature of the Lyα emission in solar flares, so as to help interpret the near-future Lyα observations made by the Solar Orbiter (Müller et al. 2019) and the Advanced Space-based Solar Observatory (ASO-S, to be launched around 2022; Gan et al. 2019). Solar Orbiter carries an Extreme Ultraviolet Imager (EUI; Rochus et al. 2019) that can observe the Sun closely at 0.28 AU in the Lyα waveband with a high spatial resolution. The Lyα Solar Telescope (LST; Li et al. 2019) on board ASO-S will obtain the full-disk images in Lyα from a Sun-synchronous orbit. These spatially-resolved observations contain much more detailed information on the Lyα emission features and the present study gives a preliminary basis for diagnosing the nonthermal and thermal origins of Lyα in solar flares. In addition, the present study can provide some constraints for the radiative hydrodynamic simulation of Lyα in solar flares, which is actually a powerful tool to interpret the optically-thick Lyα line emission. Note that the duration histogram has been cut to 120 min for a better display. Also note that the flare distance suffers a foreshortening effect and the increasing trend from disk center to the limb does not reflect a real distribution of flares per surface area. Bottom panels are the scatter plots of peak Lyα flux versus peak SXR flux, with the left one without a background subtraction and the right one with the background subtracted. In each of the panels, the red and blue colors mark the M-and X-class flares, respectively.   Figure 3. Histogram of the time difference (t p ) between the Lyα and SXR emission peaks (top) and scatter plots of t p versus peak Lyα and SXR fluxes (middle and bottom). Note that the t p histogram has been cut to a range of ±300 s for a better display. In the top panel, the orange, green, and blue colors denote the types I, II, and III flares, respectively, with the dark color marking X-class flares and the light color for M-class flares. In the middle and bottom panels, the red and blue dots represent the M-and X-class flares, respectively.    Figure 7. Histogram of the time difference (t d ) between the Lyα main peak and the peak of the time derivative of the SXR flux (top) and scatter plots of t d versus peak Lyα and SXR fluxes (middle and bottom) for all of the type I flares. Note that the t d histogram has been cut to a range of ±300 s for a better display. In each of the panels, the red and blue colors represent the M-and X-class flares, respectively.   Figure 9. The top panel shows the scatter plot of the delayed time (t ′ p ) of the gradual-phase peak of Lyα versus the impulsive-phase peak flux of Lyα for 173 type I flares (orange) and 62 type III flares (blue) in which both flare phase peaks can be well identified by eye. These two parameters have a very weak correlation with a coefficient of 0.1. The bottom panel shows the scatter plot of the delayed time (t ′ p ) of the gradual-phase peak of Lyα versus the flare loop length for 24 (out of 173) type I flares (orange) and 11 (out of 62) type III flares (blue) that are located within a circle of a radius of 500 ′′ on the solar disk. There appears to be a correlation between the two parameters, with a coefficient of 0.5.  Figure 10. EBTEL modelings for the 24 type I flares (orange) and 11 type III flares (blue) that are the same ones as shown in the bottom panel of Figure 9. The top left panel shows an example of the heating pulse (Q) with a Gaussian shape. The top right panel plots the time evolution of the average temperature (T ) and election density (n) of a single loop for an example. The bottom left panel gives the peak heating rate (Q p , diamond symbols) and the heating duration (σ, plus symbols) used in the modeling for the 35 types I (orange) and III (blue) flares. Here the flare loop length is measured from AIA images. The bottom right panel plots the relationship of the modeled cooling time (t cool ) with the loop length.