A Statistical Study of Soft X-Ray Flares on Solar-type Stars

The statistical characteristics of stellar flares at optical bands have received extensive study, but they remain to be studied at soft X-ray bands, in particular for solar-type stars. Here, we present a statistical study of soft X-ray flares on solar-type stars, which can help in understanding the multiwavelength behaviors of stellar flares. We mainly use Chandra Source Catalog Release 2.0, which includes a number of flaring stars with denoted variability, and Gaia Data Release 3, which includes the necessary information for classifying stars. We also develop a set of methods for identifying and classifying stellar soft X-ray flares and estimating their properties. A detailed statistical investigation of 129 flare samples on 103 nearby solar-type stars as selected yields the following main results. (1) The flare energy emitted at the soft X-ray band in our sample ranges from ∼1033 to ∼1037 erg, and the majority of them are superflares, with the most energetic one having energy of 6.0−4.7+3.2×1037erg . (2) The flare duration is related to its energy as formulated by Tduration,SXR∝Eflare,SXR0.201±0.024 , which is different from those derived at optical and near-IR bands, indicating distinct radiation mechanisms at different bands. (3) The frequency distribution of stellar flares as a function of energy is formulated as dNflare/dEflare,SXR∝Eflare,SXR−1.77 , which is similar to the results found at other bands and on other types of stars, indicating that the energy emitted at the soft X-ray band could be a constant fraction of the full-band bolometric energy.


INTRODUCTION
Stellar flares, which are intense energy-emitting phenomena taking place in the atmosphere of stars, are commonly argued to be caused by magnetic reconnection (e.g., Gershberg 2005;Benz & Güdel 2010;Shibata 2016).Except for long-term activities, which refer to the cyclic evolution of magnetic activities of stars and have likely an impact on the atmospheric temperature and chemical composition changes of planets (Chen et al. 2021), stellar flares, along with possible coronal mass ejections (CMEs), may cause disastrous space weather and affect the geomagnetic fields and thus the habitability of the planets orbiting these stars (e.g., Tsurutani et al. 2003).
Corresponding author: X. Cheng xincheng@nju.edu.cnWith new generations of instruments working at distinct bands, it becomes possible to conduct studies on stellar magnetic activities in detail, such as large-sample time-domain statistical surveys.Maehara et al. (2012Maehara et al. ( , 2015) ) and Shibayama et al. (2013) statistically studied stellar optical flares from solar-type stars with Kepler data.They found that flares on solar-type stars emit much more energy than solar flares, and the power-law index γ in the frequency distribution dN flare /dE flare ∝ E γ flare is roughly −1.8.With Transiting Exoplanet Survey Satellite (TESS) data, Doyle et al. (2020) and Tu et al. (2020Tu et al. ( , 2021) ) further confirmed the conclusions drawn with Kepler data.Pye et al. (2015) studied soft X-ray (SXR) flares from Tycho cool (F-M type) stars with XMM-Newton data and Getman & Feigelson (2021) studied SXR flares from pre-main sequence (PMS) stars with Chandra data.It is interesting that stellar SXR flares show a similar power-law index γ, i.e. about −1.8, to optical flares in the frequency distribu-Zhao et al.
tion.Such a law is also true for different types of stars, indicating that stellar flares may all share the same energy release process, i.e., magnetic reconnection.
Based on the magnetic reconnection theory, Maehara et al. (2015) derived that the bolometric duration of stellar flares is well associated with their energy, which can be formulated as T duration,bol ∝ E ∼1/3 flare,bol , in good agreement with solar white-light flares (e.g., Namekata et al. 2017).Such a correlation was subsequently verified by stellar optical and near-infrared (NIR) observations (e.g., Maehara et al. 2015;Tu et al. 2020Tu et al. , 2021)).However, it is unclear if such a power law is still valid and what is the slope if using flares observed on solar-type stars at the SXR band.Moreover, stellar flares often show different properties at different bands including the distributions of their duration and energy, thus it is also necessary to compare the derived statistical properties at the SXR band with that at the white-light band.
Flares could be accompanied with coronal mass ejections (CMEs) (e.g.Yashiro & Gopalswamy 2009;Aarnio et al. 2011), which are defined as large-scale expulsions of magnetized plasma from stellar atmosphere (e.g.Chen 2011; Veronig et al. 2021;Lu et al. 2022).Studies on solar flares revealed that such an association holds better for stronger and/or longer-duration (e.g.Cheng et al. 2010) flares, which indicates that the properties of stellar flares could be used to search for and diagnose potential CMEs.Type II radio bursts, i.e., slowly-drifting radio emission features, are also valuable for identifying and constraining the properties of stellar CMEs.By employing this approach, Crosley & Osten (2018) performed a multi-wavelength analysis of an M dwarf EQ Peg, and predicted the CME parameters according to solar scaling relationships.Very interestingly, Veronig et al. (2021) recently provided a strong indication for stellar CMEs based on the detection of coronal dimmings on stars.Loyd et al. (2022) also performed a coronal dimming analysis on a young K2 dwarf ϵ Eri and proposed that coronal dimmings can be used to constrain CME masses.
In order to thoroughly understand the nature of flares on solar-type stars, in this work, we conduct a detailed statistical study based on a combination of Chandra Source Catalog Release 2.0 (CSC 2.0) and Gaia Data Release 3 (Gaia DR3).The derived statistical properties are also expect to be able to help detect possible associated-coronal mass ejections.In Section 2, we first develop a set of methods for identifying and classifying stellar SXR flares and calculating their parameters.The results including flare properties and their correlations, as well as occurence frequency distributions, are shown in Section 3, in which we also make a comparison be-tween stellar and solar flares.A brief summary is given in Section 4.

Selection of flare candidates
We search for stellar flares using data in the Chandra Source Catalog Release 2.0 (CSC 2.0)1 collected by NASA's Chandra X-ray Observatory (CXO)2 through the end of 2014.CXO has an unprecedented spatial resolution of subarcsecond for most detected sources, which allows us to match their optical counterparts unambiguously.Our primary database, CSC 2.0, is the second major release of the definitive catalog of X-ray sources detected by the CXO.CSC 2.0 provides properties for 317, 167 individual X-ray sources and 928, 280 individual observation detections, allowing statistical analysis of large samples and also individual case studies.There are information about each source across 5 science energy bands (broad, hard, medium, soft, and ultra-soft) for Advanced CCD Imaging Spectrometer (ACIS) and 1 band (wide) for High Resolution Camera (HRC)3 .All of the 6 science energy bands and their effective energies are given in Figure 1.
In this work, we search for stellar flare candidates from the ACIS data by using a parameter called Variability Index (var_index) available in CSC 2.0.Variability Index is an integer in the range of [0, 10] which is calculated for each science energy band.The calculation method combines the Gregory-Loredo variability probability with the fractions of the multi-resolution light curve output by the Gregory-Loredo analysis that are within 3σ and 5σ of the average count rate4 .According to the statements on the websites of CSC 2.0, we evaluate whether the source region flux is uniform throughout the observation out of the value of this parameter.The events in CSC 2.0 with var_index ≥ 6 are commented as a definite variable, which is used for identifying possible flaring stars.
We first select the sources observed by ACIS with var_index ≥ 6 in at least one ACIS or HRC band, and obtained 35,212 events.We cross-match these events with Gaia Data Release 3 (Gaia DR3) 5 catalog, using a relatively small radius of 1 ′′ in order to restrain the proportion of false-matched results.The cross-matching Figure 1.Table 1 process is done by using the TOPCAT6 .When doing the cross-matching, the proper motion of the stars can be neglected because none of the sources have a total proper motion greater than the error of the angular distance between the two catalogs, which is calculated from the coordinates given by the two catalogs and their errors.Among the Gaia DR3 counterparts, there are 8,688 events with credible distances, i.e. the parallaxes are available and the relative parallax errors are lower than 0.2, given that the distance is a critical parameter when calculating the fluence or energy of a given stellar flare event.With parallaxes p and parallax errors σ p , one can calculate the distances d and distance errors σ d of the sources.
Next, we try to exclude astrometrically or photometrically suspectable sources, by using Gaia DR3 parameters astrometric_excess_noise, which characterises the goodness of fit of the astrometric solution to the observations, and phot_bp_rp_excess_factor, which is related to photometric reliability.We require astrometric_excess_noise to be lower than 1, and phot_bp_rp_excess_factor, which is color-dependent, to be low enough so that the over-densed region in the phot_bp_rp_excess_factor -bp_rp diagram is included and unusually high values are excluded.Such a sifting process follows the methods in the Section 4 of Tutorial: Exploring Gaia data with TOPCAT and STILTS7 (Taylor 2016).4,622 CSC events remain after this sifting process.These flare candidate events are related to 1,680 independent CXO exposures, which are identified with an independent observation identifier (ObsID) each.Besides, we download a Gaia DR3 cata-log which contains all of the sources within 100 pc, and process it by using the same method mentioned above to exclude astrometrically or photometrically spurious sources.With these numerous Gaia sources as a background, we can easily distinguish the positions of our CSC sources on a Hertzsprung-Russell diagram (HRD) and further specify their types.

Selection of solar-type stars
We used the effective temperature (T eff ) and the surface gravity in logarithmic scale (log g) available from Gaia DR3 to select solar-type stars.Following the criteria used in previous studies (e.g., Maehara et al. 2012;Shibayama et al. 2013;Maehara et al. 2015;Tu et al. 2020Tu et al. , 2021)), we set the selection criteria as follows: (1) 5100 K ≤ T eff < 6000 K and (2) 4.0 ≤ log g < 4.8.
Out of the 4,622 events, the number of events satisfying the selection criteria, i.e. flare candidates on solartype stars, is 528.

Analysis of light curves
In this work, we define target regions by encircled counts fraction (ECF).We set r 90 as the radius of each target region, which is defined by a circle around the target source inside which ECF reaches 90% and above.Besides, the background region is defined as an annulus with inner radius 2r 90 and outer radius 3r 90 .Photon count rates (denoted as CR in the following parts) in the two regions mentioned above are denoted as CR tar and CR bkg , respectively.
We then define the real CR(t) coming from the target source in each time bin as where CR bkg (t) is multiplied by a factor of 1/5 because the area of the background region is 5 times the target region.
We search for stellar flares by inspecting light curves of events from solar-type stars (see Section 2.2) among the 528 candidate events mentioned above.In this work, we plot the light curves, the CRs against the arrival times, after binning the arrival times with an appropriate time interval, which specifically means that there should be at most 200 bins in one light curve with no more than 25% of the bins containing no photons.In order to determine the times at which the photon count rate changes obviously and grasp the general tendencies of the light curves, we adopt an approach named Bayesian Blocks, which is widely used to identify and characterize flaring behavior in X-ray sources (e.g., Scargle 1998;Scargle et al. 2013;Getman & Feigelson 2021).After fitting the original light curve with Bayesian Blocks, we get a series of times when the CR has noticeable changes, i.e. changepoints, as well as a series of stepwise constant CR, as shown in Fig. 2.
Before identifying stellar flares in our data, it is needed to define the quiescent phase and flaring phase in a light curve.Here, we set the lowest value of CR of Bayesian Blocks during an event as the quiescent level (denoted as CR q ).We define a relative count rate ∆CR as (2) and then define the start time t start and end time t end of a flare as the first and last changepoints between which the ∆CR of Bayesian Blocks are consecutively higher than 5 times the photometric errors, i.e.
where σ CR denotes the photometric error with the same unit as CR.The numerical value of σ CR is estimated by the root-mean-square error of CR q , i.e.
where ∆t q is the corresponding time range of CR q .Then, the quiescent phase is defined as the time interval (t 0 , t start ) ∪ (t end , t 1 ) for events in categories F or S, or (t end,i , t start,i+1 ) as the quiescent phase between the ith and (i + 1)th flaring phases for events in categories M-F or M-S; while the flaring phase is defined as the time interval (t start , t end ) for events in categories F or S, or (t start,i , t end,i ) as the ith flaring phase for events in categories M-F or M-S, where t 0 and t 1 refer to the start and end time of the corresponding CXO exposure respectively.
All of the light curves and Bayesian Blocks fitted to them are examined both numerically and then visually, and classified into eight categories as illustrated in Fig. 2: (1) Constant (C): This category contains events whose light curves have no Bayesian Blocks with ∆CR higher than 5σ CR .This commonly occurs when an event is variable but the variability of this event is not high enough, making it difficult to determine whether this event contains a typical stellar flare.We exclude these events in the next steps of this work.
(2) Flare (F): This category contains events whose light curves show a typical shape of stellar flare, i.e. with a quick rise and a gradual decay.Meanwhile, we require the events in this category to present at least two different levels of CR, or at least one changepoint within the flaring phase.These events are most useful for our analyses, and their time parameters, peak luminosities and total energies are readily calculated.
(3) Flare with single Bayesian Block (S): This category contains events whose light curves show a clear rise and decay, but with only two major changepoints during the full exposure, i.e. there is only one single CR level of Bayesian Block within the flaring phase.The Xray durations, peak luminosities and total energies are calculated accurately for these cases.However, the rise time and decay time cannot be determined accurately and set to be equal according to our definitions (see Section 2.4).(4) Rise (R): These are events whose light curves display a quick rise, but do not fall back to the quiescent level before the exposure ends.That is to say, there is no changepoints after which ∆CR satisfies ∆CR ≤ 5σ CR consecutively.Because of the incomplete time coverage, these cases can only show the lower limits of their durations, peak luminosities and total energies.The decay times are not available in these cases.The rise times are measurable only when the peak is covered in the light curve.Otherwise, only the lower limits to the rise times are available.
(5) Decay (D): These are events which start from a level higher than the quiescent one but show a gradual decay.That is to say, there is no changepoints before which ∆CR satisfies ∆CR ≤ 5σ CR consecutively.Similarly to category R, these cases only show the lower limits of their durations, peak luminosities and total energies.The rise times are not available in these cases.The decay times are measurable only when the peak is covered in the light curve.Otherwise, only the lower limits to the decay times are available.(6) Flare with multiple flaring phases (M): This category contains events whose light curves show two or more flaring phases.More specifically, each flaring phase has count levels with ∆CR > 5σ CR ; however, between two consecutive flaring phases, there exists at least one count level with ∆CR ≤ 5σ CR .In these cases, the CR reaches an earlier peak, falling back to the quiescent level, and it then goes up again.Such a process may be repeated several times.Therefore, we regard them as showing multiple flares within one CXO exposure.These cases are classified as multiple flare events, each of which can be classified into F, S, R or D categories according to the characteristics of the flaring phase.In the following statistics, we selectively treat each flaring phase as a single flare event depending on the parameters we analyse.(7) Variable (V): Some events show strong variability (with some Bayesian Blocks with ∆CR > 5σ CR while others with ∆CR ≤ 5σ CR ), but without a typical shape of quick rise and gradual decay.In a few cases, there are even quasi-periodic variations.These phenomena are possibly caused by transits between binaries or in multiple-body systems.Apparently, it is quite difficult for us to search for stellar flare events among these cases, because even though a star flares in such systems, the flare signals are mostly obscured by a non-flare variation.Thus, we exclude these events.(8) Bad (B): There are mainly two reasons causing bad qualities of the light curves.One possible reason is that some sources are located near or at the edge of one single CCD.With the vibration of the instrument etc., not all of the photons coming from the source could be received, causing a dither effect8 .Another possible reason is that the source is too faint to emit enough photons to the detector, as a result of distance or the effect of interstellar medium etc.In the above cases, the signal cannot be clearly delineated in the light curve.Events in this category could not be used for our analyses.
Among the eight categories as elucidated above, the categories C, F, R, D and V follow the criteria adopted by Getman & Feigelson (2021) with some modifications, while the other three categories S, M and B are newly defined for the sake of covering all types of light curves retrieved from our data.
Here, we use the events in categories F and S (including F and S events extracted from category M, denoted as M-F and M-S respectively) to measure the durations, peak luminosities and total energies, and the events in categories F and M-F to measure the rise and decay times.Among the 528 flare candidates on solar-type stars, 30 events belong to category F or M-F and 99 events belong to category S or M-S.The total number of individual sources in our sample set is 103.The numbers of events (light curves) in each category are listed in Figure 3.
The 103 flare host stars are cross-matched with the SIMBAD catalog9 to check if there are any confirmed binary stars, and none of these sources are flagged as binary systems.However, there are still potential unresolved binary stars in our sample set, which is limited primarily by the photometric and spectrometric accuracies nowadays.The renormalised unit weight error (RUWE) from Gaia DR3 can serve as an indicator, since unresolved binary stars tend to show large RUWE values, e.g.> 1.4 (e.g., Stassun & Torres 2021;Chulkov & Malkov 2022).The RUWE values of the flare-host stars are listed as a column of Figure 5. Besides, those distribute above the mail-sequence band in Fig. 4c are especially likely to be binary systems, since they show higher luminosities than most main-sequence stars with the same temperatures.
The equatorial positions of the samples are illustrated in Fig. 4a; the probability distribution of the distances of the flare host stars is shown in Fig. 4b; and the samples in an HRD are plotted in Fig. 4c.

Flare properties and error analysis
The start time t start and end time t end are defined in Section 2.3.The peak time t peak of a flare is defined as the median time of the Bayesian Block with the highest CR during the flaring phase.Other useful time parameters of stellar flares, including the rise time T rise , decay time T decay and total duration T duration , are defined as follows: Time parameters t start , t peak and t end are illustrated in Fig. 2 with dash-dotted lines in green, red and blue, respectively.
In order to obtain flare parameters such as spectral indices, peak luminosities and total energies etc., we retrieve the spectra of the total 129 flare samples within each target region as mentioned in Section 2.3.The spectra are extracted with the command specextract in the software package CIAO10 by the Chandra X-Ray Center, and from 5 different types of phases: quiescent phase, flaring phase, rise phase, decay phase and    peak phase.The definitions of quiescent phase and flaring phase are already mentioned in Section 2.3.The rise phase is defined as the time interval (t start , t peak ) for events in categories F or S, or (t start,i , t peak,i ) as the ith rise phase for events in categories M-F or M-S.Similarly, the decay phase is defined as the time interval (t peak , t end ) for events in categories F or S, or (t peak,i , t end,i ) as the ith decay phase for events in categories M-F or M-S.The peak phase is defined as the time interval of the Bayesian Block with the highest CR during each flaring phase, which could be denoted as (t peak,start , t peak,end ) for events in categories F or S, where t peak,start and t peak,end are the start and end time of this Bayesian Block respectively, or (t peak,start,i , t peak,end,i ) as the ith peak phase for events in categories M-F or M-S.We set the spectrum for the quiescent phase as the background of the other spectra with the command fparkey in CIAO, thus the quiescent radiation of the star and the background radiation are deducted, leaving the pure radiation of the flare.All of the retrieved spectra are then adaptively binned and fitted.The fitting procedure is performed with NASA's X-ray spectral fitting package Xspec11 , and the model we use, i.e. the absorbed power-law spectra with pile-up, could be expressed as pileup*TBabs*powerlaw in Xspec.The pileup and TBabs components are convolved into the model to consider the pile-up effect (e.g., Davis 2001) and absorption from the interstellar medium respectively.An unabsorbed power-law spectrum is formulated as where dN photon dt dE is in the units of counts s −1 erg −1 , K is the normalization coefficient and α is the dimensionless spectral index.After doing so, we obtain the spectral indices in the flare rise and decay phases, which are denoted as α rise and α decay , respectively.
We obtain corrected fluxes within 0.5 − 8.0 keV for the flaring phase (denoted as F flare ) and peak phase (denoted as F peak ) by using the flux command to the powerlaw component.And their errors are calculated from the normalization coefficients K, i.e.
where σ ± F denotes the positive and negative errors of fluxes in linear scale.It's worthy of noticing that both F and σ ± F are corrected by a factor of 1 0.9 because the target region, i.e. the circle with the radius of r 90 , only encircles 90% of the total flux of the source (see Section 2.3).Then the average luminosities for the flaring phase or peak phase could be obtained with assuming that the radiation from the star during the flaring phase is isotropic.The average luminosity L peak during the peak phase is a good estimation of the real peak luminosity of a flare event when t peak,end −t peak,start is relatively small.In such a case, L peak is regarded as a proxy of the real peak luminosity.The errors of the luminosities are expressed as according to the error propagation formula, where d and σ d are derived in Section 2.1.Equations ( 10) and ( 11) are calculated for both the whole flaring phase and peak phase.The average luminosities and their errors during the flaring phase are denoted as L flare and σ ± L flare respectively, while the average luminosities and their errors during the peak phase are denoted as L peak and σ ± L peak , respectively.
With the average luminosity during the flaring phase L flare , we then estimate the total energy released by the flare with and the error of E flare is Note that, σ ± L peak and σ ± E flare are estimated by a combination of the error in distance and that arising from the powerlaw model in Xspec while fitting the spectra.

Estimation of flare frequencies
In order to obtain the frequency distribution of flares, we inspect a bigger data set which contains all of the solar-type stars both with and without detected flare events, and use the same rules to select our sample set as listed in Figure 3.More specifically, we identify such a data set with the same criteria mentioned in Sections 2.1 and 2.2, but without checking Variability Indices and light curves, because all of the other selection rules are unrelated to stellar activities.The total number of individual detections (per source per CXO exposure) in this data set is 8,016, with 2,748 individual sources and 2,653 independent CXO exposures.We denote the effective exposure time of each individual detection as τ n , which is obtained from the total elapsed time of the valid data of each observation, gti_elapse, available in CSC 2.0 data12 .
We use the following equation to estimate the flare frequency distribution as a function of flare energy on solar-type stars: where ∆N flare is the number of detected flares in each flare energy bin, specifically meaning the number of events in light curve categories F, S, R, D, M-F, M-S, M-R and M-D.Considering that it is impossible to accurately calculate the flare energies of flares in categories R, D, M-R and M-D, we omit these four categories, which means that the flare frequency f (E flare ) in equation ( 14) is only a lower limit.In equation ( 14), ∆E flare represents the energy range of each flare energy bin, and τ is the total effective exposure time of all of the 8,016 individual CXO detections and is calculated by We can see that such a flare frequency distribution function f (E flare ) is in the units of flares star −1 yr −1 erg −1 .After defining the frequency distribution function, we further define the corresponding normalised cumulative occurrence probability as a function of flare energy: Similarly, this equation only gives a lower limit to the actual cumulative distribution, since only F, S, M-F and M-S samples have been taken into account.

RESULTS AND DISCUSSIONS
After determining the physical quantities about stellar flares and corresponding host stars by using the methods in Section 2, we present the statistical results in this Section.Physical quantities in all figures are denoted with a subscript SXR, in order to specify that these quantities refer to the CXO soft X-ray band.
The probability density functions of some quantities in logarithmic scale, such as time parameters, peak luminosities, flare energies and spectral indices, can be well fitted by normal distributions.Note that, we apply normal distributions here only to statistically estimate the corresponding mean values and standard deviations of the key parameters without specific physical assumptions.
Physical and observational properties of flare host stars are listed in Figure 5, and the properties of corresponding stellar flare events are listed in Figure 6.Only 10 events with the highest flare energies are shown in the tables as examples, while the entire tables are available online in a machine-readable form.

Distributions of time parameters
We obtain the time parameters of flares in categories F, M-F, S and M-S.The flare rise time T rise,SXR and the decay time T decay,SXR are only measured for 30 F and M-F samples, while the flare duration T duration,SXR is measured for all of the 129 F, M-F, S and M-S samples.
Distributions of time parameters T rise,SXR , T decay,SXR and T duration,SXR of F and M-F samples are shown in Fig. 7.The flare rise time T rise,SXR spans about two orders of magnitude, ranging from ∼ 10 2.5 s to ∼ 10 4.5 s and centering at 10 3.34±0.49(1σ)s.The flare decay time T decay,SXR spans about one order of magnitude, ranging from ∼ 10 3.5 s to ∼ 10 4.5 s and centering at 10 4.04±0.40(1σ)s.The flare duration T duration,SXR spans about one order of magnitude, ranging from ∼ 10 3.5 s to ∼ 10 4.5 s and centering at 10 4.11±0.37(1σ)s.It is found that the decay phase occupies most of the entire flare duration in most of events.
Distribution of the ratio of decay time to rise time T decay,SXR /T rise,SXR for F and M-F samples is shown in Fig. 8.The ratio of decay time to rise time ranges from ∼ 1 to ∼ 10 1.5 with an average of 10 0.62±0. 32(1σ) .The majority of flares have T decay,SXR > T rise,SXR , except for only one event of F or M-F category showing T decay,SXR /T rise,SXR = 0.668.However, we suspect that this special event is a result of two or more overlapping events that are not strictly identified by our identification methods (see Section 2.3), as shown in Fig. 18 in Appendix A. We thus conclude that the SXR profile of stellar flares consists of a rise phase and a much longer decay phase, very similar to the characteristic of solar SXR flares, even though the stellar flares show a much longer duration than the solar flares, which typically last for ∼ 10 2 to ∼ 10 4 s with a median duration of ∼ 400 s (Veronig et al. 2002).For solar flares, it was found that the association between flares and CMEs increases with the class and/or duration of flares (e.g.Cheng et al. 2010).Thus the flares in F or M-F category, in particular for the longest-duration ones, are most likely accompanied by stellar CMEs.
The distribution of flare duration for all F, M-F, S and M-S samples is displayed in Fig. 9.The flare duration is found to range from ∼ 10 2.0 s to ∼ 10 4.5 s with a median of 10 3.76±0.44(1σ)s.It is obvious that, after including S and M-S samples, the distribution of T duration,SXR moves towards the shorter time scale.It indicates that S and M-S flares last for a shorter time than F and M-F flares.The difference between the two groups is mainly due to two reasons: (1) the events with a shorter duration tend to be described with fewer blocks when fitted with the Bayesian Block method, or (2) the duration of S and M-S events can be underestimated when the corresponding changepoints are regarded as t start and t end as mentioned in Sections 2.3 and 2.4.

Distributions of peak luminosities and flare energies
We obtain the flare peak luminosities and flare energies of 129 flare events in categories F, M-F, S and M-S.Distributions of flare peak luminosities L peak,SXR and flare energies E flare,SXR are shown in Fig. 10a and 10b, respectively.The flare peak luminosity ranges from ∼ 10 29 erg s −1 to ∼ 10 34 erg s −1 , centering at 10 31.07±0.68(1σ)erg s −1 .The flare energy varies from ∼ 10 33 erg to ∼ 10 37 erg with an average of 10 34.86±0.71(1σ)erg.
The energy of the strongest flare event among our samples is 6.0 +3.2 −4.7 × 10 37 erg.It is detected in ObsID = 13657 and has T rise = 3593 s and T decay = 30061 s.The corresponding host star is 0.4559 ± 0.0044 kpc from us and has a surface effective temperature of 5580 K and surface gravity of log g = 4.444.This host star is classified as a pulsating variable star according to Hartman et al. (2008), after cross-matching with SIMBAD.
One of the most energetic solar flare event in record is known as the Carrington event, which occurred on September 1st 1859 (Carrington 1859) and released a total bolometric energy of E flare,bol ∼ 5 × 10 32 erg (Cliver & Dietrich 2013).By comparison, the weakest flare from solar-type stars in our sample released an energy of E flare,SXR = 0.9 +5.0 −0.4 × 10 33 erg.It is worth mentioning that the flare energy emitted at the SXR band is only a fraction (∼ 1/10 − 1/4) of that at the optical band for events with E flare,SXR between ∼ 10 34 and ∼ 10 36 erg (Flaccomio et al. 2018) on PMS stars.Kretzschmar (2011) also suggested that, for solar flares, the energy at the optical band represents about two thirds of the bolometric energy.Thus, if the flares on stars originate from the similar physical mechanism to those on the Sun and show a similar energy fraction at the SXR band, the energy of flares at the SXR band E flare,SXR may constitute only a minor fraction of the total bolometric energy E flare,bol .Specifically, the calculated E flare,SXR here is roughly ∼ 1/10 of the bolometric energies E flare,bol .3 By comparison, the maximum bolometric energies in previous studies are listed in column (5) of Figure 11.We can see that in most cases, the maximum flare bolometric energy is 10 36 − 10 38 erg under the assumption that flares on different stars show a similar energy fraction at the SXR band, indicating that it is saturated at ∼ 10 38 −10 39 erg above which nearly no flares have been detected.The upper limit of flare bolometric energy derived by Wu et al. (2015) is 2×10 37 erg, which, however, might be underestimated as a result of selection effect of instruments or different methods of calculating flare energy.In this work, the flare energy is calculated by fitting the spectrum with an absorbed power-law spectrum and then calculating the average unabsorbed flux during the entire flaring phase (see Section 2.4), while Wu et al. (2015) calculated the flare energy at the optical or NIR bands by assuming white-light flares as blackbody emission (e.g., Maehara et al. 2012Maehara et al. , 2015;;Shibayama et al. 2013).
Moreover, solar flares are usually classified in terms of their Geostationary Operational Environmental Satellite (GOES) 1−8 Å peak SXR flux F peak,SXR .Figure 12 lists the flare classes and the corresponding GOES peak SXR flux (Cliver & Dietrich 2013), as well as the derived flare bolometric energy E flare,bol (Notsu et al. 2019;Doyle et al. 2020).The Carrington event was classified as an X45 (±5) class flare according to Cliver & Dietrich (2013).Since both GOES and CXO work at SXR bands, with the derived F peak,SXR at 1 AU, we estimate that all of the flare events in this work should be X-class ones, roughly ranging from X100 to X10, 000, 000, which are significantly greater than solar flares.Based on the fact that the Carrington event caused an enormous geomagnetic storm (Tsurutani et al. 2003), thus threatening communication or navigation systems (Maehara et al. 2015), it is reasonable to infer that the megaflares, at least superflares, as found here could have serious influences on the space weather and atmosphere of associated exoplanets if they exist.

Correlations among T duration , L peak and E flare
We further study the correlations between the time and energetic parameters of the flares, which are also fitted with the following power-law formulae: T duration,SXR ∝ E β flare,SXR . (18) Figure 6.Table 4 The scatter plot of the peak luminosity against the flare energy is shown in Fig. 13a, and the best-fitting with a power-law relation is illustrated.We use σ δ to denote the standard deviation of δ, as obtained via Monte-Carlo simulations.Only the values whose lower limits are greater than zero are taken into account when doing the Monte-Carlo simulations.The power-law index is δ = 0.76 ± 0.11(1σ δ ).The Pearson correlation coefficient is 0.77 (with a p-value of 5.6 × 10 −23 ), indicating that there is a strong correlation between L peak,SXR and E flare,SXR .
From the power-law correlation of equation ( 17  A similar scaling index of 1.16 was also reported by Wolk et al. (2005) from their Chandra survey of 41 Ktype PMS stars, although they used the average luminosity during the entire flaring phase (L flare,SXR ) rather than the peak luminosity as we use here.The similar-  ity in the scaling index between E flare,SXR and L peak,SXR and that between E flare,SXR and L flare,SXR indicates that L flare,SXR is nearly proportional to L peak,SXR .Since 1/δ is close to 1, we can say that the flare energy varies nearly linearly with the peak luminosity, which implies that the flare duration may be independent of its peak luminosity, as argued by Wolk et al. (2005) and Pye et al. (2015).To test this point, we make a scatter plot of T duration against L peak,SXR as shown in Fig. 14, and obtain a Pearson correlation coefficient of −0.095 (with a p-value of 0.32) between them, which confirms the independence between the two parameters.
The scatter plot of the duration as a function of flare energy is shown in Fig. 13b.The power-law index is β = 0.201±0.024(1σβ ), where σ β stands for the standard deviation of β, similar to σ δ as mentioned above.The Pearson correlation coefficient is 0.38 (with a p-value of Zhao et al. 2.8 × 10 −5 ), indicating that there is a weak correlation between T duration,SXR and E flare,SXR .
For comparison, the best-fitting index β of T duration as a function of E flare in logarithmic scale from previous studies are listed in column (6) of Figure 11.According to the magnetic reconnection theory of stellar flares, the flare energy originates from the release of magnetic energy and the duration is comparable to the reconnection time (e.g., Shibata 2016).Assuming that all solar-type stars have similar stellar magnetic properties (e.g., magnetic field strength B) and applying dimensional analysis, Maehara et al. (2015) derived a power-law index of β ∼ 1/3 in the relationship T duration,bol ∝ E β flare,bol , which is in good agreement with the results for solar flares (e.g., Namekata et al. 2017).It is found that the indices derived by using the data at optical and NIR bands are usually consistent with the theoretical value of β ∼ 1/3, such as 0.39 ± 0.03 obtained by Maehara et al. (2015) and 0.42 ± 0.01 by Tu et al. (2020Tu et al. ( , 2021)).
It is obvious that the value of β (0.201 ± 0.024(1σ β )) derived here is much smaller than previous values.We calculate the Bayesian Information Criteria (BIC) values of our fit and the fit where the slope is fixed to 1/3 and find that they are −211.46and −208.02,respectively, with a difference of 3.44.This reconfirms that the value of β we derived is indeed significantly different from β = 1/3.Unlike optical (white-light) flares, which well represent the overall bolometric feature of flare events, the duration of the SXR flares increases more slowly with the energy than the bolometric ones.Such a significant difference indicates that the SXR emissions of stellar flares originate from different radiation mechanisms from the optical ones, as also reflected from the distinct energy partitions at the SXR and optical bands as derived in Section 3.2.

Distributions of spectral indices
We obtain the power-law spectral indices for the SXR spectra based on the sample of 129 flares in categories F, M-F S and M-S.The rise and decay phases of the flares are analysed separately.Note that we exclude those events with non-positive power-law indices for further statistics.
After calculating the ratio of the spectral index during the decay phase to that during the rise phase and fitting the probability distribution as a function of the index ratio (in logarithmic scale), as shown in Fig. 15c, we obtain the most probable value of log 10 (α decay,SXR /α rise,SXR ) being 0.08, which is a little bit higher than zero.Some of the flare events have α decay,SXR significantly higher than α rise,SXR , implying that the spectrum for the decay phase of these events tends to be softer than that for the rise phase.It implies that for these events, more non-thermal emissions could be generated during the rise phase.

Occurrence frequency distribution of stellar flares
We use the methods mentioned in Section 2.5 to estimate flare occurrence frequencies.The results given by equations ( 14) and ( 16) are presented in the two panels of Fig. 16, respectively.
The frequency distribution as a function of flare energy, i.e. f SXR (E flare,SXR ), is fitted by a power-law function for a specific SXR energy range, as shown in Fig. 16a.The fitting gives dN flare /dE flare,SXR ∝ E −1.77 flare,SXR , where we select the energy range from 1.68 × 10 35 erg to 1.76 × 10 36 erg beyond which we lack sufficient flares  5 to ensure a reliable fitting (Maehara et al. 2012).The fitting result suggests that the average occurrence frequency of stellar flares on each solar-type star is at least once in about 2.89 years in the SXR energy range of 10 34 − 10 35 erg, once in about 16.8 years in the SXR energy range of 10 35 − 10 36 erg, once in about 98.2 years in the SXR energy range of 10 36 − 10 37 erg, and once in about 573 years in the SXR energy range of 10 37 − 10 38 erg.
The flare frequency obtained in this work is higher than those at the optical band in previous works (e.g.Maehara et al. 2012Maehara et al. , 2015;;Shibayama et al. 2013), and such a difference is mainly due to the following reason.As a stellar flare occurs, the radiation at the SXR band often increases significantly in magnitude, while that at the optical band only increases slightly, which thus makes stellar SXR flares to be detected more easily than the optical ones.In terms of physical mechanisms of solar flares, generally speaking, it is relatively difficult for the energy released in the corona to be transported to the photosphere (e.g., Benz 2017;Song & Tian 2018).
Zhao et al.It is interesting that the power-law indices we derived above are quite similar to what have been obtained from different types of stars at different bands (all roughly −1.8; e.g.Shibayama et al. 2013).More interestingly, the spectral indices for stellar flares also resemble with the values for solar activities, i.e., −1.79 for solar extreme ultraviolet (EUV) nanoflares (Aschwanden et al. 2000), −1.74 for solar SXR microflares (Shimizu 1995), and −1.53 for solar hard X-ray (HXR) flares (Crosby et al. 1993).Such a similarity implies that, even though solar and stellar flares are from different types of stars, they may share the same physical mechanism, most likely magnetic reconnection (e.g., Maehara et al. 2015;Shibata 2016;Namekata et al. 2017).Thus, the frequency distribution of stellar flares can also serve as a prediction of the probability and frequency of potential superflares occurring on our Sun.Meanwhile, the similarity of γ obtained at the SXR band and the bolometric value derived at the optical band is a good indication that the energy emitted at the SXR band is nearly a constant fraction of full-band bolometric energy, independent of the magnitude of flares.

Completeness
The flares with lower flux or shorter duration are usually more difficult to be detected, thus resulting in a selection bias or incomplete sample set.In order to check the completeness of our sample, we perform a simulation and generate a series of simulated flare light curves with different T duration and CR peak (see Appendix B).We then analyse these simulated light curves using the same methods as mentions in Sections 2.3 and 2.4, and obtain their light curve categories and time parameters.If a light curve is classified into F, S, R, D or M categories, we regard it as a "detected event", otherwise (C,  where the number of simulated events at a certain (T duration , CR peak ) is fixed to be 100 in this simulation.
The completeness distribution is shown in Figure 17.At each certain CR peak , the completeness first rises with T duration and then falls.The completeness is relatively low when T duration is very long because the exposure time could hardly cover the full flare event, thus causing CR q overestimated and then hard to identify a flaring phase (see Section 2.3).For most of the flare events we study, the completeness is higher than 50% even though only a few events have the completeness higher than 90% according to the simulation.The histograms of the time and energetic parameters in this work could be affected by the completeness issue, as the low-flux or short-duration flares often have a low completeness, or the probability of detection is low.If more low-flux or short-duration flares are detected, the distributions of time and energetic parameters could be more like power- law distributions, like flares on the Sun (e.g., Kou et al. 2020).Here, we tend to only present the observational data in the histograms because some assumptions that are made in the completeness simulation (see Appendix B) have not been verified.The simulation results can be used as a guide for determining whether the observational sample is complete but are not intended for further analysis.However, it should be addressed that the major conclusions drawn from the flares with the highest peak flux and the longest duration as derived in this work, as well as the correlations among the time, energetic parameters and the occurrence frequency distribution, are not influenced by the completeness issue.

SUMMARY
We have made a statistical study for 129 flare events on 103 nearby solar-type stars from Chandra's CSC 2.0 catalog.We presented a comprehensive set of methods for identifying and classifying SXR stellar flares, and obtained the following main results and conclusions: 1.The flare duration at the SXR band ranges from ∼ 10 2.0 to ∼ 10 4.5 s, with the majority of flares showing a longer duration than typical solar SXR flares.The duration of the decay phase occupies most of the entire duration of flares.
2. The flare peak luminosity at the SXR band ranges from ∼ 10 29 to ∼ 10 34 erg s −1 and the total flare energy at the SXR band emitted varies from ∼ 10 33 to ∼ 10 37 erg, which is much larger than the strongest solar flare ever recorded.The energy of the strongest stellar flare is found to be 6.0 +3.2 −4.7 × 10 37 erg.The majority of flares in our sample are superflares, 3 of which are megaflares (with E flare,SXR exceeding 10 36.2 erg).An upper limit seems to exist in the bolometric energy for a single event, which is ∼ 10 38 − 10 39 erg.
3. The duration of the flares is essentially independent of the peak luminosity.The duration as a function of energy at the SXR band can be best fitted with T duration,SXR ∝ E 0.201±0.024flare,SXR , where the power-law index is apparently lower than those derived at optical and NIR bands, indicating the distinct radiation processes working at the different wave bands.
4. The spectra of some of the flare rise phases are harder than that of the corresponding decay phases, to say, the rise phases of these events could be more non-thermal.
5. The frequency distribution of stellar SXR flares as a function of energy can be well fitted by dN flare /dE flare,SXR ∝ E −1.77 flare,SXR in the range from 1.68 × 10 35 erg to 1.76 × 10 36 erg.The fitted spectral index is similar to the results derived for flares from other types of stars or the Sun.We perform a simulation of light curves in order to check the completeness of our sample.The simulation methods are mainly adopted from the Appendix D of Mossoux & Grosso (2017).We consider a Gaussian-shaped flare light curve superimposed on a constant level (i.e.CR q,intrinsic ).The input, or "intrinsic" light curve (i.e.CR intrinsic (t)) can be expressed as CR intrinsic (t) = CR q,intrinsic + CR peak,intrinsic exp − 1 2 where CR peak,intrinsic is the maximum count rate of the Gaussian component and σ T stands for the standard deviation of the corresponding Gaussian distribution.σ T and CR peak,intrinsic are the input values of a simulated light curve.Other properties of a simulated light curve include the total "exposure" time T exposure , the quiescent count rate CR q,intrinsic and the peak time t peak,intrinsic .When performing the simulation, we let log 10 T exposure and log 10 CR q,intrinsic obey the normal distribution in logarithmic scale fitted by the observational data, i.e. log 10 T exposure ∼ N (µ = 4.391, σ = 0.496), (B2) log 10 CR q,intrinsic ∼ N (µ = −2.413,σ = 0.587), (B3) and let t peak,intrinisc randomly distribute between 0 and T exposure .At each (σ T , CR peak,intrinsic ), we simulate 100 events and calculate the completeness by using equation (20).Figure 19 shows an example of the simulated light curves.It's necessary to obtain the correlations between (σ T , CR peak,intrinsic ) and (T duration , CR peak ), or between the intrinsic and observed properties of the light curves when we check the completeness of our observational data.We plot the scatters of T duration − σ T and CR peak − CR peak,intrinsic of the simulated data, as shown in Figure 20, and fit the scatters with power-law relations.The best fits give T duration = 2.7956 (σ T ) 1.0443 , (B4) CR peak = 1.0076 (CR peak,intrinsic ) 1.0012 , (B5) which are employed to obtain the completeness distribution of observational data when plotting Figure 17.The best fits suggest that CR peak is nearly the same as CR peak,intrinsic , while T duration is about 2.8 times of σ T .In this case, the full width at half maximum (FWHM) of the Gaussian-shaped light curve is relatively close to the duration given by the Bayesian blocks, i.e.T duration (see Sections 2.3 and 2.4).

Figure 2 .
Figure 2. Examples of the light curves for the eight categories: (a) Constant (C), (b) Flare (F), (c) Flare with single Bayesian Block (S), (d) Rise (R), (e) Decay (D), (f) Flare with multiple flaring phases (M), (g) Variable (V) and (h) Bad (B).Binned light curves and Bayesian Blocks are shown with black and orange lines, respectively.The observation identifiers (ObsIDs), equatorial coordinates (including RA and Dec) and bin sizes of corresponding light curves are given above each panel.Three time parameters, tstart, t peak and t end , of each flaring phase are illustrated with green, red and blue dash-dotted lines, respectively.

Figure 4 .
Figure 4. Distributions of the stellar flares in our sample with light curve categories F, M-F, S and M-S.(a) Equatorial distribution.(b) Distance distribution.(c) Distribution across a Hertzsprung-Russell diagram (HRD).When plotting the HRD, we use BP-RP colour from Gaia DR3 as the x-axis, and use absolute magnitude in the G-band of Gaia (calculated from the G-band mean magnitude and distance) from Gaia DR3 as the y-axis.The probability density in panel (b) is defined as the value of the probability density function at each bin, which is normalized such that the integral over the range is 1.Such a definition also applies for the probability density in the following figures.

Figure 5 .
Figure 5. Table3 δ = 1.31 ± 0.18(1σ), which is in good agreement with the scaling relation of E flare,SXR ∝ L ∼1.2 peak,SXR derived by Pye et al. (2015) from their XMM-Newton survey of ∼ 130 flares from ∼ 70 Tycho cool (F-M type) stars.The similar scaling indices indicate that F-M type stars, including solar-type stars, may share the same mechanisms responsible for the energy relase at the SXR band.

Figure 7 .
Figure 7. Distributions of time parameters for F and events.(a)-(c) Histograms and corresponding normal distributions of Trise,SXR, T decay,SXR and T duration,SXR in logarithmic scale.The most probable values (µ) and standard deviations (σ) are listed in the panel legends.The values of µ are illustrated with orange dash-dotted lines, and the values of µ − σ and µ + σ are illustrated with blue dotted lines in each panel.

Figure 8 .
Figure 8. Distribution of the ratio of decay time to rise time for F and M-F samples.Panel legend and the meaning of line styles are the same as in Fig. 7, but for log 10 (T decay,SXR /Trise,SXR).

Figure 9 .
Figure 9. Distribution of flare durations for all F, M-F, S and M-S events.Panel legend and the meaning of line styles are the same as in Fig. 7.

Figure 10 .
Figure 10.Distributions of peak luminosities (a) and flare energies (b).Panel legends and the meaning of line styles are the same as in Fig. 7.

Figure 13 .
Figure 13.Scatter plots of (a) peak luminosity and (b) duration as a function of flare energy.Error bars are plotted for each flare.Errors of L peak,SXR and E flare,SXR are estimated from the uncertainty of distance and the uncertainty of flux derived from the model when fitting the spectra (see Section 2.4).Orange dash-dotted lines represent the best-fitting results of the power-law correlations L peak,SXR ∝ E 0.76±0.11(1σδ ) flare,SXR and T duration,SXR ∝ E 0.201±0.024(1σβ ) flare,SXR , and the power-law indices are also shown in the legends of corresponding panels.

Figure 14 .
Figure 14.Scatter plot of the duration as a function of peak luminosity.Error bars are plotted for each flare, as in Fig. 13.

Figure 15 .Figure 16 .
Figure 15.Distributions of spectral indices of SXR spectra at the rise phase (a) and decay phase (b), as well as the index ratio (c).Panel legends and the meaning of line styles are the same as in Fig. 7.

Figure 17 .
Figure 17.Completeness distribution against T duration and CR peak .The completeness of the sample set at each (T duration , CR peak ) is illustrated by the color bar, which is generated from simulated light curves (see Section 3.6 and Appendix B).Three contours with completeness of 0.50, 0.70 and 0.90 are shown with red dashed lines.The distribution of the real observed flare events in this work is shown with blue crosses.

Figure 18 .Figure 19 .
Figure 18.Light curve for the flare event with ObsID 1872, RA 86.64745 and Dec 0.14953, which contains probably two or more unidentified flare events.The meaning of line styles is the same as in Fig. 2.

Figure 20 .
Figure 20.Scatter plots of (a) T duration − σT and (b) CR peak − CR peak,intrinsic of the simulated flare light curves.The color shows the areal density of the sample: the bluer parts are denser than the redder parts.The black dashed lines show reference lines of y = x, or where the "observed" values equal the "intrinsic" values.The best-fit power-law relations are shown in grey dash-dotted lines.When fitting the scatter plots, abnormally diffuse data points are excluded.