SNAD catalogue of M-dwarf flares from the Zwicky Transient Facility Open Access

ABSTRACT

Most of the stars in the Universe are M spectral class dwarfs, which are known to be the source of bright and frequent stellar flares. In this paper, we propose new approaches to discover M-dwarf flares in ground-based photometric surveys. We employ two approaches: a modification of a traditional method of parametric fit search and a machine learning algorithm based on active anomaly detection. The algorithms are applied to Zwicky Transient Facility (ZTF) data release 8, which includes the data from the ZTF high-cadence survey, allowing us to reveal flares lasting from minutes to hours. We analyse over 35 million ZTF light curves and visually scrutinize 1168 candidates suggested by the algorithms to filter out artefacts, occultations of a star by an asteroid, and other types of known variable objects. The result of this analysis is the largest catalogue of ZTF flaring stars to date, representing 134 flares with amplitudes ranging from −0.2 to −4.6 mag, including repeated flares. Using Pan-STARRS DR2 colours, we assign a spectral subclass to each object in the sample. For 13 flares with well-sampled light curves and available geometric distances from Gaia DR3, we estimate the bolometric energy. This research shows that the proposed methods combined with the ZTF’s cadence strategy are suitable for identifying M-dwarf flares and other fast transients, allowing for the extraction of significant astrophysical information from their light curves.

1 INTRODUCTION

M-dwarf stars make up the vast majority of stars in Galaxy. As low-mass, fully convective stars, they exhibit frequent flaring events caused by powerful magnetic reconnection processes in their atmospheres (Gershberg & Pikel’Ner 1972). The study of M-dwarf flares provided key insight into stellar magnetism, high-energy phenomena, and the impacts on potential planets orbiting these stars. However, many fundamental properties of M-dwarf flares remain poorly constrained, including their occurrence rates, energy output, and relationship with stellar properties, such as age and metallicity (see Kowalski 2024 for a review).

Multiple time-domain optical surveys have been utilized for systematic M-dwarf flare search projects. Space-based missions such as Kepler (Borucki et al. 2010) and Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014) offer high regular cadence and very precise relative photometry, making their data excellent sources of stellar flares, especially those of small amplitude. Yang & Liu (2019) discovered approximately |$1.6\times 10^5$| flares on around 3400 stars in the Kepler data. Günther et al. (2020) found 8695 flares in the first TESS data release (DR), while Pietras et al. (2022) refined this number to roughly |$1.4\times 10^5$| flares over 3 yr of TESS data. While these space missions have revealed a vast number of stellar flares with good completeness, ground-based surveys could complement them. For instance, the Zwicky Transient Facility (ZTF; Bellm et al. 2019) covers sky areas that Kepler did not observe, with a bigger survey depth than TESS.

The recent advent of wide-field time-domain surveys provided new opportunities to build large statistical samples of stellar flares across a range of spectral types. Numerous systematic stellar-flare searches were performed with different ground-based surveys, occasionally observing high-amplitude flares. For example, Davenport et al. (2012) analysed hundreds of flares from SDSS and 2MASS time-domain surveys with the maximum u-band change in magnitude |$|\Delta M|$||$\sim 4.5$| mag. The ASAS-SN M-dwarf flare catalogue contains 62 flares with a maximum |$|\Delta M|$| being |$\sim 2$| mag in the V band (Rodríguez Martínez et al. 2020). Webb et al. (2021) found 96 flares from 80 stars in the DECam data, with a maximum change in magnitude of |$\sim 1.8$| mag. Aizawa et al. (2022) revealed 22 fast stellar flares from a 1-s cadence survey performed by the Tomo-e Gozen project at the Kiso Schmidt telescope in Japan, among which the largest change in magnitude was of |$\sim 3.25$| mag. Liu et al. (2023) presented a catalogue of 132 flares from 125 stars with |$|\Delta M|$| up to |$\sim 3.1$| mag found in 2-yr data of the Tsinghua University-Ma Huateng Telescopes Survey.

In this work, we use photometric data from the ZTF survey to detect transient events. ZTF runs multiple surveys, including a high-cadence survey, which provides a unique data set of minute-scale cadence with |$\sim 21.5$| limiting magnitude (Kupfer et al. 2021). This makes ZTF DRs a source of well-sampled M-dwarf flare light curves occurring in a few hundred parsecs from the Sun. Crossland et al. (2024) used ZTF high-cadence survey for a systematic search for gravitational self-lensing binaries and presented 19 candidates. However, according to the authors, most of these candidates are likely to be stellar flares, making this data set the largest ZTF stellar flare catalogue previously published. Additionally, the SNAD team has discovered few stellar flares in ZTF DRs using machine learning anomaly detection pipelines (Malanchev et al. 2020; Pruzhinskaya et al. 2023; Volnova et al. 2023).

Various algorithms have been developed for detecting stellar flares in both sparsely sampled and continuously sampled data. Traditional flare detection studies often relied on detrending light curves and using outlier detection heuristics to identify flare events. These methods focused on preprocessing the data to remove trends and then identifying unusual data points which deviate significantly from expected patterns (e.g. Osten et al. 2012; Davenport 2016b; Günther et al. 2020). Another approach is parametric fitting, which models the light curves of flares using predefined mathematical functions to identify their characteristics (Lawson et al. 2019). However, the diversity of flare shapes and the data volumes of ongoing wide-field surveys encourage the community to use machine learning techniques. For instance, Feinstein et al. (2020) developed a convolutional neural network (CNN), stella, specifically trained to find flares in TESS short-cadence data. Another machine learning algorithm presented by Vida & Roettenbacher (2018) identified flares by analysing light curves and has been successfully applied to stars such as TRAPPIST-1 and KIC 1722506. This paper explores parametric fitting and machine learning methodologies, which leads to the detection of a heterogeneous set of flaring stars. Since not all the candidates selected with these methods are stellar flares, we also conduct the dedicated one-by-one analysis to filter out bogus detections and objects of a different astrophysical origin.

This paper provides a catalogue of 134 flaring stars detected using two different methodologies. We calculate astrophysical properties of stars, such as spectral class and interstellar reddening and analyse well-sampled light curves to estimate the total energy, amplitude, and time-scale of the flares.

2 DATA

The ZTF is a wide-field astronomical survey of the entire northern sky, conducted with the 48-inch Schmidt-type Samuel Oschin Telescope at Palomar Observatory (Bellm et al. 2019). During phase I (2018 February–2020 September), ZTF performed a 3-d cadence survey of the visible northern sky and 1-d for the Galactic plane. During phase II (2020 October–2023 September), 50 per cent of the ZTF camera time is dedicated to a 2-d cadence public survey in g and r bands. Data from the public survey are released on a bi-monthly schedule as DRs. In addition, Kupfer et al. (2021) conducted a dedicated high-cadence Galactic plane survey with a cadence of 40 s.

In this work, we analyse the private and public data from ZTF DR8 (2018 March 17–2021 September 3) as target data sets for searching for red dwarf flares. We use only ‘clear’ (catflags = 0) r- and g-band observations. ZTF DRs provide unique objects in each passband and observational field, so our data set may have multiple light curves associated with a single source. In order to conduct further astrophysical analysis of the resulting candidates, we use data from additional catalogues: riz magnitudes from Pan-STARRS DR2 (Chambers et al. 2019; Flewelling 2018), geometric distances from Gaia EDR3 (Gaia Collaboration 2016, 2021; Bailer-Jones et al. 2021 ), interstellar extinction from ‘Bayestar19’ (Green et al. 2019), and Schlafly & Finkbeiner (2011) extinction maps.

3 METHODS

In this paper, we use two distinct approaches of M-dwarf flare identification. The first method is based on parametric search of flaring light curves in high-cadence subsample of ZTF Data Release 8 (DR8). The second method is based on an active anomaly detection approach applied to the full light curves of the entire ZTF DR8. We describe both methods in detail below.

3.1 Parametric fit search

As part of the ZTF survey, several different observational campaigns were conducted, including high-cadence survey (Kupfer et al. 2021). Unfortunately, bulk-downloadable ZTF DRs do not maintain the connection between each individual observation and the specific campaign to which this observation belongs. For this reason, it was necessary to develop a method for extracting high-cadence data from an entire DR.

The aim was to chunk light curves to form a data set of intra-day light curves having: (1) enough observations for further analysis, (2) high cadence, and (3) covering the time interval typical for flare duration. To achieve the first two conditions, we set the minimum number of observations to 5 and the maximum delay between two consecutive observations to 30 min. As for the minimum duration of these partial light curves, we produce two samples, one with minimum duration of 2 h (long-duration sample) and another with 30 min (short-duration sample). Such an approach allows us to detect both short and long flares.

Along with the imposed conditions to the cadence and duration, all light curves were filtered according to observed source variability. We consider an object variable if a test based on one-dimensional reduced |$\chi ^2$| statistics rejects the hypothesis concerning its non-variability (Sokolovsky et al. 2017):

where |$m_i$| is the observed magnitude and |$\delta _i$| its observational error, |$\bar{m}$| is the weighted mean magnitude, N is the number of observations, and d the number of model parameters. The final data set consists of light curves with a value of the reduced |$\chi ^2$| statistics greater than 11. The total number of intra-day r-band light curve chunks in the long-duration sample is 4027 686, and in the short-duration one is 10 351 985.

The parametric search method is based on light-curve fitting with an analytical function and subsequently selecting well-fit objects. For this purpose we adopt a semiphenomenological model of flux evolution, |$f(t)$|⁠, from Mendoza et al. (2022):

where |$f^{*}$| is stellar (quiescence) flux density, t is time, A is the normalizing factor, B is the reference time, C is the Gaussian heating time-scale, |$D_1$| is the inverse of the rapid cooling phase time-scale, |$D_2$| is the inverse of the slow cooling phase time-scale, |$F_2 \equiv (1-F_1) \cdot \exp {(\frac{C^2}{4}(D^2_1-D^2_2))}$|⁠, where |$F_1$| and |$F_2$| describe the relative importance of the exponential cooling terms, and |$\mathrm{erfc}$| is the complimentary error function. Note that A corresponds to the value used by Mendoza et al. (2022) multiplied by |$\exp {(B^2/C^2-D_1^2 C^2/4)}$|⁠. We also change the form of the equation by introducing the dimensionless function |$\alpha (t,C,D)$| for better readability and robustness of the fit.¹

A python function for fitting a light curve was implemented within the light-curve² feature extraction library (Malanchev 2021). This function chooses the optimal values of parameters A, B, C, |$D_1$|⁠, |$D_2$|⁠, |$F_2$|⁠, and |$f^{*}$| using least-squares fits provided by iminuit (Dembinski et al. 2024). For better performance of the least-squares fitting, we used best-fitting coefficients from the Markov chain Monte Carlo analysis, presented in Mendoza et al. (2022), as initial values. The mean value of the |$3\sigma$|-clipped flux is used as the initial value of stellar flux |$f^{*}$|⁠. The fit quality is evaluated using reduced |$\chi ^2$| statistics (equation 1) with |$d = 7$|⁠.

The next filtering step is to distinguish flares with enough points in the flare from ones with only a few points. For this reason, we use the OtsuSplit feature of the light-curve package (Lavrukhina, Malanchev & Kornilov 2023). This feature uses a magnitude threshold to distinguish faint and bright subsamples of a light curve based on maximization of interclass variance. To filter out candidates, we use |$\tt {lower-to-all} \gt 0.25$| (ratio of the number of points in the bright subsample to the number of points in the faint one) obtained for each object based on the determined threshold. However, this method does not guarantee perfect separation of the ‘flaring’ part of the light curve from the ‘plateau’, so some candidates with few-point flares are still present in the final sample. The total number of candidates obtained with this procedure is 308.

3.2 Machine learning method

Active anomaly detection represents a family of machine learning techniques which sequentially uses expert feedback to fine-tune an initially standard unsupervised algorithm to a particular definition of scientifically interesting anomaly. In the implementation used in this work, we employ the Active Anomaly Discovery (AAD) algorithm developed by Das et al. (2018) in the form used by Ishida et al. (2021).

The algorithm starts from a traditional isolation forest (IF; Liu, Ting & Zhou 2008). It is based on the hypothesis that objects in underdense regions of the parameter space (statistical anomalies) are more rapidly isolated from the bulk of the data than nominal ones. In the first step of AAD, a traditional IF is built and the object with the highest anomaly score is shown to a human expert, who is required to provide a binary answer (‘YES’/‘NO’) to the question: is this anomaly scientifically interesting to you? The expert makes decisions based on both light-curve behaviour and auxiliary data, such as original scientific images and catalogue cross-matches (see Section 4 for the details). If the answer is YES, the algorithm will show the second object with the highest anomaly score and pose the same question. Alternatively, if the expert answers NO, the algorithm recalculates the weights corresponding to the decision path which lead to that anomaly score. This modification is applied to the entire forest, the scores are recalculated, and the new object with the highest anomaly score is shown to the expert. After a few iterations, this procedure results in a personalized model which has a lower probability to give high anomaly scores to objects which are not in the expert’s main interest (full mathematical description is provided in Das et al. 2018). The flowchart that demonstrates the main steps of the AAD approach is presented in Fig. 1.

Since the algorithm can be adapted to the opinion of the expert, it can be used for a targeted search for objects of a given type (see Pruzhinskaya et al. 2023). Therefore, in this analysis, a human expert considers only M-dwarf flare candidates as anomalies; all other objects proposed by the algorithm are rejected by the expert as nominals.

The minimum number of observations per light curve was selected to be 300. To avoid non-variable objects, we select |$\chi 2 \gt 3$| (equation 1, |$d=1$|⁠). This selection resulted in 21 469 857 light curves.

We run AAD on the obtained data set and visually inspect 860 objects. Among those, we responded ‘YES’ to 35 objects.

4 SAMPLE SELECTION

Visual inspection by human experts is a part of the sample selection process for both methods: as a final check of the parametric fit outputs, and at each iteration step of the AAD algorithm. All candidates were inspected using the SNAD ZTF Viewer,³ a special web interface developed by the SNAD team to facilitate expert analysis (Malanchev et al. 2023). For each object, the Viewer displays its multicolour light curves and enables easy access to the individual exposure images and to the Aladin Sky Atlas (Bonnarel et al. 2000; Boch & Fernique 2014). Moreover, it provides cross-matches with various catalogues of stars and transients, including simbad and VizieR data bases (Wenger et al. 2000), AAVSO VSX (Watson, Henden & Price 2006), Pan-STARRS DR2 (Flewelling et al. 2020), Gaia DR3 (Gaia Collaboration 2023), and ZTF alert brokers.^4,5,6

When deciding whether an event belongs to a red dwarf flare, we applied the following criteria:

• The light-curve profile should be typical for a red dwarf flare: sharp increase in brightness followed by a smooth decrease. For example, the short-duration plateau observed during a ‘flare’ is typically explained by the occultation of a star by an asteroid (left panel of Fig. 2).

• Absence of an artefact at the frame, like satellites, frame edges, defocusing, ghost, bright star nearby, cosmics or hot pixels (right panel of Fig. 2).

• The object should not be a known variable star, whose observed changes in brightness are caused by other processes.

Based on these criteria, we scrutinized 343 objects of both the parametric fit and the machine learning candidates. This resulted in a final data set of 134 flaring M-dwarf stars. Only two flares (762109400005614 and 764114400003060) were found in ZTF g band, the rest were identified in ZTF r band. The remaining 209 candidates turned out to be artefacts (75 objects), instances of a star being occulted by an asteroid (128 objects) and parts of short-periodic star light curves (six objects).

The final sample of M-dwarf flares and their main characteristics are presented in Table A1. For M-dwarfs with recurrent flares, only the flare with the largest number of photometric measurements is listed. The first column contains the ZTF DR object identifiers (OIDs). The equatorial coordinates (⁠|$\alpha$|⁠, |$\delta$|⁠) are presented in the second and third columns, respectively. The fourth column contains the geometric distance of each object estimated with Gaia EDR3 and its uncertainty (Bailer-Jones et al. 2021). In the fifth column, the object’s |$A_r$| extinction is given (Schlafly & Finkbeiner 2011; Green et al. 2019). The sixth, seventh, eighth, and ninth columns contain peak time, full width at half maximum (FWHM), amplitude, and number of points in the flare (see Section 5.2). The 10th column presents our spectral class prediction based on Pan-STARRS photometric colours. The last column contains notes. The next section provides detail of the analysis which led to columns four to nine.

5 ANALYSIS

5.1 Flare energy

For further energy analysis, we use geometric distances derived from Gaia EDR3 (Bailer-Jones et al. 2021). Many candidates were found to be long-distant stars, for which distances have high uncertainties. We calculate the flare energy for the subsample of candidates whose Gaia DR3 parallax was measured with uncertainty not higher than 20 per cent. We also keep candidates with enough points in the flare for higher quality of flare profile fitting. The final subsample consists of 13 flares, which were used for the energy calculation.

We assume that flare radiation could be described by an optically thick blackbody with a temperature of |$T_\text{flare} = 9000$| K (Hawley & Fisher 1992). However, we acknowledge that this assumption is rather simplistic, as the optical continuum spectrum is significantly more complex, varying from flare to flare and from the impulsive to the gradual phases of each flare (Kowalski et al. 2013, 2019; Brasseur et al. 2023). Additionally, there is some |$\textrm {H}\,\alpha$| line contribution to the r band during the gradual decay phase (Hawley & Pettersen 1991). None the less, we use this model to maintain consistency with the early assumptions made for calculating flare energies in Kepler (Shibayama et al. 2013), so the bolometric luminosity, |$L_{\textrm {flare}}(t)$|⁠, is given by

where |$\rm {\sigma _{SB}}$| is the Stefan–Boltzmann constant, and |$\mathcal {A}_\text{flare}(t)$| is the flare surface area that changes over time.

Now our objective is to estimate |$\mathcal {A}_\text{flare}(t)$| from the observed flare flux. First, we introduce |$\mathcal {A}_\perp$|⁠, a projection of |$\mathcal {A}_\text{flare}$| in the picture plane. Then, the spectral flux density of the flare |$F_{\nu } = \Omega B_\nu$|⁠, where |$\Omega = \mathcal {A}_\perp /d^2$| is the solid angle of the flare as observed from distance d, and |$B_\nu$| is the blackbody intensity–Planck function.

We did not directly observe the spectral flux density with photometric surveys. Instead, its value was averaged over the passband transmission of the photometric filter in use. The averaged spectral flux density in ZTF r band, |$F_r$|⁠, is given by Koornneef et al. (1986):

where R is the filter transmission function.

Using equation (6), we got |$\mathcal {A}_{\perp }(t)$|⁠:

where |$B_r$| is the blackbody intensity averaged over r passband. Being combined with equation (5) it gives the final expression of bolometric luminosity:

Since we do not know the shape of the flare, we introduce the geometric factor |$\mathcal {A}_\text{flare}/\mathcal {A}_\perp$|⁠. To be consistent with previous studies (Shibayama et al. 2013; Yang et al. 2017), we assume that this factor does not change with time and equals one.

Finally, we arrive at the expression for the bolometric energy:

This integral value depends only on the filter-averaged spectral flux density |$F_r(t)$|⁠, which can be derived from observed magnitudes m. We convert observed magnitudes m at time moments |$t_i$| to fluxes taking into account extinction |$A_r$| given by the 3D map of Milky Way dust reddening (Green et al. 2019):

where |$m_0$| is the AB-magnitude zero-point. To make the integral value more robust to observation uncertainties, we fit observations with the parametric function (2) and use these fitted model fluxes as a proxy to the filter-averaged spectral flux density: |$F_r(t) = f(t)-f^{*}$|⁠.

The energy calculation results are given in Table 1.

5.2 Flare fit parameters

We use the same parametric function from (2) to estimate the amplitude, FWHM, and number of points in the flare (Table A1). First, a manual operation of flares’ observations subtraction from the full light curve was conducted – each stars’ light curve was trimmed to capture the flare itself and some amount of quiescent points which are necessary for the further fitting operation.

Secondly, for each flare, its profile was fitted using a parametric function (equation 2) to get a continuous representation convenient for further analysis. For complex light curves with several flare peaks, only the peak with the maximum amplitude (main peak) is considered for further fitting. The amplitude is calculated as the difference between the maximal and minimal model flux, which is achieved over a time interval. As a next step, the FWHM is measured as the difference between time points where the model curves possess values equal to half of the defined amplitude.

In order to define the number of points in the flare, we adopted the following criterion: all points of the light curve whose observed flux exceeds |$f^{*}$| (quiescence stellar flux obtained from the model) by |$3\sigma$|⁠, where |$\sigma$| is the mean observational error of the selected part of the light curve, were considered as belonging to the flare.

5.3 Spectral class

Due to the lack of spectral information for most of our objects, we employ a photometric approach to estimate their spectral types. We use stacked magnitudes from Pan-STARRS DR2 (Chambers et al. 2019; Flewelling 2018) to build the (⁠|$r-i$|⁠, |$i-z$|⁠) colour–colour diagram of our sample. The method is adopted from Kowalski et al. (2009) and allowed us to define the spectral subtypes of M-dwarfs based only on their photometric colours.

First, we correct the colours for the galactic reddening. We use extinction values and coefficients from two different sources: the three-dimensional map of Milky Way dust reddening ‘Bayestar19’ (Green et al. 2019), and, if this map does not contain the object, the map of Galactic Dust Reddening and Extinction presented by Schlafly & Finkbeiner (2011). Based on these maps, we calculate the final extinction value using

where |$\boldsymbol{A} = (A_r, A_i, A_z)$| represents the final extinction value in each filter, |$\boldsymbol{R} = (R_r, R_i, R_z)$| is the ‘extinction vector’, relating a scalar reddening to the extinction in each passband and E is the reddening in the dust-map specific units. The extinction vector is |$\boldsymbol{R} = (2.617, 1.971, 1.549)$| for the 3D map ‘Bayestar19’, and |$\boldsymbol{R} = (2.271, 1.682, 1.322)$| for the 2D map of Galactic Dust Reddening and Extinction. For the 3D map, we use Gaia EDR3 geometrical distance as discussed above.

Kowalski et al. (2009) proposed a table with the best-fitting parameters of the two-dimensional Gaussian distribution for each M-dwarf spectral subtype in the (⁠|$r-i$|⁠, |$i-z$|⁠) colour space. We calculate the probability of belonging to the corresponding spectral subclass following their proposal,

where m is an M-dwarf spectral subclass index, from 0 to 7, |$\boldsymbol{\mu }_m$| and |$\Sigma _m$| are Gaussian mean values and covariance matrix for this subclass, |$\boldsymbol{x} = ({r-i}, {i-z})$| are deredded colours of the studied star.

For each star, we assign the spectral subclass m corresponding to the maximum probability |$p_m$|⁠. We show the objects in the colour–colour diagram (Fig. 3) with subclass corresponding to the point estimations of the pair of object’s colours. However, the deredded colour values may have large uncertainties associated with Pan-STARRS photometry error, distance estimation error (applicable only for the 3D dust maps), and extinction errors.

Finally, we apply a colour cut based on West et al. (2008) and Bochanski et al. (2007), to construct a more pure M-dwarfs sample. Kowalski et al. (2009) presented the following limits, which take extinction into account: |$(r-i) \gt 0.42$|⁠, |$(i-z)\gt 0.23$|⁠.

Starting from the 134 stars found, we determine the spectral subclasses of 132 of them (see Table A1, Fig. 3). Five objects that lie outside of the M-dwarf limits might have a different spectral type. Their OIDs are 435211200068171, 540215200069194, 687207100049742, 687214100050598, and 768209200100383, with corresponding |$A_r$| values of 1.14, 5.11, 1.15, 2.82, and 2.67, respectively. It is possible that the interstellar reddening for these objects has been overestimated: four of them do not have Gaia distances, and for another, the distance is measured with high uncertainty, which may cause the object to appear shifted towards bluer colours.

For the 16 brightest objects, Gaia DR3 provides effective temperature estimations (Gaia Collaboration 2023). We use these estimations to define a spectral class based on the relationship between the stellar effective temperature and the spectral class for main-sequence stars adopted from Malkov et al. (2020). In all cases, this method confirms the spectral class M for these objects. For seven objects, we have an exact agreement for the spectral subclass between this method and the method based on Pan-STARRS colours. There is also one object (719206100008051) with available spectrum from SDSS DR18 (Almeida et al. 2023) data and the corresponding spectral classification, which is consistent with our colour-based classification. See the individual object information in the ‘note’ column of Table A1.

It should be noted that using Gaia geometric distances, we can estimate the absolute magnitude of our objects. The objects having much brighter absolute magnitude than expected for main-sequence M-dwarf stars (Pecaut & Mamajek 2013) also have large parallax uncertainties. We mark these objects with crosses in the distance column of Table A1. Although most of those with more confident distance estimations are consistent with the main sequence, the outliers may be explained by system multiplicity or other systematic errors.

6 DISCUSSION

6.1 Flares morphology

The flaring stars we study in this paper vary in the number of observational points per flare, flare recurrence, and light-curve profiles. The light curves of all the flares are provided in the supplementary material.

For some of the flaring M-dwarfs multiple flares are observed: either multiple flares in a single passband (Fig. 4a), or a single flare simultaneously observed in g and r bands (Fig. 4b).

Additionally, we distinguish between classical and complex flare events on the basis of their temporal structure (e.g. Kowalski et al. 2013; Davenport et al. 2014). Classical flares have a single-peak profile, characterized by fast rise and exponential decay (Fig. 4c). However, the majority of our flares display a much more complex structure. This complexity ranges from relatively simple flares (see Fig. 4d) to highly complex flares consisting of multiple components (see Fig. 4e). According to Davenport (2016a), studying such flare complexity could clarify their origins, as it remains uncertain whether they are produced by a single active region or by triggering separate nearby regions. In the first scenario, Davenport (2016a) anticipates that the most energetic flare occurs first, followed by a sequence of less energetic events. However, this is not supported by some flares in our sample, where a less energetic peak precedes the most energetic one (Fig. 4f).

Despite our best efforts, we recognize that there is still some inherited uncertainty in the classifications presented here, especially when only one point is available. They could be associated with hotspots on a stellar surface, self-lensing binaries or other types of stars that flare (e.g. Kowalski 2024). As it was recently stated by Crossland et al. (2024), a possible scenario would be self-lensing detached binaries, containing a stellar-mass neutron star or a black hole. The brightening occurs when the compact object transits the companion star. In that case, a symmetric light curve should be observed. Among our flares, there are a few candidates which satisfy this criteria (e.g. Fig. 4g).

6.2 Parametric fit versus AAD

Searching for flares on M-dwarfs is a complex task. Each data set, whether it comes from different surveys or not, is unique. That is why it is so important to explore various searching strategies. In this paper, we apply two different methods for M-dwarf flares search. Below we compare both approaches.

First of all, the parameters associated with the discovered flares are different. Due to the ZTF’s sporadic observation schedule, flares found with AAD are sampled from incomplete light curves, often missing the peak brightness. In contrast, a parametric fit was applied to high-cadence data, thereby providing a well-defined light-curve profile. Consequently, flares found with AAD algorithm have systematically smaller amplitudes in comparison with the ones from the parametric fit (see Fig. 5a). The number of points in AAD flares is significantly lower, which makes them less reliable and more difficult to confirm. Also, recurrent flares were more often found by the AAD approach, since it had access to long-duration light curves (e.g. Fig. 4a).

Second, regarding the spatial distribution of flares, the parametric fit search was limited by a ZTF high-cadence coordinate cut, resulting in flares that are located within the Milky Way plane. In contrast, the AAD does not have any coordinate restriction, yet the observed bias towards higher galactic latitudes in the Fig. 6 can be attributed to the selection effect, i.e. a fewer runs of the algorithm in fields within or close to the galactic plane, since this method requires expert evaluation at every step of its operation. Moreover, flares obtained by the parametric fit method have a larger distance to Earth according to Gaia EDR3 (Fig. 5b). This can be explained by the use of a reduced |$\chi ^2$| metric to evaluate goodness of fit, which includes photometric errors, therefore lower metric quantities are systematically associated with fainter and more distant sources.

Finally, since AAD is an anomaly detection algorithm, it was not originally intended for flare detection. However, by treating flares as anomalies, we are able to successfully adapt it to search for flares. It is interesting to note that only one over a few high-cadence flares identified by AAD exhibit exceptionally high amplitudes when compared to typical flare energy outputs (Fig. 4, see also Rodríguez Martínez et al. 2020; Gorbachev & Shlyapnikov 2022). On the other hand, the parametric fit technique excels in isolating well-sampled flares within high-cadence amplitudes, enabling the observation of complex substructures within the light curves (e.g. Fig. 4e). It allows performing a comprehensive analysis of individual flares (Davenport 2016a).

7 CONCLUSIONS

This paper is devoted to the search for M-dwarf flares in the eighth data release of the ZTF survey. We explore two different approaches: a parametric fit search and a machine learning method.

We visually scrutinized 1168 candidates, filtering out artefacts and known variables of other types, to identify 134 M-dwarf flares. This constitutes the largest sample of ZTF M-dwarf flares identified to date. The associations with the M spectral class are confirmed through the (⁠|$r-i$|⁠, |$i-z$|⁠) colour diagram analysis, though some classifications may be incorrect, and we note opportunities for identifying exotic events like self-lensing binaries. For 13 objects, we calculate the flare energy, ranging from |$\sim 7\times 10^{33}$| to |$\sim 404\times 10^{33}$| erg, which is consistent with the higher end of the energy distribution reported in the literature (Yang et al. 2017).

The comparison between the two approaches shows that each of them identifies flares of different parameters and distribution in the sky. For example, the parametric fit search found fainter flares, while AAD, despite lagging in flare light curve quality, identified recurrent flares. Additionally, the highest amplitude flare in the sample was discovered using AAD. Since each method has its own limitations, diverse strategies for flare detection are necessary to form a comprehensive picture of these phenomena.

Although the ZTF survey is not specifically designed for fast transients due to its 2–3 d cadence, it conducts private high-cadence observational campaigns. Such campaigns are also envisaged by the observational strategy of the Vera Rubin Observatory Legacy Survey of Space and Time with Deep Drilling Fields program (LSST Science Collaboration 2009). For the search and study of red dwarfs, we should not rely solely on dedicated surveys; instead, we must learn to extract necessary information from surveys not originally intended for this purpose. Therefore, developing methods for data filtering and flare identification is highly relevant.

ACKNOWLEDGEMENTS

AL, MK, AV, and TS acknowledges support from a Russian Science Foundation grant 24-22-00233, https://rscf.ru/en/project/24-22-00233/. Support was provided by Schmidt Sciences, LLC. for KM. VK acknowledges support from the youth scientific laboratory project, topic FEUZ-2020-0038. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under grant no. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation grant no. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation. This work has made use of data from ZTF, supported by the National Science Foundation under grants no. AST-1440341 and AST-2034437 and a collaboration including current partners Caltech, IPAC, the Oskar Klein Center at Stockholm University, the University of Maryland, University of California, Berkeley, the University of Wisconsin at Milwaukee, University of Warwick, Ruhr University, Cornell University, Northwestern University and Drexel University. Operations are conducted by COO, IPAC, and UW.

DATA AVAILABILITY

The ZTF light-curve data underlying this article are available in the NASA/IPAC Infrared Science Archive.⁷ The final table with selected flares and their parameters are given in Appendix  A.

Footnotes

REFERENCES

APPENDIX: TABLE OF M-DWARF FLARES

We show here the table of all found flares and the corresponding stars with their main characteristics.

Author notes