© The Authors 2025

1 Introduction

Stellar flares are stochastic and intense bursts of electromagnetic radiation and charged particles that are triggered by magnetic reconnection in the chromosphere of stars. This phenomenon occurs when the magnetic field lines in the chromosphere realign, which releases energy across the electromagnetic spectrum. This can generate flares of varying intensity and duration and can significantly increase the stellar luminosity (Benz & Güdel 2010). The stellar magnetic activity significantly influences the frequency and energy of stellar flares. M dwarfs are known for their high levels of magnetic activity and for flaring more frequently than earlier-type stars because their dynamo processes are more vigorous and yield a higher incidence of magnetic reconnection events (Kowalski 2024). These stars can experience more frequent and energetic flares that sometimes reach energy levels that are several orders of magnitude more intense than solar flares (Bruevich et al. 2001; Howard et al. 2018). The effects of stellar flares on planetary habitability are various. On the one hand, stellar flares can damage planetary atmospheres: X-ray and ultraviolet radiation, along with charged particles that are emitted during flares, can deplete atmospheric layers, especially on planets with weak magnetic fields or located near their host star. This atmospheric erosion can compromise the potential for habitability by reducing the protection against harmful radiation and altering the planetary climate (Ribas et al. 2016; Airapetian et al. 2017). On the other hand, the energy that is radiated during stellar flares may play a role in triggering prebiotic chemical processes: The increased ultraviolet radiation during flares can trigger chemical reactions that are crucial for the emergence of life. This suggests that flares may be necessary for prebiotic chemistry to happen on planets orbiting M dwarfs (Scalo et al. 2007; Rimmer et al. 2018).

It is essential to understand the formation mechanisms, frequency, and energy distribution of stellar flares for several reasons. It provides insights into stellar magnetic activity and its variation with spectral type, especially in comparison to our Sun. By investigating the flare activity of various stellar types, it is additionally possible to refine models of stellar and planetary evolution (Güdel 2007; Lammer et al. 2012). TESS (Transiting Exoplanet Survey Satellite) and Kepler have provided valuable data on flare occurrences with cadences of minutes to a dozen seconds. The 20-second cadence of TESS revealed that complex flare structures are frequent, including multi-peak profiles and quasi-periodic pulsations (QPPs). These complex structures suggest that these non-classical flare shapes may result from emission cascades within a single active region or from sympathetic flares from neighbouring regions (Davenport 2016; Howard & MacGregor 2022). In this context, the study of smaller-scale flare events becomes increasingly important for a better characterisation of stellar variability and the formation mechanisms of flares. Extending photometric precision and observing cadence down to a few seconds could greatly improve our understanding of these phenomena and might reveal even finer flare structures. This increased sensitivity may enable us to detect micro- and nano-flares, which are small-scale flaring events with energies ranging from 10²² to 10²⁷ erg that are thought to play a crucial role in coronal heating (Parker 1988). No concluding evidence supports micro- and nano-flares against other possible explanations, however, mainly because it is hard to gather precise statistics about these flares because they are hidden in the instrumental noise. It is therefore still unclear whether micro- and nano-flares follow the same trends as more powerful flaring events. Furthermore, flare formation scenarios that diverge from self-organised criticality propose that flares might consist of an avalanche of several combined micro- and nano-flares. This would have implications for the current understanding of the impact of flares on exoplanetary atmospheres (Audard et al. 2000; Sheikh et al. 2016; Aschwanden 2019).

The space telescope CHaracterising ExOPlanet Satellite (CHEOPS) is primarily designed to measure exoplanet transits. It offers exceptional photometric precision. With a precision of approximately 20 ppm over a 6-hour integration period for a V ∼ 9 star and down to 150 ppm for shorter timescales, CHEOPS surpasses other photometric instruments in the visible spectrum, such as TESS and current ground-based telescopes (Benz et al. 2021; Oddo et al. 2023). The exposure time of CHEOPS can also reach 0.001 seconds for very bright stars (V ∼ 6), which is several orders of magnitude faster than those of TESS and Kepler (Borucki et al. 2010; Ricker et al. 2015). Although CHEOPS is limited to on-target observations, it provides high photometric precision and a fast accessible observing cadence. These are crucial for capturing the short and stochastic nature of low-energy flares.

To enhance the detection of flares in the CHEOPS light curves, denoising techniques based on the wavelet transform can filter out noise from the signal. The wavelet transform decomposes time-series data into different frequency components, which facilitates the isolation and reduction of noise. Unlike the Fourier transform, which only provides frequency characteristics of the signal, the wavelet transform recovers both temporal and frequency information, making it more efficient in analysing transient signals such as stellar flares. The discrete wavelet transform (DWT) is commonly used as a base for time-series denoising algorithms, and it has many applications outside of astronomy (Pasti et al. 1999; Polat & Özerdem 2018; Jang et al. 2021). The main downside of the DWT is that it introduces a downsampling of the input signal by definition, which can lead to a loss of details due to its non-redundant nature (Donoho & Johnstone 1994). In contrast, the stationary wavelet transform (SWT) avoids this by upsampling the wavelet coefficients at each decomposition level. Although this involves a higher computational cost, the SWT provides the advantage of maintaining the original length of the signal at each decomposition level and providing a redundant, shift-invariant representation of the data (Rhif et al. 2019; Kumar et al. 2021). In the context of astronomy, techniques based on the SWT have been developed for noise filtering and demonstrated significant improvements in the signal-to-noise ratios (S/N) and the ability to recover faint astrophysical signals (Starck & Murtagh 2006). These methods emphasise the importance of adapting the current wavelet-based denoising techniques to the specific characteristics of astronomical time-series and show that these denoising algorithms could enhance the detection of faint transient phenomena such as low-energy flares.

An efficient wavelet-based denoising algorithm is crucial for optimising the extraction of meaningful flare statistics in our analysis pipeline using CHEOPS data. We explore the performance of the combination of the photometric precision of CHEOPS and an SWT-based denoising algorithm in enhancing the flare detection in M dwarf light curves. In Section 2 we describe the target selection, the data reduction, the SWT parameter exploration, the flare detection, and flare-breakdown algorithms. In Section 3, we quantify the improvement in flare recovery that is achieved by denoising, we analyse the obtained flare catalogues, and we compare our results to the literature. In Section 4, we discuss the limitations of this study and suggest several avenues for improvement. We conclude in Section 5.

2 Methods

2.1 Target selection and data reduction

Among the data gathered by CHEOPS since the start of its operations, several observational programmes have targeted late-type stars. CH_PR100018 (PI I. Pagano) is a prominent programme that has been actively investigating the variability of late-K and M dwarf stars since the beginning of scientific observations in April 2020. This programme is particularly focused on characterising the photometric variability of these stars, including flares and star spots. CH_PR100010 (PI G. Szabó) operated from July 2020 to December 2023. While its primary aim was to study the debris disks around various spectral types, this programme included observations of several late-type stars to explore the interaction between stellar activity and circumstellar dust. CH_PR130057 (PI G. Szabó) was another significant programme, running from August 2022 to September 2023, and specifically targeted the young M dwarf star AU Mic. Known for its strong stellar activity and prominent debris disk, AU Mic is an ideal subject for examining the correlation between stellar flares and its surrounding environment.

The nominal observation mode of CHEOPS produces sub-arrays (200-pixel diameter) every few tens of seconds by co-adding multiple short sub-exposures onboard. In parallel, smaller cutouts called imagettes (60-pixel diameter) are down-linked at a higher cadence and can be used to construct light curves with finer temporal resolution. While individual sub-exposures can be as short as 0.001 seconds, the fastest cadence for downlinked subarrays is 22.65 seconds, set by telemetry and onboard processing constraints. In contrast, imagettes can be downlinked with cadences as short as 2.8 seconds. All programmes mentioned above adopted a 3-second cadence for imagettes to balance S/N, avoid saturation, and maximise duty cycle (Boldog et al. 2023; Bruno et al. 2024). This cadence remains well-suited for detecting and analysing rapid, transient events like stellar flares and offers higher time resolution than the 20-second cadence of TESS. To take advantage of the high precision and fast observing cadence provided by CHEOPS, we retrieved all publicly available imagettes data on main sequence M dwarfs from the CHEOPS archive¹. Due to the proprietary period associated with CHEOPS programmes, imagettes data are not yet available for all observations. We therefore completed our sample with the raw light curves provided in Bruno et al. (2024). In total, we obtained 66.65 days of on-target time from the three programmes, comprising observations of 110 different M dwarfs. For each target, we retrieved the Gaia G-band magnitude (G_mag) and effective temperature (T_eff) from the CHEOPS file metadata. We searched for H_α equivalent width values from Schöfer et al. (2019), as well as $\text{log\hspace{0.17em}}{R}_{HK}^{\prime }$ values from Astudillo-Defru et al. (2017) and Boro Saikia et al. (2018), but found H_α EW values for only 64 of the 110 stars and $\text{log\hspace{0.17em}}{R}_{HK}^{\prime }$ values for 71 stars. This limited availability prevented us from providing a comprehensive activity indicator distribution for the entire sample. Instead, we used the mean rotational velocity (V sin i) values from Table A.1 of Bruno et al. (2024) as an alternative indicator of stellar activity. We calculated the distance to each target using its parallax from Gaia DR3, or Gaia DR2 if not available. We obtained the radius of each star from the TESS Input Catalogue. The parameter distributions in our sample are shown in Figure 1, and the full list is available in Table A.1. As expected from the focus of the three observing programmes, our sample is heavily biased towards AU Mic, with more than 340 visits observed, while all other stars have been observed during less than 70 visits each.

The CHEOPS Data Reduction Pipeline (DRP) does not automatically process imagettes datasets (Hoyer et al. 2020). We hence relied on PIPE², a photometric extraction package developed by the CHEOPS consortium to complement the official DRP, to reduce imagettes and extract Point Spread Function (PSF) photometry (Brandeker et al. 2024). We optimised the PSF parameters for each target by allowing the software to freely explore the impact of the radius of the fitting region, the number of components for the fit, a simultaneous fit for the background, dark current correction, and static image removal through a gradient-like search algorithm. In the case of the light curves recovered by Bruno et al. (2024), this imagettes reduction was already performed using an ad hoc tool, the Data Reduction Tool (DRT), described in Morgado et al. (2022) and Fortier et al. (2024). We consistently found marginal differences of less than 50 ppm when comparing raw light curves extracted through both PIPE and the DRT.

Fig. 1 Parameters of our stellar sample. From top to bottom: Histograms for the Gaia magnitude, Teff, log(v sin i), radius, and distance.

2.2 Detrending

Variability in fast-cadence light curves can arise from multiple sources (Sulis et al. 2020; Kálmán et al. 2025). In the case of CHEOPS, white noise primarily results from photon counting, while time-correlated variations originate from both instrumental and astrophysical effects. To mitigate instrumental systematics, we used the package PyCHEOPS³ (Maxted et al. 2022, 2023), a tool developed by the CHEOPS consortium to remove correlations between flux and various instrumental parameters. Specifically, we removed linear, quadratic, and cubic flux correlations with time, centroid position, sine and cosine of the roll angle, background flux, contamination, and smear. Among astrophysical contributors, granulation in M dwarfs is expected to be negligible at CHEOPS precision (Sulis et al. 2023). The dominant source of residual variability can thus be attributed to surface inhomogeneities associated with magnetic active regions (such as spots, faculae, pores, etc) that induce rotational modulation. To account for this, we applied a Savitzky-Golay (Savgol) filter with a 20-minute window, corresponding to the duration of the longest flare detected in Bruno et al. (2024). This approach effectively removes trends on timescales longer than 20 minutes while preserving shorter-term variations caused by flares. Since long-duration flares are typically prominent and unlikely to have been missed in Bruno et al. (2024), this method ensures minimal loss of flare signal. A surface variability on timescales shorter than 20 minutes may remain in the light curves, however, and contribute to false positives. To optimise the detrending process, we followed Bruno et al. (2024) and tested polynomial degrees ranging from 2 to 10 in the Savgol filter. For each light curve, we selected the polynomial degree that minimised the Akaike Information Criterion (AIC) in the residuals, ensuring that more active stars, potentially exhibiting stronger modulation from surface activity, were detrended using higher-degree polynomials. Finally, we sigma-clipped outliers at each iteration at a 5σ threshold of the flux level in order to minimise the impact of large flares on the detrending process. We also explored other detrending methods, including non-windowed polynomial fitting, Gaussian processes with a Matérn kernel, and M-estimators using Tukey’s biweight function (Hippke et al. 2019). However, we found these approaches to be less effective at removing time-correlated variations. In particular, Gaussian processes performed the worst, likely due to the short duration of CHEOPS light curves preventing the model from converging. Ultimately, while allowing polynomial orders to vary slightly improved residual scatter, the most significant improvement came from using a windowed detrending approach. After visual inspection, we removed light curves exhibiting excessively high RMS values (>10 000 ppm), likely resulting from a failed photometric extraction and/or detrending, which accounted for ∼ 6% of the sample. Among the remaining light curves, 84% were observed at a cadence of 3 seconds, 14% at 5–7 seconds, and the remaining 2% at cadences longer than 12 seconds. Since the next steps of the analysis involve a denoising process, we refer to the detrended and normalised CHEOPS light curves as ‘original’ throughout the rest of the paper (as opposed to ‘denoised’).

2.3 Flare detection

CHEOPS orbits the Earth with a stable period of 99 minutes, but the available science time per orbit is reduced due to Earth occultations, the South Atlantic Anomaly, and other observational constraints. The uninterrupted science time of an observation is referred to as a visit, and the duration of a visit depends on the target’s position in the sky. Targets located near the ecliptic plane can benefit from longer uninterrupted observations, sometimes extending across multiple CHEOPS orbits, while targets near the ecliptic poles usually have shorter visits with frequent interruptions (Benz et al. 2021). To provide an overview of the actual on-target time composing our sample, Table 1 presents the number of stars per spectral subtype, the corresponding number of visits, and the total science time. A flare candidate cannot be considered as a candidate when it spans two separate visits because data points belonging to the flare profile would be missing from the resulting light curve. Recovering only part of a flare would lead to underestimating the flare energy and biasing the derived flare statistics. Our flare detection algorithm therefore only considered flares whose start and stop times occurred within the same visit. Figure 2 shows the distribution of the duration of the visits constituting our sample. Although most of the visits are as expected shorter than 99 minutes, we found a few visits with longer durations of up to 4.5 hours corresponding to targets easily observable by CHEOPS. Since M dwarf flares can have durations of up to ∼ 10 hours (Davenport 2016; Pietras et al. 2022), conducting flare detection with CHEOPS implies that the population of long-duration flares will not be recovered in this study.

The photometric signature of flares is typically characterised by a sharp rise (usually associated with bremsstrahlung radiation) followed by an exponential decay (corresponding to radiative cooling) (Kowalski et al. 2013; Davenport et al. 2014). An effective way to detect flares is by monitoring the increase in flux they produce in the stellar light curve. For this purpose, we employed AltaiPony⁴, a flare detection tool developed for Kepler, K2 and TESS observations (Davenport 2016; Ilin 2021). It implements a sigma-clipping algorithm based on the flare detection criteria of Chang et al. (2015): N₁, the minimum positive flux deviation from the median normalised by local scatter; N₂, a corresponding minimum deviation accounting for photometric errors; and N₃, the minimum number of consecutive data points meeting these conditions. For CHEOPS, we adopted N₁ = 3, N₂ = 2, and N₃ = 5, requiring a minimum of 5 consecutive points above the 3σ threshold. This setting reduces false positives but excludes flares shorter than 15 seconds.

Since the detection threshold of a sigma-clipping method for flare detection depends on the light curve noise properties, we intend to use a wavelet-based algorithm to denoise our data. This method is expected to lower the effective noise floor while preserving the amplitude and shape of flares, thereby increasing the S/N of low-energy flares and improving their recovery rate. To validate the advantages of this approach, we performed a comparative analysis with the PEAKUTILS⁵-based detection method used in Bruno et al. (2024), which combined a peak detection algorithm with a flare profile fitting. In our comparison, we applied both AltaiPony and PEAKUTILS, with and without wavelet-based denoising, using a consistent 3σ detection threshold. A summary of the results is presented in Table 2.

While PEAKUTILS detected a slightly higher number of flares overall, it also exhibited a higher false positive rate, especially in the denoised data. Flare injection-recovery tests showed that PEAKUTILS detects weak flares more readily, but also retrieved more spurious detections below ∼10²⁷ erg. In contrast, AltaiPony offers a more conservative detection by integrating a photometric uncertainty criterion. Its lower false positive rate and sharper detection threshold yield a more homogeneous flare sample. The application of wavelet denoising consistently improved detection rates for both algorithms. However, integrating it within AltaiPony provides a more robust sample for accurately assessing the impact of wavelet-based denoising on flare recovery rates.

Consequently, our specific workflow was to:

1. determine the noise profile of a ‘prototypical’ M dwarf light curve observed by CHEOPS.

2. generate a corresponding synthetic light curve in which to inject flares of different energies, durations, and amplitudes.

3. determine the denoising parameters that optimise the recovery of the injected flares.

4. apply the wavelet-based denoising algorithm with the optimal parameters on the original CHEOPS light curves.

5. conduct flare detection on the original and denoised CHEOPS light curves using AltaiPony.

A flow diagram of the entire procedure we followed is presented in Figure 3.

Fig. 2 Distribution of the visit duration within our sample.

Fig. 3 Flow diagram of the procedure we followed.

2.4 Estimating noise profiles

To determine the typical noise profile of a CHEOPS light curve, we computed the average noise profile for each spectral subtype by averaging the photometric fluctuations from the median value of each light curve in our sample, following the method of Jess et al. (2019); Dillon et al. (2020) and Grant et al. (2023). These profiles therefore quantify the averaged white noise and time-correlated scatter at timescales shorter than 20 minutes for each spectral subtype. We discarded the M6V subtype since only one star was included in this group (see Table 1). We therefore included targets of spectral type ranging from M0V to M5V. Figure 4 displays the averaged fluctuations for each spectral subtype and compare them against a normal distribution. The resulting noise distributions for each spectral subtype appear thinner and taller than a normal distribution, with more pronounced tails, indicating a leptokurtic nature. In statistics, the kurtosis is a measure of the shape of a probability distribution, specifically how the tails and peak compare to a normal distribution. A leptokurtic distribution (one with positive kurtosis) has a sharper peak and heavier tails than a normal distribution, indicating more extreme deviations from the mean. This is confirmed by the quantile– quantile (Q–Q) plots for each spectral subtype in the second row of Figure 4, which compare the observed data distribution to a theoretical normal distribution. If the data followed a normal distribution, the points would align along the diagonal. However, the deviation in the tails indicates an excess of outliers (see Figure 3 of Jess et al. 2019 for a collection of different distribution types). A leptokurtic noise distribution has several implications for the light curves and the flare detection. First, the presence of heavier tails means that extreme fluctuations occur more frequently than expected under a normal distribution. This could mimic or obscure low-amplitude flare events, complicating their detection. The shapes of the noise distributions indicate that the photometric noise in the light curves is not purely Gaussian and that the remaining time-correlated scatter contributes to a higher frequency of outliers. Consequently, the denoising process must accurately account for these leptokurtic features to efficiently distinguish between data points caused by flares and those arising from noise artefacts. Table 3 compiles the mean value, standard deviation, skewness, and Fisher kurtosis of the averaged noise distribution for spectral subtype. Understanding these specific noise characteristics allows us to refine the denoising and flare detection algorithms, reducing false positives and improving flare detection reliability.

Fig. 4 Averaged noise representation for each spectral type. In the top row, the dashed red curves correspond to a normal distribution, and the black histograms correspond to the noise distribution histograms in the CHEOPS light curves. In the bottom row, the dashed red lines correspond to the quantiles of a normal distribution, and the black curves correspond to the quantiles of the noise distribution in the CHEOPS light curves.

Fig. 5 Injected flare energy matrix as a function of the injected flare amplitude, FWHM, and spectral subtype.

2.5 Injecting synthetic flares

To create synthetic light curves for each spectral subtype, we first injected into empty light curves the noise distribution corresponding to each spectral subtype using the parameters compiled in Table 3. We then injected flares of known energies, amplitudes, and durations in the light curves corresponding to each spectral subtype. For this, we started by determining the stellar parameters of a ‘prototypical star’ from M0V to M5V. The effective temperatures were determined by calculating the average temperature of the stars in our sample within each spectral sub-type. The radius values for each spectral subtype were sourced from Pecaut & Mamajek (2013). The quiescent luminosities were calculated using the Stefan-Boltzmann law with the obtained effective temperatures and radii. We compile the average effective temperature, radius, and luminosity for each spectral subtype in Table 4.

Flares were simulated using the package Llamaradas Estelares⁶ (Tovar Mendoza et al. 2022). This package simulates the morphology of white-light flares with a template derived from a set of flares observed in GJ 1243 based on two parameters: relative amplitude and FWHM (Full Width at Half Maximum, serving as a proxy for flare duration). We could simulate flares with different energies by using a known range of parameters for the amplitude and the FWHM. Following Davenport (2016), flare energies were calculated by multiplying the Equivalent Duration (ED) by the stellar quiescent luminosity. Bruno et al. (2024) found that the sensitivity of CHEOPS could detect flares with energies down to ∼ 10²⁷ erg. We determined that relative flare amplitudes ranging between [10⁻⁴–10⁻²] and FWHM between [0–60] seconds could simulate flares of energies within [10²⁴–10³⁰] erg, depending on the quiescent luminosity of the star. Simulating flares within this energy range would enable us to test the impact of denoising on the recovery rate around this detection threshold. Figure 5 displays the injected flare energy as a function of the flare amplitude, FWHM, and spectral subtype of the star.

2.6 Determining the optimal denoising parameters

The Discrete Wavelet Transform (DWT) is commonly used to decompose and denoise time-series data. This is done by scaling and translating a mother wavelet across the data, resulting in a decomposition into different frequency bands through low-pass and high-pass filtering, followed by a downsampling. Stochastic noise, typically in high-frequency components, is thresholded and removed before reconstructing the data by using the inverse transform. Unlike other denoising methods based on the Fourier Transform, the DWT retains both temporal and frequency information, which is beneficial for recovering non-stationary signals (Mallat 2009). However, the Stationary Wavelet Transform (SWT) is often preferred for denoising as it avoids downsampling, maintaining the signal length at each decomposition level through an à trous algorithm. This approach offers a redundant, shift-invariant representation that enhances noise reduction accuracy and effectively preserves signal features, resulting in more consistent and reliable denoising results (Percival & Walden 2000; Fors 2006).

The threshold value in wavelet-based denoising is key to distinguish between noise and signal by determining which wavelet coefficients to suppress or retain. A well-chosen threshold effectively reduces noise while preserving essential signal features. However, an overly high threshold can erase important details, while a low threshold may fail to remove enough noise. Determining the threshold value typically involves statistical methods, such as the universal threshold proposed by Donoho & Johnstone (1994). Its logarithmic scaling with data length ensures effective noise filtering for larger time-series while preserving signal features (Abramovich & Benjamini 1995; Percival & Walden 2000). We adapted this universal threshold by replacing its dependence on the standard deviation of the signal by the median absolute deviation, which is more robust to outliers. This adjustment effectively balances noise reduction and signal preservation, especially in the presence of non-Gaussian noise (del Ser et al. 2018).

The decomposition level determines the number of times the data is iteratively decomposed into wavelet components, with higher levels providing finer resolution of the signal and stronger denoising. Each decomposition level corresponds to different frequency scales in the data, enabling multi-resolution analysis that typically separates noise from significant signal features. The decomposition level is therefore generally chosen based on the number of data points in the input data. An appropriate decomposition level allows us to effectively separate noise from the signal while preserving key signal characteristics. A level that is too low may fail to remove enough noise, while a level that is too high may suppress important signal details and introduce artefacts from the mother wavelet function (Yang et al. 2016). Since each decomposition level halves the data frequency, consequently capturing more detail at lower frequencies, wavelet decomposition follows a binary scaling (Mallat 1989). The maximum number of useful decomposition levels J therefore scales logarithmically with the number of data points N, such that 2^J ≤ N. This ensures effective noise isolation while preserving the signal integrity (Daubechies 1992). CHEOPS visits vary in the number of data points based on visit duration and observing cadence. The decomposition level is tailored to the number of data points in each visit, ensuring shorter visits are decomposed at a lower level than longer ones. This optimisation improves the denoising process and ensures effective use of the SWT across the different visit durations in our sample.

The choice of the mother wavelet is crucial in wavelet-based denoising, as it affects the accuracy of noise separation from the signal. Wavelets with short support, such as the Haar wavelet, are efficient at detecting sharp discontinuities but may introduce high-frequency noise, whereas wavelets with longer supports, like the Daubechies and Symlets, are better at capturing smooth variations but may blur sharp features (Haar 1910; Mallat 1989; Daubechies 1992). Vanishing moments refer to the ability of a wavelet to represent polynomials up to a certain degree. Wavelets with higher vanishing moments are more effective at filtering noise but are computationally demanding and may lead to over-fitting (Donoho & Johnstone 1994). For the SWT, discrete and biorthogonal wavelets are used to ensure accurate data reconstruction, non-redundancy, and minimal introduction of artefacts (Chui 1992). To identify the optimal parameters for denoising CHEOPS light curves, we followed a systematic approach, illustrated in Figure 6:

1. We created synthetic light curves for each spectral subtype and injected flares with energies from 10²⁴ to 10³⁰ erg, simulating various flare profiles. For each energy bin, spanning 10²⁵ erg, we simulated 1000 flares using the amplitude-FWHM parameters described in Figure 5. We separated flare peaks by 45 minutes to ensure that the flux had returned to a quiescent level before the beginning of the next flare.

2. We injected a noise profile corresponding to each spectral subtype according to the distribution parameters compiled in Table 3.

3. We conducted a SWT-based denoising using the previously mentioned threshold for each discrete and biorthogonal mother wavelet included in the package PyWavelets (Lee et al. 2019). This included wavelets from the Haar (haar), Daubechies (db), Symlets (sym), Coiflets (coif), and Biorthogonal (bior) families, with decomposition levels from 1 to 8. Detailed properties of all the considered wavelets are compiled in Table B.1.

4. We detected flares in the resulting light curves using the flare detection algorithm described in Section 2.3.

As an example of the output of this process, Figure 7 presents the improvement in flare recovery rate across all spectral subtypes after denoising with the Haar wavelet at a decomposition level of 2. The primary limitation for flare recovery is the amplitude, with a sharp decline in recovery rates for very low amplitude flares. This behaviour is expected, as even after denoising, low amplitude flares are unlikely to exceed the 3σ detection threshold above the noise floor and fulfil the N₁ criterion. In addition, recovery rates also drop for very low FWHM, likely due to the flares not satisfying the N₃ criterion due to insufficient consecutive data points above the detection threshold. This suggests that even though flares of different impulses can have the same energy (see Figure 5), high-impulse flares generally have a higher recovery rate compared to low-impulse ones, provided they have durations longer than N₃ data points. Overall, the denoising process improves recovery rates across all spectral subtypes, enabling the detection of flares with amplitudes smaller by approximately half an order of magnitude than those identified in the original light curves.

To compare the gain in recovery rate obtained after denoising with different sets of parameters, we averaged the recovery rate across spectral type for each wavelet and decomposition level. Figure 8 displays the average flare recovery rate obtained for all flares between 10²⁴ and 10³⁰ erg after denoising with different wavelets and decomposition levels. Recovery rates corresponding to each spectral subtype can be found in Figure C.1. Although a single wavelet did not produce a significant improvement in recovery rate compared to others, Coiflets and Daubechies wavelets consistently provided the highest recovery rate across several vanishing moments. Given their frequent use in analysing astrophysical datasets, we selected Daubechies wavelets for denoising the CHEOPS light curves (Ojeda González et al. 2014; Souza-Feliciano et al. 2018; Bolzan et al. 2020). Specifically, we selected the db18, which provides a balance between a high and low number of vanishing moments.

An important parameter of the SWT is the coefficient thresholding method. Hard thresholding sets coefficients below the threshold to zero, preserving signal features but potentially introducing discontinuities and artefacts in the reconstructed signal. Soft thresholding, on the other hand, not only sets coefficients below the threshold to zero but also shrinks the remaining coefficients towards zero, resulting in smoother data with fewer artefacts. However, it can slightly attenuate preserved features. Generally, hard thresholding is better for retaining sharp features, while soft thresholding improves signal smoothness and continuity. We compared the impact of both thresholding methods on the recovery rate and found only a marginal difference of less than 1% within each energy bin. We therefore chose to use a hard thresholding to denoise the CHEOPS light curves for several reasons. The signal of low-amplitude flares could be affected by the shrinkage introduced by soft thresholding, which would lead to the non-detection of the flares reaching just above the detection threshold. Additionally, the shrinkage could potentially lead to inaccurate energy estimations of the detected flares. Finally, hard thresholding is less computationally intensive than soft thresholding because it requires one operation less on each coefficient value. This makes it a more efficient method when it is applied to large datasets.

Our flare detection algorithm essentially functions as a binary classifier, distinguishing between ‘flaring’ and ‘non-flaring’ states. To evaluate its performance, we generated 1000 flare-free light curves for each 10²⁵ erg energy bin and applied the detection algorithm to measure its classification efficiency. Specifically, we constructed confusion matrices for each injected flare energy bin, defining four possible outcomes: true positive (TP), where a flare was injected and correctly retrieved; false positive (FP), where no flare was injected but one was erroneously retrieved; false negative (FN), where a flare was injected but not retrieved; and true negative (TN), where no flare was injected and none was retrieved. We then calculated the true positive rate (TPR), false positive rate (FPR), false negative rate (FNR), and true negative rate (TNR) in function of the injected flare energy.

In Figure 9, we present the evolution of these rates, averaged across all spectral types, for both the original and denoised light curves. We observed that the false positive rate (FPR) remained independent of the injected flare energy but instead depended on the noise level in the light curves. Notably, the FPR remained consistently low across the entire energy range: below 0.025 for non-denoised light curves and below 0.010 for denoised light curves.

The denoising process significantly enhanced the true positive rate (TPR), effectively shifting the detection threshold towards lower injected flare energies by nearly half an order of magnitude. To quantify these improvements, we identified four key thresholds:

• t₁: where the TPR for denoised light curves exceeds 10% (1.0 × 10²⁷ erg),

• t₂: where the TPR for original light curves exceeds 10% (3.5 × 10²⁷ erg),

• t₃: where the TPR for denoised light curves reaches 90% (1.1 × 10²⁸ erg),

• t₄: where the TPR for original light curves reaches 90% (4.2 × 10²⁸ erg).

Figure 10 illustrates the confusion matrices at these thresholds. Consistent with the trend shown in Figure 9, the denoising process significantly improves the TPR between t₁ and t₄, enabling a more complete detection of flares at lower energies compared to the original light curves.

Fig. 6 Step-by-step representation of the flare injection, noise injection, wavelet denoising, and flare recovery process. Only 20 flares are shown to better visualise the flux variations in the light curve. In this example, the wavelet-based denoising was conducted with the Haar wavelet and a decomposition level of 2. Only 16 out of 19 flares are recovered in this portion of the light curve. The recovered and lost flare candidates are highlighted with a green and red background, respectively.

Fig. 7 Flare recovery rate as a function of the injected flare amplitude, FWHM, and spectral subtype after denoising with the Haar mother wavelet and a decomposition level of 2. The top row corresponds to the recovery rates before denoising, and the bottom row shows the recovery rates after denoising.

Fig. 8 Flare recovery rate, averaged across all spectral subtypes, as a function of the selected mother wavelet and decomposition level.

Fig. 9 Evolution of confusion matrix metrics as a function of the injected flare energy. The TPR is shown in blue, the FPR in orange, the TNR in green, and the FNR in red. Metrics for the original light curves are represented by dashed lines, while those for the denoised light curves are shown as solid lines. The vertical black dotted lines indicate the four detection thresholds.

Fig. 10 Confusion matrices obtained from the flare injection and recovery process for the original (left) and denoised (right) light curves. The matrices are shown at the four detection thresholds, corresponding to TPR rises above 10% (t1 and t2) and reaching 90% (t3 and t4), as defined in the text.

2.7 Flare detection in denoised light curves

We then denoised all CHEOPS light curves (further referred to as ‘denoised’), separated by visit and performed the flare detection described in Section 2.3. Following the method of Raetz et al. (2020); Bruno et al. (2024) and Fortier et al. (2024), we determined the quiescent luminosity of each star in the sample by adopting the CHEOPS zero-point and effective wavelengths available on the filter profile service of the Spanish Virtual Observatory (Rodrigo et al. 2012) as follows: $L_0\text{\hspace{0.17em}=\hspace{0.17em}4}\pi d^2\times Q_{cheops}\text{(}{F}_{vega}\times {10}^{-0.4\times G}\text{)}$(1) where d is the distance of the star, Q_cheops is the quantum efficiency of the CHEOPS sensor integrated over the 330–1100 nm wavelength range, F_vega is the CHEOPS flux of Vega (CHEOPS zero-point), and G is the Gaia G-band magnitude of the star. We also derived several flare parameters, such as the flare duration and FWHM, flare impulse (amplitude/FWHM), and flare peak luminosity (amplitude × quiescent luminosity).

2.8 Complex flare breakdown

To break down complex flares displaying several visible peaks, we consider that such events can be described as a superposition of several simple flares, which can be modelled by the empirical flare template of Llamaradas Estelares. For a given complex flare event, the objective is therefore to automatically determine the number of flares and their parameters (namely t_peak, amplitude, and FWHM) that best describe it. For this, we developed an algorithm to iteratively fit up to 5 flares for every flaring event detected using AltaiPony. The flare region to fit for each event included 15 seconds before the start time and 45 seconds after the stop time of the flare to ensure the flux fully returned to the quiescent level within the fitting window. Following the procedure of Davenport et al. (2014), we imposed a first guess and boundaries for the fitted flare parameters. We seeded the parameters of the fitting by using the peak flux amplitude and time and 15% of the flare full duration for the FWHM. The relative amplitude of each component was required to be larger than 10⁻⁴, being the flare amplitude after which the recovery rate dropped to 0% across all spectral types (see Figure 7). We also required t_peak to occur within the boundaries of the flare window, and FWHM to be larger than the observing cadence of the light curve and smaller than 50% of the flare total duration. No priors on the relations between flare amplitude, FWHM, and t_peak were included. Following Bruno et al. (2024), we used the AIC to determine the best fit, which took the form $AI{C}_n=M\times ln(\frac{{\chi }^2}{M})+2k_n$(2) where M is the number of data points that fell within the flare window, and k_n is the number of degrees of freedom of the n-th model. This statistic determines the improvement of the fit (decreased χ²) while penalising the increasing number of free parameters used in the model. We chose to use the AIC instead of the BIC like Davenport et al. (2014) because the AIC penalises less the number of components. The AIC therefore tends to favour more complex models with a higher number of parameters, which would enable further breakdown. Finally, the best solution fit was selected as the n-th model with the smallest AIC to have decreased by at least 6 units from the previous (n − 1)th model. The choice of a 6-unit improvement threshold was determined in Bruno et al. (2024) by visual inspection of complex flare events to ensure the algorithm did not overfit them.

3 Results

We began our analysis by applying the flare detection algorithm to the detrended CHEOPS light curves, initially identifying 215 flares with energies ranging from 4.2 × 10²⁶ erg to 8.9 × 10³⁰ erg. To optimise the detection process, we employed a SWT-based denoising algorithm using the db18 wavelet base, as detailed in Section 2.6, to reduce noise within each CHEOPS visit. We then reran the flare detection algorithm on the denoised light curves and detected a total of 291 flares with energies rang- ing between 3.7 × 10²⁶ erg and 8.9 × 10³⁰ erg. After running the flare-breakdown algorithm, 101 individual components were recovered from 48 complex events in the original light curves, with energies ranging from 5.8 × 10²⁶ erg to 4.1 × 10³⁰ erg. 271 individual flare components were recovered in the denoised light curves from 123 complex events, with energies ranging from 1.9 × 10²⁶ erg to 4.1 × 10³⁰ erg. The denoising process resulted in a significant increase in the overall flare recovery of ∼ 35% and improved the recovery of individual flare components by ∼ 64%. For complex flares, the number of recovered components increased by ∼ 168%. Eight events (less than 2% of the denoised flare sample) had recovered energies lower than 10²⁷ erg, therefore belonging to the upper end of the micro-flare energy range.

In Figure 11, we present examples of the improved flare recovery enabled by wavelet-based denoising. As shown in panel b, the denoising allowed the recovery of a low-amplitude flare from EV Lac by smoothing out the local scatter caused by the light curve noise, which had previously prevented the flare from meeting the N₃ criterion. Flares slightly below the detection limit therefore become detectable as they fulfil the N_1,2,3 criteria. As shown in panel d, the denoising also allowed for the identification of an additional flare component in a complex flare from AD Leo, which had been merged into a single one due to the light curve noise. This increased sensitivity might allow us to recover a larger population of complex flares, with much lower amplitudes and durations than those detected before applying the denoising.

We also observed several events with well-known flare morphologies. Panel a of Figure 11 shows a pre-flare dip from AU Mic, with an amplitude comparable to the one of the subsequent flare, while panel c presents a QPP candidate with a sub-minute period from V1054 Oph. Such pre-flare dips and QPPs in flaring light curves have been extensively discussed in Bruno et al. (2024). We present more pre-flare dip detections and QPP candidates from GJ 317, GJ 3323, AU Mic, and AD Leo in Figure D.1. Furthermore, we detected many flare candidates with the typical ‘peak-bump’ morphology, which were broken down by our algorithm into two distinct components. While some of these events are explicable as a random superposition of two sympathetic flares, Tovmassian et al. (2003) claims that most peak-bump events are caused by a single flare described by a two-phase underlying emission mechanism, as energy is re-radiated by the stellar photosphere after the peak phase. The complete profile should only be observable when the emission site is close to the centre of the visible stellar disk, however. We also found a few occurrences where the ‘bump’ occurs before the ‘peak’, however, which seems to contradict this hypothesis. Finally, a few cases of ‘flat-top’ structures were also found, although this could be due to a typical ‘peak-bump’ shape with a lower contrast between the ‘peak’ and the ‘bump’ phases. Examples of such flare morphologies from Gl 841 A, V1054 Oph, and Ross 733 are provided in Figure D.2.

Figure 12 displays the number of peaks observed per flaring event. We find that ∼42% of flares have complex structures composed of multiple peaks, in agreement with previous studies on M stars (Davenport et al. 2014; Howard & MacGregor 2022; Bruno et al. 2024). In Figure 13, we show the number of components detected in flares that were observed both before and after denoising. In most cases, denoising allows us to identify additional components and enables us to break down complex flare structures in more detail. We compare the percentage of peaks per flaring event with those obtained by Bruno et al. (2024), who conducted a similar analysis on CHEOPS light curves from the same programmes. We find that our method retrieved less 3-to 5-peak flares than Bruno et al. (2024), which we attribute to two factors. First, a minor difference in determining the ‘best-fit’ within the breakdown algorithm: we required each (n − 1)th model to improve the AIC by at least 6 units before exploring the n-th model, similar to the approach used by Davenport et al. (2014) to prevent overfitting. Second, as shown in Figure 14, denoising facilitates the detection of new flares that were not identified in the original sample. These newly detected flares are generally of low amplitude and short duration, making them more challenging to decompose into multiple components, and contributing to a higher proportion of single- and two-peak flares in the denoised sample.

In Figure 15, we present a corner plot of the parameters of the recovered flares before and after denoising. In both cases, the energy of recovered flares spans from ∼10²⁶ erg to ∼10³¹ erg. Most flares exhibit relative amplitudes between 10⁻³ and 10⁻¹. Flare durations range between 10 seconds and 20 minutes. Most flares display impulses between 10⁻⁵ and 10⁻³ s⁻¹. The distributions of flare energy, duration, amplitude, and impulse are notably similar between flares detected in the original and denoised CHEOPS light curves, suggesting that the population of recovered flares remains consistent across both datasets. We note a small secondary peak in the flare peak luminosity distribution after denoising, however, that resembles a bimodal pattern that was absent before denoising. We attribute this to the small sample size.

The denoising process led to a higher recovery rate across the entire energy range, with the most significant improvements observed for flares weaker than 10²⁹ erg. This is consistent with the results of the injection-recovery test, which suggested that the flare recovery rate started to drop at this energy (see Figure 9). As can be seen from the amplitude and duration histograms, the denoising allows us to recover an increased proportion of low amplitude and short duration flares, while those with relative amplitudes higher than 10⁻² and durations longer than ∼4 minutes are recovered equally well. This confirms that observing cadence and photometric precision is the main bottleneck for low-energy flare detection. By optimising the S/N of weak flares, the denoising improves one of these factors and facilitates the recovery of a greater proportion of low-amplitude flares. Additionally, denoising appears to increase the flare recovery across most of the impulse range. Contrarily to Hawley et al. (2014), we did not find that individual components of complex flares have longer durations than simple flares.

We provide a tentative linear fit to evaluate correlations between components. A power-law relation is visible between the flare duration and recovered energy, with longer flares generally releasing more energy. Similarly, the recovered energy logically scales with amplitude. On the other hand, the recovered energy does not seem to be correlated to the flare impulse, implying that the flare detection and breakdown algorithm efficiently recovers flares of various impulsiveness. A slight trend is visible between the flare amplitude and duration, where the flare duration appears to scale inversely with the amplitude. This might be due to a detection bias and implies that within this energy range, the detection of flares of low amplitude and short durations (i.e. the flares at the lower end of the energy range) is incomplete. This is expected considering the results of the injection-recovery process (see again Figure 9). Finally, as expected, flare peak luminosities scale with flare amplitudes, and FWHMs scale with durations, as these properties act as proxies for each other.

Considering the denoised light curves, flares were detected on 62 out of the 110 stars in our sample. In Figure 16, we compare the stellar parameters of the flaring stars to the ones of the entire stellar sample. No clear trends are evident in G_mag, T_eff, radius, v sin i, distance or spectral type. Previous studies have attributed a higher proportion of flaring stars to fast-rotating stars, which are generally younger and possess more active magnetic dynamos than their slower-rotating counterparts (Howard et al. 2020). This trend has been observed to decline significantly in ultra-fast rotators (P_rot < 0.2 days), however, which might be due to the high rotation speed that constrains the magnetic field lines in the stellar chromosphere, thereby inhibiting magnetic reconnection (Günther et al. 2020; Doyle et al. 2022; Ramsay et al. 2022). Previous large-scale studies have also attributed a higher rate of flaring stars on late Ms compared to their earlier counterparts (Günther et al. 2020; Pietras et al. 2022). Our results do not provide sufficient evidence to support any of these hypotheses, likely due to the limited sample size.

Figure 17 shows statistics on the flares separated by spectral subtype. The majority of flares were detected on M1V stars, likely due to the substantial observing time dedicated to AU Mic in our sample. M5V stars produced the fewest flares, likely due to the small number of stars in this group. The mean flare energy appears consistent across spectral subtypes but slightly decreases towards late-type stars. This effect is expected as the average quiescent luminosity decreases with spectral type (a similar trend is visible in Figure 5). These mean flare energies remain lower by more than an order of magnitude than the mean flare energy recovered by Bruno et al. (2024) using combined TESS and CHEOPS data, however, and they are lower by several orders of magnitude than the flare samples observed by TESS (Günther et al. 2020; Pietras et al. 2022), Kepler (Davenport 2016; Yang et al. 2017), and Evryscope (Howard et al. 2019). When comparing the mean equivalent duration, we find that it tends to increase with later spectral subtypes, suggesting that late M dwarfs produce flares that are more energetic relative to their quiescent luminosities. This trend could be attributed to enhanced magnetic activity in fully convective stars. Flare rates appear to be consistent across spectral subtypes, except for M3V stars, which exhibit a mean of 30 events per day. Most distributions are heavily skewed by a few frequently and intensely flaring stars, however, as indicated by mean values exceeding the third quartile. This also highlights the limited statistical confidence due to the small sample size.

To better compare the energy range of our flare sample with that of Bruno et al. (2024), we display in Figure 18 the histogram of the energy distribution of our original and denoised flares, separated into simple flares and individual components of complex flares. First, we find that denoising preserves the number of simple flares while significantly increasing the number of complex flares, suggesting that it helps to break down initially simple flares into several components. Likewise, the number of high-energy complex flares remains similar, indicating that denoising helps to further resolve them into finer components. When comparing the probability density of the flare energies with those in Figure 13 of Bruno et al. (2024), the mean flare energy we recovered is lower by more than an order of magnitude for simple and complex flares, which is likely due to the contribution of flares detected by TESS composing the majority of their sample.

To explore potential relations between flare parameters and the host star properties, we analysed correlations between flare parameters and those of the corresponding star. We show some of these correlations in Figure 19. A slight exponential trend is observed between the targets’ G-band magnitude and the relative amplitude of the recovered flares. This trend may be attributed to low-amplitude flares being more easily detectable on brighter stars and suggests that a portion of such flares remain undetected in the noise floor of the light curves of dimmer stars. Jackman et al. (2021) and Pietras et al. (2022) found a correlation between the flare amplitude and the effective temperature of the star, with flares detected on cool stars having, on average, higher amplitudes because of the greater contrast between the flaring and quiescent emission. We find a similar trend, although less pronounced, for our simple and complex flares. Moreover, a power-law trend is visible between the flare duration and the stellar rotational velocity, indicating that longer flare durations are usually associated with faster-rotating stars. This may be related to fast-rotating stars often having stronger magnetic fields generated by the dynamo effect, which can lead to more intense reconnection events. Finally, an exponential trend is observed between the flare energy and the distance of the star. This is likely due to a detection bias, where low-energy flares become increasingly difficult to detect as the distance to the star increases. We present all remaining correlations, classified as simple versus complex flares, in Figure E.1, and provide the same correlations categorised as flares detected in the original versus denoised light curves in Figure E.2.

The relation between the flare energy and occurrence frequency is generally assumed to be characterised by a power-law described as dN/dE ∼ E^−α, where α is the power-law index governing the frequency dN/dE of flaring events with energy E. To compare the impact of denoising on our recovered flare samples, we constructed Flare Frequency Distributions (FFDs) for our flare catalogues recovered in original and denoised light curves. The FFDs represent the cumulative rate of flares per day, indicating the frequency of detecting a flare of a given energy or higher. This approach excludes inactive stars from the analysis. In Figure 20, we compare the FFDs of flares recovered in original and denoised light curves. For the denoised sample, we also compare FFDs of simple flares and individual components of complex flares. Following Yang et al. (2023), we also separate flares from fully and partially convective stars by setting the separation threshold at M3V. We observe that simple flares are more frequent before denoising, as are individual components of complex flares with energies above 10²⁹ erg. We attribute this to the denoising facilitating the identification of additional complex flare components, as shown in Figure 13 and in the increased percentages in Figures 12 and 18. Thus, some simple flares are broken down into multiple components, and high-energy components are split into smaller ones after denoising, explaining the increase in observed frequency. When comparing the full sample of original versus denoised flares, we see that the frequency of recovered flares is higher after denoising across the entire energy range. Even so, the denoising process only marginally extends the lower limit of the detected energy range. Similarly, when comparing flares from fully and partially convective stars, the denoising process generally allows us to recover an increased frequency of flares across the entire energy range. This is not the case for flares from fully convective stars above 10²⁹ erg, which could be caused by the smaller sample size compared to partially convective stars (see the comparison of the number of flares per spectral subtype in the first panel of Figure 17). Finally, in the denoised light curves, individual components of complex flares appear more frequently at the same energy compared to simple flares. This suggests that highly energetic flares are more likely the result of multiple weaker, sympathetic flares originating from neighbouring active regions, rather than single, powerful events. Consequently, within this energy range, flares appear to occur more frequently in groups rather than as isolated occurrences.

For each obtained FFD, we determined the corresponding power-law index α by using an MCMC fit based on the Bayesian flare frequency predictor outlined in Wheatland (2004) and Ilin et al. (2021). The fit range began at 4.4 × 10²⁸ erg for the original light curves and 1.1 × 10²⁸ erg for the denoised light curves, corresponding to the energy from which 90% of flares are expected to be recovered (see t₃ and t₄ from Figure 9). We performed a Kolmogorov-Smirnov test on each fit to assess whether the power-law assumption should be rejected, using a significance level of 0.05 (Weidner et al. 2009). In all cases, the power-law assumption was accepted with p-values below 0.05. We compile the obtained α values in Table 5. The obtained values correspond to the higher end of the general α range reported in the literature (Günther et al. 2020; Pietras et al. 2022; Yang et al. 2023). Specifically, we find α = 1.99 ± 0.10 for the entire flare population recovered in the denoised light curves. This value is close to the crucial threshold of α_H > 2, above which smaller flares are thought to play a dominant role in coronal heating (Hudson 1991). Caution is advised when interpreting the obtained α indexes, however, because the portion of the energy range used for the fit is limited and only represents a small subset of events. For the flares detected in the original light curves, only about 34% had energies above t₄ and could be included in the power-law fit. Similarly, for the denoised light curves, just 33% of flares had energies above t₃ and were used in the fit. This highlights that extending the detection limit would allow to improve the detection completeness and obtain a more accurate extension of the FFDs towards lower energies.

A power-law distribution is generally assumed when fitting FFDs in order to determine an α index. A power-law implies that flares are governed by a scale-free process, where flares arise from self-organised criticality. This is the current understanding of magnetic reconnection events, where magnetic loops can accumulate energy and release it in a scale-invariant manner, with a size distribution following a power-law. Verbeeck et al. (2019) and Sakurai (2022) suggested, however, that a preliminary step should be to verify that the FFD indeed follows a power-law and could not be better described by an alternative distribution. A log-normal distribution would indicate that high-energy flares result from several smaller independent events combining multiplicatively, rather than from a critical build-up and sudden release of energy. This could happen if flares resulted from many independent magnetic reconnection sites or gradual energy accumulation processes. If flares followed a log-normal distribution, it would suggest that the processes driving flare energy release operate within a characteristic energy range, with fewer instances of extreme events, which could explain the ‘trun-cated power-law’ distribution observed by Howard et al. (2019) at high flare energies.

To test this, we followed the procedure of Verbeeck et al. (2019) and Bruno et al. (2024) and compared our obtained FFDs to both a power-law and a log-normal distribution. Figure 21 displays the Complementary Cumulative Distribution Functions (CCDFs) of the flare energies recovered in our denoised light curves, and compares the fitted power-law obtained in Figure 20 to a log-normal distribution. Again, the fitting range started at the energy of 90% recovery (t₃ = 1.1 × 10²⁸ erg). For each denoised category, we identified the best-fitting distribution by calculating the log-likelihood ratio R between the power-law and log-normal fits. A positive R indicates a preference for the power-law, while a negative R favours the log-normal. Our results show that a log-normal distribution better fits the data in all cases, although the p-value of complex flares and flares from fully convective stars is found above 0.05, showing a lack of statistical significance for these subsets. This further emphasises that extending flare observations to lower energies could help clarify the best-fitting distribution for FFDs, which may have important implications for understanding the mechanisms driving stellar flares.

Fig. 11 Comparison of the flare recovery before (left) and after (right) denoising. In each panel, the original light curve (black data points) is compared to the denoised light curve (red data points). The vertical bars indicate the photometric uncertainty associated with each flux measurement, as estimated during the photometry extraction. The star on which each event was detected is mentioned in the corresponding panel title. The green shaded areas correspond to the data points flagged as belonging to a flare candidate by the flare detection algorithm, and the red shaded areas correspond to a non-detection in the same portion of the light curve. The solid lines represent the fitted flare components. Below each light curve, we show the residuals after subtracting the best-fit flare model. In the top right corner of each panel, the ${\tilde{\chi }}^2$ calculated on the portion of the light curve we show is given as indication. Each panel group corresponds to a specific case: (a) Simple flare detected in the original and denoised light curves. (b) Simple flare detected only after denoising. (c) Complex flare with all components recovered in the original and denoised light curves. (d) Complex flare with an additional component recovered after denoising.

Fig. 12 Percentages of peaks observed per flaring event. The black and red bars correspond to flares recovered in original and denoised light curves, respectively. We compare our obtained percentages to the ones of Bruno et al. (2024), displayed as dashed blue bars.

Fig. 13 Comparison of the number of peaks detected in each flare before and after denoising. Each point corresponds to a single flare. The red dashed line indicates the 1:1 relation, where the number of peaks remains unchanged. Points above the line represent flares that gained components after denoising, while those below lost components.

Fig. 14 Amplitude vs duration of newly detected flares identified only after denoising. Each point represents a flare, with colours indicating the number of components retrieved from the flare breakdown.

Fig. 15 Corner plot of the different flare properties, namely relative amplitude, duration, measured energy, FWHM, impulse, and peak luminosity. The black data points represent metrics of individual flare components detected before denoising, while red data points correspond to those detected after denoising. A linear fit is provided for reference in each distribution. The discrete bins for flare duration and FWHM are due to the observational cadence equal to 3 seconds.

Fig. 16 Comparison of Gmag, Tef f, radius, log vsini, distance and spectral type for flaring stars (red bars) against the entire stellar sample (white bars).

Fig. 17 Flare statistics by spectral subtype for the stars in our sample. The top panel shows the mean number of flares per star, the middle panel displays the mean flare energy, and the bottom panel presents the mean flare equivalent duration. Vertical bars represent the interquartile range (1st to 3rd quartiles), while horizontal bars indicate the median values. In the M5V bin of the bottom panel, only one star is included, resulting in equal values for the 1st quartile, median, and 3rd quartile.

Fig. 18 Histograms of the number of flares observed per energy bin. The grey shaded histogram represents simple flares in denoised light curves, and the red shaded histogram represents individual components of complex flares. The dashed black and red histograms correspond to the same in the original light curves.

Fig. 19 Correlations between flare parameters and stellar parameters. The blue and orange data points correspond to simple flares and individual components of complex flares, respectively. Each panel includes a linear fit for reference.

Fig. 20 Flare frequency distributions for the flares recovered in our original (black) and denoised (red) samples. We separate FFDs into simple flares (top left), individual components of complex flares (bottom left), flares recovered on partially convective stars (top middle), flares recovered on fully convective stars (bottom middle), and the entire flare samples (top right). We also provide a comparison between the simple flares (blue) and individual components of complex flares (orange) recovered in the denoised light curves (bottom right). Solid lines indicate linear fits to the double-logarithmic FFD, extrapolating into regimes not directly observed.

4 Discussion

The use of CHEOPS allowed us to study an unexplored region of the white-light flare energy range in main sequence M dwarf stars. Moreover, we employed a tailored denoising algorithm based on the SWT to optimise the detection rate of low-energy flares. Several potential improvements could extend this work, however.

The primary limitation encountered in the flare observation stems from the constraints of the detection algorithm employed. The algorithm we used was initially developed for standard flare detection and optimised for TESS and Kepler light curves. Since low-energy flares are generally detected just above the 3σ threshold, many flares with variations below this level are missed. Lowering the N₁ criterion could help capture more low-energy flares, but would significantly increase the false positive rate, thereby affecting the overall flare statistics. We therefore highlight the limitations of a classical sigma-clipping detection method when applied to low-energy flares. Emerging detection algorithms, such as those based on hidden Markov models, use non-stationary time-series and have shown improved efficiency in detecting low-energy flares and accurately measuring their energy (Esquivel et al. 2025; Zimmerman et al. 2024). Flare detection algorithms employing deep neural networks and random forest classifiers have also demonstrated improved recovery rates and increased sensitivity to low-energy events in Kepler and TESS data, outperforming classical detection techniques (Vida et al. 2021; Lin et al. 2024). Integrating such algorithms could improve the sensitivity to low-energy flares and increase the confidence in the derived flare statistics.

Another limitation of the flare detection algorithm is the required number of consecutive data points for identifying a flare candidate. To minimise false positives, we set the minimum consecutive data points for detection to N₃ = 5. Given that the observing cadence for most light curves we studied is 3 seconds (see Figure 2), this results in a minimum detectable flare duration of 15 seconds. Several studies using observing cadence of 1 second or faster have shown, however, that flares can have durations as brief as a few seconds (Schmitt et al. 2016; Aizawa et al. 2022). As illustrated in Figure 15, flare energy typically scales with duration, which is supported by extensive statistical analyses conducted by TESS and Kepler (Davenport 2016; Howard et al. 2019; Pietras et al. 2022). Enhancing the observing cadence while conserving the same photometric precision could enable the detection of these shorter-duration, less energetic flares and potentially extend the flare detection toward the lower end of the energy range. CHEOPS can reach an exposure time as fast as 0.001 seconds for very bright stars (V<6) (Broeg et al. 2013; Benz et al. 2021). Conducting a dedicated observing programme targeting such bright stars could provide high-cadence light curves of active stars, thus decreasing the minimum detectable flare duration.

An increased cadence would capture more data points, allowing for more precise energy estimates (Clarke et al. 2024). Denoising low-cadence time-series can cause flare peaks to be treated as outliers and pulled toward the quiescent level, particularly for very short flares with limited data points. This effect leads to a higher proportion of lost energy and increased relative error in energy calculations, as the flare equivalent duration is determined from the area under the curve. Accelerating the observing cadence could help alleviate this issue by providing more data points and minimising the impact of denoising on flare peak values.

Moreover, the flare profile used for the injection and recovery process in Section 2.5 is based on GJ 1243 flare data observed by Kepler (Davenport et al. 2014; Tovar Mendoza et al. 2022). Seli et al. (2025) identified a correlation between the average flare shape and spectral type, however, which this profile does not account for. Additionally, Pietras et al. (2022) demonstrated that this model overestimates the decay phase of short-duration flares, while Dillon et al. (2020) and Grant et al. (2023) found that low-energy flares tend to exhibit more symmetric profiles with balanced rise and decay times. This suggests that the profile used here may not be ideal for modelling low-energy flares, potentially leading to the selection of suboptimal denoising parameters. As a result, the choice of mother wavelet might be biased toward identifying longer, higher-energy flares while being less effective for shorter, low-energy ones. A more representative flare profile in the injection and recovery process could improve the determination of optimal denoising parameters, requiring a detailed analysis of how flare morphology varies with energy.

Furthermore, the injection and recovery process in Section 2.5 were limited to simple flares. At high cadence, however, flares often appear as stacked events, adding complexity to their light curve profiles (Davenport et al. 2014; Bruno et al. 2024). Even low-energy flares can exhibit diverse morphologies, showing the need of separating individual events through improved photometric sensitivity and cadence. A more refined injection-recovery process incorporating both simple and complex flares would ensure that the chosen denoising parameters effectively capture both types. Extending the empirical flare model of Tovar Mendoza et al. (2022) to include well-known morphologies (e.g., ‘peak-bump,’ ‘flat-top’) could also help misclassification of complex flare events.

Additionally, recent advancements in wavelet transform applications for time-series denoising offer promising avenues for improvement. Polyharmonic wavelets, which use functions able to represent multiple harmonics, provide better frequency localisation and improved signal decomposition across multiple scales, reducing artefacts in the reconstructed signal (Kounchev et al. 2010). Although not in an astronomy context, Ji et al. (2016) and Alvarado et al. (2023) have also explored hybrid approaches that combine different wavelets or transformation methods to optimise denoising for varying signal complexities. Investigating how these methods impact flare recovery rates compared to standard SWT denoising could provide valuable insights and help refine denoising parameters for optimal performance.

Most importantly, the main limitation of this study is the relatively small sample size, with only 62 flaring stars identified and less than a thousand flaring events catalogued. This constraint arises from the limited data gathered by CHEOPS due to its on-target observing mode. While precise for focused targets, this mode is less effective at quickly monitoring a large number of targets than survey-operating missions like TESS and Kepler. As a result, flare statistics derived from CHEOPS data carry more uncertainty than those from TESS and Kepler, which observed thousands of stars (Davenport 2016; Pietras et al. 2022). FFDs are generally constructed using hundreds or thousands of flaring events. In our study, however, the average number of flares per star was around 10 for each spectral sub-type. This small sample size likely impacts the reliability of the derived FFDs, leading to substantial uncertainties in the calculated α values and possibly explaining the better performance of log-normal fits over power-law fits. This limitation emphasises that, while CHEOPS provides the photometric precision and observing cadence needed for low-energy flare observations, its on-target mode makes it less suited for large-scale flare studies. Future missions combining high photometric precision, fast cadence, and survey-mode capabilities will be crucial for advancing our understanding of low-energy flare statistics.

PLATO (PLAnetary Transits and Oscillations of stars), scheduled for launch at the end of 2026, stands out as the most promising candidate for this matter. Although its primary objective is to observe transits of Earth-like exoplanets around Sun-like stars, it will also collect photometric data on active late-type stars as part of its stellar variability programme. PLATO is anticipated to achieve a photometric precision comparable to CHEOPS, with a fast observing cadence of up to 2 seconds in its high-cadence mode. It will monitor two 49°× 49° regions of the sky for two continuous years each (Rauer et al. 2024). The combination of observational capabilities and a survey-mode strategy can significantly increase the detection of low-energy flares and improve the precision of their statistics for late-type stars. Additionally, for the brightest targets, PLATO will provide simultaneous photometry in two bandpasses (red and blue), allowing for direct colour measurements of flares. This will enable more accurate blackbody temperature estimates and extend energy modelling towards low-energy flares, reducing uncertainties in flare energetics. Moreover, the long-term uninterrupted observations from PLATO could help track flare activity cycles and investigate correlations with stellar rotation and magnetic activity over multi-year timescales.

Furthermore, the recently announced TESS Extended Mission 3 presents great capabilities for low-energy flare studies. TESS has been of prime importance in gathering flare statistics in recent years and is set to extend its observing cadence to 2 seconds from 2025 to 2028. While its photometric precision remains lower than those of CHEOPS and PLATO, the enhanced cadence will improve the detection of lower-energy flares compared to previous TESS surveys. This would help bridge the gap between large-scale studies of regular flares and studies of low-energy flares, which remain limited by small sample sizes (Günther et al. 2020; Pietras et al. 2022; Feinstein et al. 2024).

Finally, the type of denoising algorithm we presented might improve the recovery of low-energy flares in JWST exoplanet spectra and facilitate the understanding of the impact of flares on exoplanetary atmospheric retrievals (Howard et al. 2023). During occultations, certain atmospheric elements (e.g. CO₂) can produce spectral profiles similar to flare activity, potentially biasing the measurements of emission features of planetary atmospheres (Lustig-Yaeger et al. 2019). Likewise, flares themselves can mimic the appearance of emission lines from various atmospheric elements, complicating the accurate identification of chemical signatures (Lim et al. 2023). Flares are known to produce strong chromospheric emission in key spectral lines such as H_α, CaII H&K, and HeI 10 830 Å, which are also used as tracers for planetary atmospheres (Jensen et al. 2012; Yan et al. 2019; Sanz-Forcada et al. 2025). This can lead to false detections of atmospheric species in transmission spectra if the flare signal is not properly accounted for (Konings et al. 2022). Similarly, in emission spectra, flares can alter the measured stellar flux during a planetary occultation, introducing variability in the observed planetary signal. Additionally, the intense UV and X-ray radiation from flares can heat the upper layers of planetary atmosphere, enhancing atmospheric escape and affecting long-term emission measurements. This effect is particularly important for planets orbiting M dwarfs due to their enhanced activity. Since both transmission and emission spectra can be contaminated by flare-induced spectral features, it is essential to develop techniques to separate true planetary atmospheric signals from stellar activity (Thompson et al. 2024). Applying a similar denoising algorithm could help differentiate between genuine atmospheric features and transient flare signals, thereby reducing the reduce the risk of misidentifications.

Fig. 21 Complementary cumulative distribution functions for recovered flare energies in denoised light curves. We separate flares into simple events (blue) and individual components of complex flares (orange), as well as flares from partially (green) and fully convective stars (purple). The full flare sample is displayed in black. Observed distributions, power-law fits, and log-normal fits are displayed with solid lines, dashed lines, and dotted lines, respectively. The grey dashed vertical lines indicate the energy of 90% recovery (t3) and the beginning of the energy range used for fitting. In each panel, the legend indicates the log-likelihood ratio R between the power-law and log-normal fits and its associated p-value.

5 Conclusion

We used CHEOPS to study white-light flares in an unexplored region of the energy range. We employed a denoising algorithm based on the SWT to enhance the detection of low-energy flares in the light curves of a sample of M dwarfs. We optimised the denoising parameters by performing a flare injection and recovery in synthetic light curves that replicated the noise profile of the CHEOPS instrument. The CHEOPS light curves were then denoised using a db18 wavelet and a decomposition level tailored to the number of data points in each visit. Each flare was analysed by a breakdown algorithm to identify substructures. We finally compared the flare populations detected in the original and denoised light curves.

After denoising, our sample was composed of 291 flaring events, with energies ranging between 3.7 × 10²⁶ erg and 8.9 × 10³⁰ erg. We identified a portion of ∼ 42% of flares that had complex structures. 123 complex flares were divided into individual components by using the flare-breakdown algorithm. The denoising process significantly improved the flare detection. It increased the overall recovery rate by ∼ 35% and the identification of individual flare components by ∼ 64%. An improvement in recovery was observed across the entire energy range, but the highest gain was obtained for flares with energies lower than 10²⁹ erg. We identified several trends between flare parameters that are likely due to a detection bias, which suggests that the flare population we studied is still not entirely recovered at low energies.

The limited sample size prevented us from establishing any clear relation between the flaring frequency and the spectral type, or between the proportion of flaring stars and the spectral type. We observed a trend, however, in which the mean flare energy decreased from early-to late-M stars, while the mean equivalent duration increased. We constructed FFDs and applied power-law fits within the full recovery regime to determine the α indexes. For the complete flare sample, we obtained α = 1.99 ± 0.10, which is at the higher end of the values typically reported in the literature. We found no statistically significant differences in the α indexes between the distributions of simple and complex flares or between flares from partially and fully convective stars.

We also explored an alternative scenario to describe the observed flare distribution by fitting a log-normal distribution and comparing it to the traditional power-law. The log-likelihood-ratio tests consistently favoured the log-normal fit over the power-law, although statistical confidence was not achieved for all subsets. If flares were indeed better represented by a log-normal distribution, it might indicate that high-energy flares arise from the superposition of multiple weaker sympathetic flares that erupt from nearby active regions. This shift could have strong implications for understanding the flare-formation mechanisms, for the impact of flares on exoplanetary atmospheres, and for assessing long-term stellar activity.

Finally, we reviewed the main limitations of this study. First, the small sample size introduces considerable uncertainty in the derived flare statistics, which leads to a low confidence in the α indexes extracted from the FFDs. Second, the flare-detection algorithm and the model we used to inject flares were both designed for regular flare studies and may bias the denoising parameter selection. This might favour the recovery of regular over low-energy flares. Third, although CHEOPS outperforms most other instruments in photometric precision and observing cadence, which allows it to detect flares with energies corresponding to the upper end of micro-flares, it still cannot capture all low-energy flares. This creates a detection bias in the flare catalogue because the detections are incomplete. We proposed several methods for mitigating these limitations and commented that upcoming missions, such as PLATO and the third mission extension of TESS, could significantly enhance low-energy flare statistics.

Acknowledgements

The authors thank the anonymous reviewer for their insightful comments, which greatly contributed to improving the manuscript. This work was partially supported by the Spanish programme MICIN/AEI/10.13039/501100011033 and by the “European Union Next Generation EU/PRTR” through grant PCI2022-135049-2 and by the “ERDF A way of making Europe” by the European Union through grants PID2021-125627OB-C31 and PID2022-136828NB-C41, by the programme Unidad de Excelencia María de Maeztu CEX2019-000918-M to the ICCUB and CEX2020-001058-M to the ICE-CSIC, and by the Generalitat de Catalunya/CERCA programme. I.R. acknowledges financial support from the European Research Council (ERC) under the European Union’s Horizon Europe programme (ERC Advanced Grant SPOTLESS; no. 101140786). This research used the new computing cluster Nyx of the ICCUB. CHEOPS is an ESA mission in partnership with Switzerland with important contributions to the payload and the ground segment from Austria, Belgium, France, Germany, Hungary, Italy, Portugal, Spain, Sweden, and the United Kingdom. CHEOPS public data analysed in this article are available in the CHEOPS mission archive. We also used the Data Analysis Center for Exoplanets (DACE), dedicated to extrasolar planets data visualisation, exchange, and analysis. DACE is a platform of the Swiss National Centre of Competence in Research (NCCR) PlanetS, federating the Swiss expertise in Exoplanet research. Additional resources include the SIMBAD database and VizieR catalog access tool from CDS, Strasbourg Astronomical Observatory, France (Wenger et al. 2000; Ochsenbein et al. 2000). Software: Python (Van Rossum & Drake 2009), AltaiPony (Ilin 2021), AstroPy (Astropy Collaboration 2022), AstroQuery (Ginsburg et al. 2019), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), Pandas (Wes McKinney 2010), PIPE (Brandeker et al. 2024), Powerlaw (Alstott et al. 2014), PyCHEOPS (Maxted et al. 2023), PyWavelets (Lee et al. 2019), and SciPy (Virtanen et al. 2020).

Appendix A Target list

Appendix B Wavelet list

Appendix C Flare recovery heatmaps

Fig. C.1 Flare recovery rate as a function of the selected mother wavelet and decomposition level, separated per spectral subtype.

Appendix D Additional flare detections

Fig. D.1 Additional flare detections similar to Figure 11. This figure shows two pre-dip flares (panels a and b) and two QPP candidates (panels c and d).

Fig. D.2 Additional flare detections similar to Figure 11. Panel a) shows a flare with a ‘peak-bump’ profile, panels b) and c) show flares with a ‘bump-peak’ profile, and panel d) shows a flare with a ‘flat-top’ profile.

Appendix E Flare-star parameters correlations

Fig. E.1 Correlations between the denoised flare parameters and the ones of its star, for simple flares (blue) and individual components of complex flares (orange). Linear fits are shown for indication.

Fig. E.2 Correlations between the flare parameters and the ones of its star, for flares recovered in original (black) and denoised (red) light curves. Linear fits are shown for indication.

References

All Tables

All Figures

	Fig. 1 Parameters of our stellar sample. From top to bottom: Histograms for the Gaia magnitude, Teff, log(v sin i), radius, and distance.
In the text

	Fig. 2 Distribution of the visit duration within our sample.
In the text

	Fig. 3 Flow diagram of the procedure we followed.
In the text

Fig. 4 Averaged noise representation for each spectral type. In the top row, the dashed red curves correspond to a normal distribution, and the black histograms correspond to the noise distribution histograms in the CHEOPS light curves. In the bottom row, the dashed red lines correspond to the quantiles of a normal distribution, and the black curves correspond to the quantiles of the noise distribution in the CHEOPS light curves.

In the text

	Fig. 5 Injected flare energy matrix as a function of the injected flare amplitude, FWHM, and spectral subtype.
In the text

Fig. 6 Step-by-step representation of the flare injection, noise injection, wavelet denoising, and flare recovery process. Only 20 flares are shown to better visualise the flux variations in the light curve. In this example, the wavelet-based denoising was conducted with the Haar wavelet and a decomposition level of 2. Only 16 out of 19 flares are recovered in this portion of the light curve. The recovered and lost flare candidates are highlighted with a green and red background, respectively.

In the text

	Fig. 7 Flare recovery rate as a function of the injected flare amplitude, FWHM, and spectral subtype after denoising with the Haar mother wavelet and a decomposition level of 2. The top row corresponds to the recovery rates before denoising, and the bottom row shows the recovery rates after denoising.
In the text

	Fig. 8 Flare recovery rate, averaged across all spectral subtypes, as a function of the selected mother wavelet and decomposition level.
In the text

	Fig. 9 Evolution of confusion matrix metrics as a function of the injected flare energy. The TPR is shown in blue, the FPR in orange, the TNR in green, and the FNR in red. Metrics for the original light curves are represented by dashed lines, while those for the denoised light curves are shown as solid lines. The vertical black dotted lines indicate the four detection thresholds.
In the text

	Fig. 10 Confusion matrices obtained from the flare injection and recovery process for the original (left) and denoised (right) light curves. The matrices are shown at the four detection thresholds, corresponding to TPR rises above 10% (t1 and t2) and reaching 90% (t3 and t4), as defined in the text.
In the text

Fig. 11 Comparison of the flare recovery before (left) and after (right) denoising. In each panel, the original light curve (black data points) is compared to the denoised light curve (red data points). The vertical bars indicate the photometric uncertainty associated with each flux measurement, as estimated during the photometry extraction. The star on which each event was detected is mentioned in the corresponding panel title. The green shaded areas correspond to the data points flagged as belonging to a flare candidate by the flare detection algorithm, and the red shaded areas correspond to a non-detection in the same portion of the light curve. The solid lines represent the fitted flare components. Below each light curve, we show the residuals after subtracting the best-fit flare model. In the top right corner of each panel, the ${\tilde{\chi }}^2$ calculated on the portion of the light curve we show is given as indication. Each panel group corresponds to a specific case: (a) Simple flare detected in the original and denoised light curves. (b) Simple flare detected only after denoising. (c) Complex flare with all components recovered in the original and denoised light curves. (d) Complex flare with an additional component recovered after denoising.

In the text

	Fig. 12 Percentages of peaks observed per flaring event. The black and red bars correspond to flares recovered in original and denoised light curves, respectively. We compare our obtained percentages to the ones of Bruno et al. (2024), displayed as dashed blue bars.
In the text

	Fig. 13 Comparison of the number of peaks detected in each flare before and after denoising. Each point corresponds to a single flare. The red dashed line indicates the 1:1 relation, where the number of peaks remains unchanged. Points above the line represent flares that gained components after denoising, while those below lost components.
In the text

	Fig. 14 Amplitude vs duration of newly detected flares identified only after denoising. Each point represents a flare, with colours indicating the number of components retrieved from the flare breakdown.
In the text

Fig. 15 Corner plot of the different flare properties, namely relative amplitude, duration, measured energy, FWHM, impulse, and peak luminosity. The black data points represent metrics of individual flare components detected before denoising, while red data points correspond to those detected after denoising. A linear fit is provided for reference in each distribution. The discrete bins for flare duration and FWHM are due to the observational cadence equal to 3 seconds.

In the text

	Fig. 16 Comparison of Gmag, Tef f, radius, log vsini, distance and spectral type for flaring stars (red bars) against the entire stellar sample (white bars).
In the text

Fig. 17 Flare statistics by spectral subtype for the stars in our sample. The top panel shows the mean number of flares per star, the middle panel displays the mean flare energy, and the bottom panel presents the mean flare equivalent duration. Vertical bars represent the interquartile range (1st to 3rd quartiles), while horizontal bars indicate the median values. In the M5V bin of the bottom panel, only one star is included, resulting in equal values for the 1st quartile, median, and 3rd quartile.

In the text

	Fig. 18 Histograms of the number of flares observed per energy bin. The grey shaded histogram represents simple flares in denoised light curves, and the red shaded histogram represents individual components of complex flares. The dashed black and red histograms correspond to the same in the original light curves.
In the text

	Fig. 19 Correlations between flare parameters and stellar parameters. The blue and orange data points correspond to simple flares and individual components of complex flares, respectively. Each panel includes a linear fit for reference.
In the text

Fig. 20 Flare frequency distributions for the flares recovered in our original (black) and denoised (red) samples. We separate FFDs into simple flares (top left), individual components of complex flares (bottom left), flares recovered on partially convective stars (top middle), flares recovered on fully convective stars (bottom middle), and the entire flare samples (top right). We also provide a comparison between the simple flares (blue) and individual components of complex flares (orange) recovered in the denoised light curves (bottom right). Solid lines indicate linear fits to the double-logarithmic FFD, extrapolating into regimes not directly observed.

In the text

Fig. 21 Complementary cumulative distribution functions for recovered flare energies in denoised light curves. We separate flares into simple events (blue) and individual components of complex flares (orange), as well as flares from partially (green) and fully convective stars (purple). The full flare sample is displayed in black. Observed distributions, power-law fits, and log-normal fits are displayed with solid lines, dashed lines, and dotted lines, respectively. The grey dashed vertical lines indicate the energy of 90% recovery (t3) and the beginning of the energy range used for fitting. In each panel, the legend indicates the log-likelihood ratio R between the power-law and log-normal fits and its associated p-value.

In the text

	Fig. C.1 Flare recovery rate as a function of the selected mother wavelet and decomposition level, separated per spectral subtype.
In the text

	Fig. D.1 Additional flare detections similar to Figure 11. This figure shows two pre-dip flares (panels a and b) and two QPP candidates (panels c and d).
In the text

	Fig. D.2 Additional flare detections similar to Figure 11. Panel a) shows a flare with a ‘peak-bump’ profile, panels b) and c) show flares with a ‘bump-peak’ profile, and panel d) shows a flare with a ‘flat-top’ profile.
In the text

	Fig. E.1 Correlations between the denoised flare parameters and the ones of its star, for simple flares (blue) and individual components of complex flares (orange). Linear fits are shown for indication.
In the text

	Fig. E.2 Correlations between the flare parameters and the ones of its star, for flares recovered in original (black) and denoised (red) light curves. Linear fits are shown for indication.
In the text