Statistical Signatures of Nanoflare Activity. II. A Nanoflare Explanation for Periodic Brightenings in Flare Stars observed by NGTS

Several studies have documented periodic and quasi-periodic signals from the time series of dMe flare stars and other stellar sources. Such periodic signals, observed within quiescent phases (i.e., devoid of larger-scale microflare or flare activity), range in period from $1-1000$ seconds and hence have been tentatively linked to ubiquitous $p$-mode oscillations generated in the convective layers of the star. As such, most interpretations for the observed periodicities have been framed in terms of magneto-hydrodynamic wave behavior. However, we propose that a series of continuous nanoflares, based upon a power-law distribution, can provide a similar periodic signal in the associated time series. Adapting previous statistical analyses of solar nanoflare signals, we find the first statistical evidence for stellar nanoflare signals embedded within the noise envelope of M-type stellar lightcurves. Employing data collected by the Next Generation Transit Survey (NGTS), we find evidence for stellar nanoflare activity demonstrating a flaring power-law index of $3.25 \pm 0.20 $, alongside a decay timescale of $200 \pm 100$ s. We also find that synthetic time series, consistent with the observations of dMe flare star lightcurves, are capable of producing quasi-periodic signals in the same frequency range as $p$-mode signals, despite being purely comprised of impulsive signatures. Phenomena traditionally considered a consequence of wave behaviour may be described by a number of high frequency but discrete nanoflare energy events. This new physical interpretation presents a novel diagnostic capability, by linking observed periodic signals to given nanoflare model conditions.


INTRODUCTION
Magnetic reconnection is a process occurring throughout the outer solar atmosphere, often visible in the form of solar flares. Energies associated with flares express a wide range of magnitudes and frequencies; from very large, but infrequent X-class flares (with X-ray flux exceeding 10 −4 W/m 2 at the Earth, or ∼ 10 31 ergs per event; Maehara et al. 2015), down to micro-and nano-flares, each with energies on the order of 10 −6 and 10 −9 , respectively, of a typical X-class flare, but with occurrence rates that are orders of magnitude more frequent than the large-scale events. Stellar flares with energies similar to, and exceeding those of our own Sun have also been observed in many observations of stellar sources (e.g., Lacy et al. 1976;Maehara et al. 2012;Shibayama et al. 2013;Jackman et al. 2018), predominantly occurring in stars with convective atmospheres, which is required to generate the magnetic fields responsible for reconnection to take place (Pedersen et al. 2017).
The relationship between flare energy and the frequency of occurrence is commonly described by a power-law (Aschwanden et al. 2000), which applies at both low and high flare energies (Aschwanden 2019). The power-law exponent governs the frequency, dN/dE, of flaring events with an as-sociated energy, E, through the relationship, where α represents the power-law index. Low energy solar and stellar flares have long been a topic of wide interest. The power-law relation dictates that low-energy flares will be many many times more frequent than larger events. Parker (1988) proposed the power-law index is an indicator of the role of magnetic reconnection in maintaining the multimillion degree solar corona, with α > 2 allowing low energy (but highly frequent) nano-flares to supply sufficient thermal energy to the outer solar atmosphere to maintain its elevated temperatures.
Low energy stellar flares have been investigated by a number of authors (e.g., Hudson 1991;Robinson et al. 1995Robinson et al. , 1999Kashyap et al. 2002;Güdel et al. 2003;Güdel 2004;Welsh et al. 2006;Reale 2016, to name but a few). Much like their solar counterparts, there has been no clear consensus on the flaring rates of small-scale stellar flares, with the proposed power-law indices in the aforementioned studies spanning the range 1.5 ≤ α ≤ 2.7.
A review by Güdel (2004) suggested that power-law indices with α > 2 may be present in M-dwarfs. Butler et al. (1986) reported the presence of small-scale microflares in dMe flare stars observations that had previously been considered quiescent. Other authors (e.g., Brasseur et al. 2019) have investigated near-ultraviolet (NUV) flare events, with powerlaws of α = 1.72 ± 0.05 uncovered. The authors concluded that NUV flare mechanics are governed by the same physical processes as captured in solar events. Optical microflare signatures on M-dwarfs have also exhibited short time-scale variability as discussed by Schmitt et al. (2016), finding flare rise timescales on the order of seconds, with flare signatures of comparable brightness to its quiescent B-band luminosity. These studies highlight the growing interest in small-scale flare events, and demonstrate the synergy between stellar and solar observational and modeling efforts.
However, there is a gap in the current literature, with few studies investigating the role of nanoflares on other stellar sources. Falla & Potter (1999) examined the production of nanoflare energies in the X-ray emission of RS CVn systems. The authors concluded that while nanoflares may be produced in these stars, current observational limits would prohibit the direct detection of nanoflare events in the Xray band. True to the predictions of Falla & Potter (1999), currently the lowest energy stellar flares that have been directly observed are on the order of 10 28 ergs Benz & Güdel 2010), which are orders of magnitude above the traditional range of individual nanoflare energies. It is generally predicted that the flare occurrence rate will be higher on magnetically active stars, such as dMe flare stars (Walkowicz et al. 2011). As such, nanoflares may be even more frequent on these stellar sources when compared to the Sun, thus producing power-law indices substantially larger than estimates for the solar case.
Direct observation of solar nanoflares has also remained a challenging endeavor, with their signals lying below the noise floor of current generation instrumentation. As a result, researchers have had to turn their attention to other approaches, such as spectroscopic techniques to compare the scaling between kinetic temperatures and emission measures of coronal plasma (e.g., Klimchuk & Cargill 2001;Sarkar & Walsh 2008;Sarkar & Walsh 2009;Bradshaw et al. 2012), comparisons drawn between EUV and X-ray emission (e.g., Sakamoto et al. 2008;Vekstein 2009), or the examination of the time delays between different temperature-sensitive EUV imaging channels (e.g., Viall & Klimchuk 2011, 2016. In addition, Terzo et al. (2011) and Jess et al. (2014) employed statistical techniques to provide evidence of solar nanoflares. These statistical approaches are further developed in the recent work by Jess et al. (2019), who infer the presence of nanoflares in a seemingly quiescent solar dataset by comparing intensity fluctuations extracted from high time resolution imaging with those from Monte Carlo synthetic lightcurves designed to replicate the presence of small-scale nanoflare events. Jess et al. (2019) further suggested that similar nanoflare statistical techniques could also be directly applied to high time resolution observations of stellar sources, i.e., to modernize the work of Audard et al. (1999) and Kashyap et al. (2002) through the comparisons of intensity fluctuations with nanoflare-specific simulations.
On the contrary to the flare frequencies predicted by the dN/dE power-law relationship, several studies have documented evidence for 'periodic' brightness variability through the examination of stellar intensity fluctuations, with periods ranging between 1 − 1000 s , (Andrews 1989;Rodríguez et al. 2016;McLaughlin et al. 2018). These periodic brightenings are of uncertain origin, but are believed to be linked to ubiquitous p-mode oscillations or other magnetohydrodynamic (MHD) wave behavior (e.g., Aschwanden et al. 1999;Nakariakov & Verwichte 2005;Nakariakov et al. 2010;McLaughlin et al. 2018) generated in the convective layers of stars. The link to p-mode oscillations is due to a comparable period range (1 − 1000 s), in addition to them being observed during periods of quiescence (i.e., no associated macroscopic flaring signatures).
In a number of publications, Andrews (1989Andrews ( , 1990a examined dMe flare stars across a range of conditions, from immediately after large-scale flare events, to during relatively long periods of quiescence, and found that the dMe flare stars exhibited small periodic brightenings, on a scale of seconds to minutes. The author interpreted these periodic signals as a likely consequence of MHD wave behavior, as the periodic signals were observed during times of quiescence, with no impulsive activity witnessed in the time series. A follow up study by Andrews & Doyle (1993) investigated whether flaring events can re-produce signals with 1 − 1000 s periodicities, and suggested that while individual small-scale flares may contribute to such signatures, they were unable to provide sufficient evidence to directly link flaring events to the periodic signals.
However, flare-related variability giving rise to periodic phenomena has been documented across a range of solar observing sequences. McLaughlin et al. (2018) discuss selfoscillatory flaring (perhaps due to magnetic dripping, as discussed by Nakariakov et al. 2010), which can produce a periodic signal, despite non-periodic driving. Additionally, Arzner & Güdel (2004) discussed flare clustering, and the relationship between the mean flaring interval and expected count rates. This led Jess et al. (2019) to speculate that smallscale flaring may have a quasi-periodic nature, due in part to the power-law governing its occurrence rates. With this in mind, the superposition of hundreds or thousands of (quasi-) periodic nanoflare signatures each second may give rise to a periodic brightness signal, without any 'flare-like' impulsive signatures seen in the corresponding stellar lightcurve. By combining the statistical parameterization techniques developed for solar nanoflare detection with a novel Fourier spectral analysis, here we investigate stellar nanoflare signals and their potential role in periodic brightenings found in stellar lightcurves.

OBSERVATIONS WITH NGTS
The impulsive rise and subsequent decay phase for solar nanoflares are on the order of tens to hundreds of seconds (Jess et al. 2019). Stellar flare decay rates on UV Ceti-type stars are around one order of magnitude shorter than for the Sun, leading to even faster signal evolution on the order of tens of seconds (Gershberg 1975). As a result, high-frequency resolution and a short temporal cadence are required to fully capture these dynamic signals. Jackman et al. (2018Jackman et al. ( , 2019aJackman et al. ( ,b,c, 2020 employed the high cadence of the Next Generation Transit Survey (NGTS; Wheatley et al. 2018) to apply techniques developed for solar flare analysis to stellar flare oscillations, inspiring our use of the NGTS to extend statistical solar nanoflare techniques to stellar lightcurves. The NGTS is a ground-based array of 12 telescopes that scan the sky in the optical domain searching for transiting exoplanet signals, but has also become a platform for stellar flare analyses. The NGTS has cadence of ≈12 s, providing a Nyquist frequency of ≈41.6 mHz, with the observations spanning up to hundreds of thousands of frames for a single star.
When searching for signatures of nanoflare activity we extracted lightcurves for M-type stars. M-dwarf flares have a higher contrast due to their lower quiescent background flux than is typically seen on G and K stars (Günther et al. 2020). This increased contrast is essential to capture nanoflare signals below the noise floor. Additionally, M-star flares have a strong contribution in white-light (Walkowicz et al. 2011), ideal when utilizing datasets from optical surveys, i.e., the NGTS. These benefits outweigh an increased photometric noise level, which is itself minimized by leveraging the large number statistics of statistical nanoflare analysis. Finally, as these are flare-active stars , flare occurrence rates will be higher than in 'solar-like' stars. This means M-dwarf stars are likely to provide the best conditions for the manifestation of detectable nanoflare signals. Specifically, the stars NGTS J030047.1-113651, NGTS J030415.6-103712, and NGTS J031800.1-212036 were chosen as each of these had more than 10 5 datapoints available for study, hence maximizing the available number statistics for our analyses.
As our scientific analyses revolves around flare-active Mtype stars, it was deemed important to also examine nonflare active stars, which can act as a control test to ensure our data analysis techniques are not incorrectly mistaking residual systematic signals as evidence for stellar nanoflares. Since A-type stars are absent of a convective zone, their resulting lack of flare-like behavior provides an ideal set of complementary data products. Some recent studies (Balona 2012;Fossati et al. 2018;Balona 2020) do suggest A-stars are capable of flaring, but this has also been disputed (Pedersen et al. 2017). If the observed signals are indeed A-star flares, then only extremely energetic flares have been observed; Balona (2020) discuss A-star flares with energies in the range 10 35 − 10 36 ergs, 10 orders of magnitude above traditional nanoflare activity. If only highly energetic events can rise above the high background luminosity on A-stars, this would explain the rarity of A-star flare observations. As such, low-energy nanoflaring would be entirely lost within the lightcurves of these stars due to the minimal contrast invoked, meaning A-type stars would appear quiescent at small-scale flare energies in the NGTS datasets, regardless of their true flaring behavior. This meant A-type stars cannot exhibit a signal consistent with nanoflares. As such, we examined the A-type stars NGTS J025840.5-120246, NGTS J030958.4-103419, and NGTS J030129.4-110318. It is important to note that A-type stars have a very different spectral energy distribution to M-dwarf stars, so are not a conventional choice for relative photometric comparison. To ensure robust null testing, we also examined low-activity K stars which have a more comparable SED to M-type stars (i.e., choosing similar spectral types as is standard for photometric comparison, e.g., Amado et al. 2000). The lowactivity K-type stars were chosen over low-activity M-dwarf stars due to their higher luminosity, leading to decreased low-energy flare contrast when compared to the M-types. While the low-activity K-type stars could theoretically have some weak nanoflaring signature present, it would be minimized compared to the M-types, so this still serves as a valid null test. We used the K2V type stars NGTS J030000.7-105633 , NGTS J030848.9-112217 , and NGTS J030538.9-114145. These K-stars were low-activity and had no macroscopic flare events in their observed timeseries.
All of the A-type stars , K-type stars , and two of the three Mtype stars were obtained from the same observational field (NG0304-1115) and camera (809), hence ensuring consistency across the processed A-, K-, and M-type data sequences. NGTS J031800.1-212036 was from a different field (NG0313-2230), but had noise statistics, magnitude, and stellar parameters consistent with the other M stars used in the present study.
The magnitudes of the stars employed were comparable (see Table 1). This was important to ensure the noise statistics were consistent across the stars. The majority of the stars were around mag 13. At this magnitude, the dominant noise source is photon noise (see Figures 3 & 14 of Wheatley et al. 2018), with scintillation noise only becoming dominant at the highest frequencies in the data, which are beyond the typical p-mode periodicities we are investigating (Osborn et al. 2015). The A star NGTS J030129.4-110318 was the brightest, with an NGTS magnitude of 11.69. At this magnitude, scintillation became a dominant source of noise. This allowed us to investigate the effect of increased scintillation noise on our analysis techniques. We utilized the stellar parameters from the TESS Input Catalog Version 8 (TIC V8) (Stassun et al. 2018), along with the initial spectral classification provided via Spectral Energy Distribution (SED) fitting performed by the NGTS pipeline (see section 5.1.1 in Wheatley et al. 2018) to assign the spectral types. See Table A4 in the Appendix, for this and other observational parameters (i.e., GAIA Source ID, RA, Dec, mass, radius, luminosity, distance, approximate macroscopic flare rates, and the log Lx L Bol ratio, where L x and L Bol are the x-ray and bolometric luminosities, respectively).
The only M star with an x-ray luminosity measurement was NGTS J030047.1-113651 which had a x-ray flux measurement available from the 4XMM XMM-Newton Serendipitous Source Catalog (Webb et al. 2020). This corresponded to an x-ray luminosity of 6.47 × 10 28 ergs s −1 . The ratio of x-ray luminosity to the bolometric luminosity is an indication of the activity rate of the star. We find log Lx L Bol = −3.09 ± 0.21, which compares to the literature values for a young and active M-type star, with saturated x-ray emission of log Lx L Bol ∼ −3 (Kastner et al. 2003;López-Santiago et al. 2010).
The lightcurves were corrected for background and flatfielded according to the NGTS data reduction pipeline described in Wheatley et al. (2018). This pipeline provides a relative error in the flux at each point in the time series. These error bars are affected by cloudy weather and high airmass. Any fluctuations in this error exceeding 1σ above the mean value were removed, resulting in ∼10% of each time series being omitted. This removed any data that had statistically significant increases in its associated flux uncertainties, therefore preventing any large flux errors (largely due to poor seeing conditions) from contaminating the final time series. Next, the lightcurves extracted for each observing sequence were examined for the presence of macroscopic flare signatures, something which occurred in ∼0.2% of the remaining M-type time series (i.e., following the removal of datapoints exceeding 1σ in their relative flux errors). To isolate the macroscopic brightenings, each lightcurve was searched for emission signatures exceeding 3σ above the mean value, which lasted continually for a minimum of 1 minute (5 datapoints). Based on a normal distribution, the probability of this occurring by chance is 2 × 10 −13 , and hence allowed for the robust detection of intensity fluctuations resulting from macroscopic flaring activity. Once the larger scale flare signatures had been identified, they were cut from the time series using an interval of ±5 minutes (25 datapoints) from the first and last detection above the 3σ threshold. For consistency, the same processing steps were applied to the A-type and K-type stellar lightcurves, but no macroscopic brightenings were found for these sources. The number of macroscopic flares removed were used to calculate approximate flare rates for the M stars. These are listed in Table A4 in the Appendix. The flare rates were of a comparable magnitude for the three M stars, with rates of 0.012, 0.027 and 0.003 flares removed per hour for NGTS J030047.1-113651, NGTS J030415.6-103712, and NGTS J031800.1-212036 respectively. Combining this with the x-ray luminosity of NGTS J030047.1-113651 being that expected for a young and active M-star, we extrapolate that all three M stars are macro-flare active, with roughly comparable activity levels. Upon completion of the lightcurve filtering, the lowest number of datapoints remaining was 97 060. To ensure consistency across all subsequent analyses, each of the other eight M-K-and A-type time series were cropped to the same 97 060 datapoints.
Once the macroscopic flare signatures had been extracted from the time series, each of the remaining lightcurves were normalized (night by night) by subtracting a linear line of best fit that was derived from the corresponding time series. Next, the lightcurves were divided by their respective standard deviations, σ N , providing time series of fluctuations around a common mean that can be readily cross-compared with other star types and data products. This statistical treatment resembles common Z-score testing, which is a statistical technique regularly employed in physical and social sciences (Sprinthall 2012). To ensure that the output data products did not contain any long-term and/or instrumental trends that are not accounted for using the initial preparatory routines, we subsequently detrended these data products using low-order polynomial fits.

ANALYSIS AND DISCUSSION
As documented by Terzo et al. (2011) and Jess et al. (2014, time series commonly referred to as 'quiescent' may in fact contain a wealth of small-scale nanoflare signatures that are embedded within the inherent noise of the photometric signals. It is possible to uncover these signatures through statistical analyses of the intensity fluctuations. We employ the same techniques described by Jess et al. (2019) to attempt to recover nanoflare signatures in our M-dwarf lightcurves.
The dominant source of noise in seemingly 'quiescent' NGTS lightcurves will be shot noise, which follows a Poisson distribution . The fluctuations will be random, and in the limit of large number statistics, will demonstrate equal numbers of positive and negative fluctuations about the time series mean (Frank 2009). Therefore, plotting a histogram of the inherent shot noise fluctuations for a truly quiescent time series would produce a symmetric distribution, with the mean and median centered at zero. Any subtle offsets and/or asymmetries to this idealized case may be interpreted as signatures of impulsive events, with subse- Figure 2. Histograms of intensity fluctuations, each normalized by their respective standard deviations, σN , for the NGTS J025840.5-120246 (A-type; above) , NGTS J030000.7-105633 (K-type; middle) , and NGTS J030047.1-113651 (M-type; below) lightcurves. A standardized Gaussian profile is overplotted in each panel using a red dashed line for reference. The M-type distribution has a negative median offset with respect to the Gaussian, in addition to elevated occurrences at ∼ 2 σN , which is consistent with the statistical signatures of nanoflare activity. On the other hand, the A-type and K-type intensity fluctuations provide no signatures of flare activity, with the resulting distribution remaining consistent with the presence of photon-based shot noise. Zoomed insets highlight the ranges spanning −0.4 ≤ σN ≤ 0.0 and 1.7 ≤ σN ≤ 2.2, where M-type negative median offsets and occurrence excesses, respectively, are found. The blue and gold lines display the derived distributions, The M-type exhibited a small dip below the idealized Gaussian at around −0.90 σN , which are not seen in the Aand K-star. We believe this is connected to the negative median offset signal, which is causing a consequential dip elsewhere in the statistical distribution, but the exact nature of the signal is unknown. quent exponential decays, embedded within the noise floor of the lightcurve (Terzo et al. 2011).
As discussed by Jess et al. (2019), nanoflares give rise to two distinct signals in the resulting intensity fluctuation histograms. The first is a negative median offset, whereby the median value of the histogram is < 0 σ N . This is a characteristic signal associated with an exponentially decaying signature, i.e., the decay phase following an impulsive deposition of energy occurs over a longer timescale, hence providing more fluctuations that are beneath the elevated signal mean caused by the impulsive event. The second signature is an excess of fluctuations at ∼ 2 σ N , which is caused by the impulsive nature of the nanoflare intensity rises, and gives rise to an asymmetric distribution that can be benchmarked using Fisher skewness coefficients. As the evolution of a nanoflare produces an almost discontinuous increase in the lightcurve intensity, a distinct positive peak manifests in the resulting histogram of intensity fluctuations. Therefore, a seemingly quiescent lightcurve exhibiting both of these signals is a strong candidate to contain embedded nanoflare signatures. Additionally, we benchmark the shape and widths of the distributions through calculation of the histogram kurtosis values, in addition to the ratio of its full-width at eighthmaximum to that of its full-width at half-maximum (i.e., Table 2. Characteristics of the intensity fluctuation histograms associated with the A-, K-, and M-type NGTS sources. Note that a standard Gaussian distribution will demonstrate ζ = 1.73, hence deviations from this provide an indication of the intensity fluctuation occurrences taking place close to, and far away from the time series mean.

NGTS
, which is defined as 'ζ' for simplicity (see Jess et al. 2019, for a more thorough overview of this key statistical parameter). Note that a standard Gaussian distribution will have ζ = 1.73, hence deviations from this provide an indication of the intensity fluctuation occurrences taking place close to, and far away from the time series mean. .5-120246, NGTS J030000.7-105633, and NGTS J030047.1-113651, respectively. As expected, the non-flare active A-type star and low-activity K-star show little variation from the standardized Gaussian distribution (dashed red lines in Fig. 2), with median offsets of 0.000 ± 0.004 σ N and no visible excess at ∼ +2 σ N . This suggests that the A-type and K-type stars have no embedded nanoflare characteristics, and therefore reiterates their importance as a control test for subsequent M-type star analysis. On the other hand, the M-type star displays both of the characteristic nanoflare signatures, with a negative median offset equal to −0.050 ± 0.004 σ N , and a visible occurrence excess at ∼ +2 σ N , culminating in an associated positive Fisher skewness value of 0.031 ± 0.008 that is above the expectations of a pure Gaussian distribution.

NGTS Datasets
The other candidate stars exhibited consistent signals, with the M-type stars showing histogram signatures consistent with nanoflare activity, while the A-and K-type stars showed no indication of impulsive behavior beneath the noise floor. A-type star NGTS J030129.4-110318 exhibited a small positive skew of 0.003 ± 0.008, but the associated uncertainty makes this less definitive when compared to the positive skewness values exceeding 0.040 for some M-type sources.
Furthermore, NGTS J030129.4-110318 also demonstrated zero median offset, remaining inconsistent with a distribution comprised of impulsive events followed by gradually decaying tails. This star had a much larger deviation from Gaussian statistics, evidenced by a kurtosis value of 0.688±0.016, and ζ ratio of 1.814 ± 0.015. This deviation from Gaussian statistics is due to the increased brightness of this star (see Table 1) compared to the other candidates, resulting in scintillation becoming a more significant noise source (Osborn et al. 2015;Wheatley et al. 2018). It is important to note that while the scintillation noise produces statistics offset from a Gaussian, it is still distinct from the characteristic signatures of nanoflaring. This highlights the robustness of the statistical nanoflare analysis. The characteristics derived for all 9 stellar sources are documented in Table 2. The M-type stars exhibited a small dip below the Gaussian around −0.90 σ N , which was not seen in the A-or K-stars. We believe this is connected to the negative median offset signal, which is causing a dip elsewhere in the statistical distribution, but the exact nature of the signal is unknown. Future investigation could uncover the source of this dip, and potentially use it as a further diagnostic.
Employing the high time resolution and long duration imaging sequences of the NGTS data products has enabled us to provide the first tentative evidence of nanoflares occurring on stellar sources (see, e.g., Fig. 2 and Table 2). However, while the statistical signatures derived for the NGTS M-type lightcurves resemble those expected for nanoflare activity, they do not provide any indication of the specific underlying plasma conditions at work.
As previously demonstrated by Andrews (1989), Rodríguez et al. (2016), andMcLaughlin et al. (2018), smallscale brightenings -here hypothesized to be the result of nanoflare activity -often give rise to periodic signatures in the corresponding lightcurves. This has also been observed . It can be seen that theA-and K-type PSD are relatively flat, with no clear power enhancements,apart from slight enhancement in the K-type star, in the range of 1 − 10 mHz, indicative of the expected p-mode oscillations seen in Solar-like stars. Contrarily, the M-type PSD has a primary power peak at ≈0.8 mHz, followed by decreasing spectral power exhibiting a spectral slope of β = −0.30 ± 0.05, followed by numerous power peaks in the range of 3 − 10 mHz, consistent with previous links to stellar p-mode spectra.
in the case of small-scale solar activity (Terzo et al. 2011). To investigate the manifestation of periodicities in the stellar lightcurves, time series were extracted for each star that contained the maximal number of successive frames, where no breaks resulting from problematic flux calibrations, macroscopic flare events, or day/night cycles were present, i.e., the longest consecutive number of frames consistent across the 9 stars. The lowest number of viable consecutive frames was 2316 from M-type star NGTS J031800.1-212036. As such, each of the remaining five lightcurves were cropped to an identical 2316 datapoints (≈ 27 800 s duration) so that the final A-,K-, and M-type time series had identical lengths, helping to ensure consistency between both the Nyquist frequency and frequency resolution in the subsequent analyses.
Each of the six extracted NGTS lightcurves were passed through a Fast Fourier Transform (FFT) to determine whether power exists at frequencies synonymous with a typical p-mode spectrum (often in the range of 1 − 1000 s; Kjeldsen et al. 1995;Guenther et al. 2008;Handler 2013;Di Mauro 2016). The input data resulted in a Nyquist frequency of ≈41.6 mHz being complemented by a frequency resolution, ∆f = 0.0356 mHz, in the corresponding FFTs. However, it must be pointed out that a strictly periodic wave signal would not manifest as median offsets and/or asymmetries in the fluctuation histograms documented in Figure 2, since the evolution of a purely sinusoidal wave signal is symmetric about its given mean. The resulting Fourier power spectra were transformed into power spectral densities (PSDs) following the methods defined by Welch (1961) and Vaughan (2013).
Following the generation of PSDs from the nine NGTS lightcurves, we find that the A-, K-, and M-type sources exhibit consistent and distinct features in their corresponding PSDs, with examples depicted in Figure 3. The upper , middle, and lower panels of Figure 3 display the PSDs for the A-, K-, and M-type stars NGTS J025840.5-120246,  Figure 3 that the A-type and K-type spectra are relatively flat across all frequencies with no evidence of distinct peak frequencies. The K-type does show some slight power enhancement between ≈ 1 − 10 mHz, consistent with stellar p-mode oscillations, as have been previously observed in K-type Solar-like stars (e.g. Chaplin et al. 2009) . The Mtype PSD exhibits more pronounced fluctuations across the frequency domain. In the lower panel of Figure 3, the solid red line highlights the presence of a primary power peak at ≈0.8 mHz, followed by a gradual decline in power as the frequency increases. This reduction in power, as a function of frequency, can be represented by a spectral slope, β, following the form f β . In the lower panel of Figure 3, the spectral slope is calculated to be β = −0.30 ± 0.05. For each M-type star the position of the primary peak, and its associated spectral slope, were calculated. The primary peaks (or 'Turning Point') were found in the range of 0.6 − 0.9 mHz, with the corresponding spectral slopes calculated to span −0.30 ≤ β ≤ −0.26. Once the spectral slopes had been calculated, they were subsequently subtracted from each PSD to better highlight power fluctuations above the background level (similar to the processing undertaken by Krishna Prasad et al. 2017). Following the detrending of the PSDs, the frequency demonstrating maximal power above the background was subsequently extracted, and found to reside in the range of 2.58 − 3.09 mHz for the M-type stellar sources, which is consistent with previous interpretations related to the presence of p-mode oscillations (Andrews 1989(Andrews , 1990a. The specific characteristics derived from the M-type PSDs are displayed in Table 3.The Figure B8 in Appendix B plots all the PSD trendlines on one plot for clarity. With the lightcurve intensity fluctuations statistically benchmarked, and the corresponding power spectra uncovered, we now generate Monte Carlo nanoflare simulations that have been tailored for stellar sources. This will enable direct comparisons to be made between the observed and simulated time series (for both the statistical fluctuations and the power spectra features), which will help quantify the specific plasma parameters at work in each of the stellar sources. The modeled time series will be cropped to the same length as the NGTS time series, i.e., 97 060 data points (at a cadence of ∼ 12s) for the statistical analysis, and 2316 data points for the PSDs, to ensure consistency in their number statistics.

Stellar simulations
We adapt the Monte Carlo simulations described by Jess et al. (2019) to synthesize the intensity time series expected for a broad range of initial plasma conditions. The adaption process first necessitated altering the 'area' over which the simulations took place. In the work of Jess et al. (2019), a two-dimensional image was generated to simulate data acquired by the Atmospheric Imaging Assembly (Lemen et al. 2012) onboard the Solar Dynamics Observatory (Pesnell et al. 2012), where the pixels had an area of approximately 10 15 cm 2 . However, considering our one-dimensional lightcurves contain no resolvable spatial information of the stellar sources, we need to increase the modeled area to represent the entire Earth-facing surface area of the star. This meant setting the pixel area to around 10 21 cm 2 , or approximately 10% of the surface area of the Sun, which corresponds to the surface area of a typical M-dwarf stellar source (Reid & Hawley 2005). Note that we use the entire Earth-facing surface area. While larger flares require specific high energy magnetic conditions (e.g., large-scale spots that may only cover a small proportion of the stellar surface and be more aligned with the stellar equator), it is expected that nanoflares can effectively occur anywhere across the stellar atmosphere, requiring only small-scale magnetic activity to trigger them. This area, along with an exposure time (10 s) and final cadence (12 s) matched to the NGTS observations, was used to re-compute the number of flaring events expected (following Equation 1 and the work by Aschwanden et al. 2000;Parnell & Jupp 2000), for a given power-law index, α, and across a specific time interval. The quiescent flux of the M-stars was used to generate the underlying Poisson noise in the flare models and the nanoflare energies were then calibrated to this noise level, following the steps taken in Jess et al. (2019) .
The flare energies included in our model spanned 10 22 − 10 25 ergs, placing them within the energy regime synonymous with solar nanoflares. dMe flare stars are not 'solarlike'; arguably, the energy span of M-type stellar nanoflares may be orders of magnitude larger than for the solar case, due to the increased flare energies associated with M-type stars. However, various authors (Falla & Potter 1999;Robinson et al. 1999;Güdel et al. 2002) have applied the solar energy span derived by Aschwanden et al. (2000) directly to stellar investigations, hence we follow the same convention for consistency. The conversion of flare energies to peak detector counts, DN , is performed via a direct one-to-one scaling relationship. According to Yang et al. (2017), the flare energy is linear with area, which is linear with flux, assuming a constant black-body emission temperature. As flares emit primarily in optical and UV wavelengths (Neidig 1989;Woods et al. 2006;Schmitt et al. 2016), whitelight observations are likely to capture the resulting nanoflare emission (Kretzschmar 2011), particularly for M-dwarf flares which emit strongly in white light (Walkowicz et al. 2011), resulting in a DN ∝ E relationship. This is similar to the pulse-heating model proposed by Jess et al. (2019), whereby DN ∝ E 4/3 . For the energy range relevant to nanoflares (i.e., spanning only 3 orders-of-magnitude; 10 22 − 10 25 ergs), the differences between the linear scaling and pulse-heated models is relatively small. However, if accurate modeling and replication of full-scale flaring events (i.e., 10 22 − 10 31 ergs) is required, then more precise whitelight emission models would need to be developed. (Procházka et al. 2018).
Due to the large spatial integration (≈ 10 21 cm 2 ), the simulations are more computationally intensive than described by Jess et al. (2019). As a result of integrating over the entire stellar disk, the generation and superposition of hundreds of thousands of independent nanoflare events becomes a more time consuming endeavor, requiring approximately 300 s on a 2.90 GHz Intel Xeon processor to generate a synthetic NGTS time series incorporating 97 060 individual frames (∼13.5 continuous days of data at a cadence of 12 s). An example depicting the generation of a synthetic NGTS lightcurve is shown in Figure 4. Here, the lightcurve is cropped to a 36 000 s interval to more clearly reveal its constituent components. The upper panel of Figure 4 displays (black line) modeled flaring events using a power-law index α = 3.25 and a decay timescale (i.e., reflecting the e-folding time of the flare decays) of τ = 245 ± 24.5 s. Note that the decay timescale varies by ±10% (i.e., τ = 245 ± 24.5 s) to allow for subtle variations in the mechanisms responsible for cooling in the immediate aftermath of the flaring events (Antiochos & Sturrock 1978). This nanoflare time series is the superposition of individually generated flare events. The red dots represent the background shot noise, which follows a Poisson distribution. According to the limits of large number statistics, this Poisson profile will transform into a Gaussian distribution, with ≈68.3%, ≈95.5%, and ≈99.7% of the noise fluctuations contained within the invervals of ±1σ N , ±2σ N , and ±3σ N , respectively. It is visible from the upper panel of Figure 4 that even larger flaring events, e.g. occurring at ∼200 s and ∼500 s, are contained within the noise envelope. Once the shot noise contributions have been added to the synthetic flaring signals, the resulting time series (lower panel of Fig. 4) mimics very closely typical stellar lightcurves (i.e. the NGTS lightcurves in Fig. 1), with the original nanoflare signal now indiscernible from the embedded noise. Figure 4 documents the steps taken to generate a synthetic lightcurve for a specific power-law index (α = 3.25) and e-folding timescale (τ = 245 ± 24.5 s). However, in order to more accurately constrain our observational findings using our synthesized models required us to repeat the processing steps documented in Figure 4 using a dense grid of nanoflare input parameters. Specifically, power-law indices spanning 1 ≤ α ≤ 4 (in intervals of 0.05) and e-folding timescales ranging across 5 ≤ τ ≤ 500 s (in steps of 5 s, consistent with previous estimations for solar nanoflares; Terzo et al. 2011;Jess et al. 2014) were employed. This produced 6100 final synthetic NGTS lightcurves, each with 97 060 datapoints to remain consistent with the observational NGTS time series, ensuring identical number statistics and allowing direct comparisons to be made between the observations and simulations.

Comparing Simulation to Observation
Each synthetically generated lightcurve was treated in an identical manner to that of the NGTS observations, whereby each of the 6100 simulated time series were detrended and normalized by their respective standard deviations, before generating their intensity fluctuation distributions and subjecting them to FFT analyses. It must be noted that there were no instances in any of the 6100 simulated time series where a sequence of 5 successive time steps exceeded +3σ N above the mean, hence highlighting the consistency between the simulated lightcurves and the final time series extracted from the NGTS observations. First, to compare the observational intensity fluctuation distributions depicted in Figure 2 to those extracted from the dense grid of simulation input parameters, we generated a number of statistical maps (Fig. 5) where the parameter values extracted from the intensity fluctuation histograms are displayed as a function of the power-law index, α, and the corresponding decay timescale, τ . These statistical benchmarks are the same as those calculated for the NGTS stars in Table 2, only now graphically displayed in a two-dimensional format to aid visual clarity.
The measured output parameters depicted in Figure 5 allows us to cross-correlate the observational signatures to those synthetically generated via the Monte Carlo modeling work, hence allowing us to estimate the specific plasma conditions (i.e. the α and τ values) responsible for the observational signatures. Importantly, the synthetic stellar lightcurves are consistent with those expected from solar modeling efforts (Terzo et al. 2011;Jess et al. 2014Jess et al. , 2019, whereby a negative median offset is coupled with an increase in the Fisher skewness value. From Figure 5 it can be seen that the majority of nanoflare conditions produce a negative median offset and positive Fisher skewness in the resulting statistical intensity fluctuation distribution, despite the presence of seemingly quiescent lightcurves (see, e.g., the lower panel of Fig. 4). Similar dips below the idealized Gaussian at approximately −0.90σ N (as were seen in the M-type stars) were exhibited in the simulations, suggesting these are linked to the embedded nanoflare signals.
When comparing the intensity fluctuation statistical outputs for the M-type stars to those derived from the Monte Carlo simulations, we found overlap in the median offset, Fisher skewness, kurtosis, and ζ ratio corresponding to two distinct plasma conditions governed by the flare power-law index, α, and the associated decay timescale, τ . The first set of self-similar parameters corresponded to α = 3.25 ± 0.15 and τ = 200 ± 100 s, while the second set of parameters consisted of α = 2.00 ± 0.15 and τ = 200 ± 100 s. These values highlight the fact that the observational M-type NGTS lightcurves show remarkable agreement with the statistical signals derived from Monte Carlo synthetic lightcurves consisting of nothing but nanoflare signals embedded in characteristic shot noise. Contrarily, the A-type and K-type stellar parameters do not map consistently onto the statistical parameters depicted in Figure 5, reiterating our interpretation that the A-type and K-type sources do not exhibit nanoflare signatures.
In order to further examine the link between nanoflare activity and periodic variability in the synthetic lightcurves, we generated PSDs for each of the 6100 simulated time series, which could then be compared directly with the PSD features found in the NGTS observations. To remain consistent with the observational PSDs depicted in Figure 3, we cropped the synthetic time series to 2316 datapoints to ensure the frequency resolution was maintained at ∆f = 0.0356 mHz. As the comparison between the observed and modeled intensity fluctuation distributions revealed a self-similar set of  Table 2 and Figure 2) compare to the modeled statistical distributions with overlapping parameters corresponding to α = 3.25 ± 0.15 and τ = 200 ± 100 s, in addition to α = 2.00 ± 0.15 and τ = 200 ± 100 s. statistical parameters corresponding to a power-law index α = 3.25 ± 0.15 and a decay timescale τ = 200 ± 100 s, we provide example PSDs for α = 3.25 and τ = 245 ± 24.5 s in Figure 6. Such Fourier analysis offers an additional paramaterization of the nanoflare signal, allowing us to resolve any ambiguities arising through examination of the statistical signatures alone.
The upper panel of Figure 6 shows the corresponding PSDs for both the raw nanoflare (red crosses) and Poisson-based shot noise (blue squares) signals. The solid black and gold lines in the upper panel of Figure 6 depict the trendlines for the nanoflare and shot noise signals, respectively, established over ±6 frequency elements (±0.427 mHz). It can be seen that at lower frequencies ( 5 mHz) the nanoflare signal dominates over the corresponding noise profile, while at higher frequencies the noise becomes dominant and begins to mask the frequency-dependent signals of nanoflare activity. The lower panel of Figure 6 displays the PSDs extracted from the final simulated lightcurve, where the nanoflare signal has been embedded within the synthetic noise profile.
To remain consistent with the lower panel of Figure 3, the black crosses represent the individual frequency-dependent power measurements, while the solid red line depicts a trendline established over ±6 frequency elements (±0.427 mHz). The similarities between the lower panels of Figures 3 & 6 are remarkable, exhibiting similar primary power peaks at ∼1 mHz, followed by a decrease in spectral power with increasing frequency, before finally demonstrating a number of power peaks within the range commonly associated with p-modes. It must be remembered that the A-type and K-type stellar sources provided flat and relatively featureless spectra, with no spectral slopes visible in their corresponding PSDs. Hence, the A-and K-type PSDs (see, e.g., the upper and middle panel of Figure 3) show no agreement with the synthetic PSD depicted in Figure 6, and serves as a further indicator that there is no nanoflare activity present on our A-and Ktype stellar samples.
In a consistent manner with how the M-type stellar PSDs were processed, each of the 6100 synthetic lightcurves were examined and their corresponding primary frequencies, spec- tral slopes, and dominant frequencies (following detrending by the computed spectral gradients) were calculated. In order to more readily display these sets of measured parameters, we display them in Figure 7 in a two-dimensional format as a function of the power-law index, α, and the corresponding decay timescale, τ . This is similar to the intensity fluctuation statistical measurements depicted in Figure 5, only Figure 7 now displays the corresponding parameters extracted from the analysis of the synthetic PSDs. Figure 7 documents interesting behavior of the key Fourier-based parameters as a function of the power-law index, α, and the corresponding decay timescale, τ . As the power-law index increases, the spectral slopes (upper-left panel of Fig. 7) associated with the PSDs begin to flatten. This is likely a consequence of increased energy being spread across the entire frequency spectrum as a result of the larger power-law indices (Jess et al. 2020). Previous work on turbulent cascades have revealed spectral slopes within the range of −2 ≤ β ≤ −1 in both solar and stellar plasma (Podesta 2011;Huang et al. 2017), believed to be a feature of wave behavior. The spectral slopes found in our simulation outputs varies largely within this range (−1.85 ≤ β ≤ 0.00), but is a result of pure nanoflare signals, with no presence of strictly wave-based signatures. An explanation could be that nanoflares are individually low-energy events, but they occur very frequently all over the surface of a star. They may come together with an additive effect to form (quasi-)periodic signals, as opposed to the breaking effect of a wave cascade. This cascade-like signal has been documented previously by Hudson (1991), wherein solar nanoflare simulations produced a similar power spectrum cascade, but here we present the first evidence in stellar-specific simulations. This cascade signal is also similar to the 'inverse magnetic cascade' process discussed in Christensson et al. (2001), who found a reverse turbulence effect in 3-D MHD simulations, lending support to an inverse cascade signal generated by magnetic behaviour.  (upper-right), and the percentage of nanoflare power above the noise floor in the range of 1 − 5 mHz (lower-right), displayed as a function of the power-law index, α, and the decay timescale, τ , used to generate the synthetic time series. The observational PSD characteristics (see Table 3 and Figure 3) compare to the modeled PSDs in the range of α = 3.3 ± 0.2 and τ = 200 ± 100 s. Furthermore, the primary frequency (lower-left panel of Figure 7) is sensitive to the nanoflare decay timescale, rising from ∼0.7 mHz at the longest e-folding times (≈500 s) to ∼1.4 mHz at the most rapid decay timescales (≈10 s). Interestingly, once the PSDs have been detrended by their corresponding spectral slopes, the dominant frequencies (upperright panel of Figure 7) present are within the range of 2.7 − 4.2 mHz. This frequency range is often synonymous with the presence of p-mode waves (Andrews 1989(Andrews , 1990a, even though our simulations contain no strict wave activity. An interesting metric to benchmark how significant the power peaks are within the range of 1 − 5 mHz involves the calculation of the percentage of the nanoflare spectral power equal to, or greater, than the corresponding power found in the synthetic noise PSD (lower-right panel of Figure 7). We find that the spectral power arising from strictly nanoflare signatures is 10 − 100% greater than the corresponding (flat) noise power arising from a Poisson-based shot noise distribution. This can be seen in the upper panel of Figure 6, whereby the power arising from nanoflare signals is above that corresponding to the noise floor.
Comparing the simulated PSD features to the M-type stars (see Table 3), we find overlaps with the two-dimensional maps shown in Figure 7 for a power-law index α = 3.3 ± 0.2 and a nanoflare decay timescale τ = 200 ± 100 s. These values are consistent with the first set (α = 3.25 ± 0.15 and τ = 200 ± 100 s) of plasma conditions extracted from the intensity fluctuation statistical distributions. Importantly, we do not find self-similar PSD results substantiating the second set (α = 2.00 ± 0.15 and τ = 200 ± 100 s) of plasma conditions extracted from the intensity fluctuation statistical distributions. This demonstrates the usefulness of employing both statistical and Fourier-based benchmarking of the observational and synthetic time series, since it has allowed us to alleviate a potentially ambiguous result found using just a single analysis method.
M-type stars are nearly or fully convective, with more powerful magnetic activity than the Sun, leading to increased flare activity. However, this alone cannot explain a higher power-law index, as a general boost to activity levels would enhance all frequencies and energies, thus preserving the same power-law index. Instead, it is possible that smallscale nanoflare energies in the range 10 22 − 10 25 ergs are boosted disproportionately in these flare active stars. While low energy flares are likely to be governed by the same underlying physical processes (Lu & Hamilton 1991), and the power-law relationship is scale-free (applying to both small and large flares; Aschwanden 2019), Robinson et al. (1995) and Vlahos et al. (1995) suggested that a discontinuity in the power-law indices of high and low-energy flare events would be an inherent feature of the self-organized criticality model of flaring (wherein small magnetic reconnections occur very frequently, each with the potential to setoff another reconnection nearby, causing an avalanche effect, and following a power-law distribution of energies). They suggest that while high energy flaring would exist at power-law indices of α = 1.8, the power-law index of low energy (i.e. micro and nanoflares) would range around 3 ≤ α ≤ 4. This is in agreement with our stellar nanoflaring power-law index of α = 3.25 ± 0.20. Another explanation for this enhanced rate of small-scale flare activity in M-dwarf stars could lie in the reconnection process itself. Tsuneta & Katsukawa (2004) suggested that low energy (pico-/nano-)flares may occur more favorably via Sweet-Parker reconnection (instead of Petschek processes). If such flare stars have lower Lundquist numbers (i.e., higher plasma resistivity) with respect to the Sun, then this may help explain the enhanced nanoflare rates found in our present study. The mostly convective atmosphere of these flare stars may be able to modify the underlying Lundquist number, allowing for enhanced low-energy nanoflare rates via Sweet-Parker reconnection, but not modifying the rates of the higher energy events that will proceed (as normal) via Petschek reconnection processes. This enhanced nanoflaring may also be linked to the dynamo in these stars. The M-stars in this study sit on the boundary of fully convective atmospheres. While the spectral sub-type where full convection begins is still under debate, estimates are in the range M3 and above (Wright & Drake 2016) to more recent studies suggesting M2.1 to M2.3 (Mullan & Houdebine 2020). Fully convective stars lack the tachocline between convective and radiative zones which powers the solar dynamo . A dynamo powered by helical turbulence is believed to operate in these fully convective stars (Durney et al. 1993;Browning 2008;Pipin, V. V. & Seehafer, N. 2009). This may operate in tandem with the enhanced Sweet-Parker reconnection, through altering the Lundquist number. Investigating the power-law indices of nanoflaring signatures for stars either side of this convective boundary (i.e., M1 and M5) would allow us to test this theory in a future study.

CONCLUSIONS
We have employed a combination of statistical and Fourier-based analysis techniques to search for evidence of nanoflare activity in M-type stars observed by NGTS. The intensity fluctuation distributions of the M-type stars revealed both negative median offsets and positive Fisher skewness values, highlighting the presence of impulsive intensity rises, followed by exponential decays, trapped within the noise envelop of their corresponding lightcurves. To validate these signatures, we examined complementary A-type non-flare active and K-type low-activity stars, which demonstrated zero median offsets, alongside very minor Fisher skewness values, highlighting the more symmetric composition of these A-and K-type time series that are devoid of nanoflare signatures.
Previous studies have observed periodic phenomena in M-type stars that has been interpreted as evidence of pmode wave activity. To investigate whether nanoflare signatures, which are governed by a power-law index, may contribute to similar (quasi-)periodicities we calculated power spectral densities (PSDs) of the NGTS time series. Longduration successively acquired time series (2316 individual datapoints) were employed to maximize the frequency resolution. We found contrasting spectral features between the A-K-and M-type time series. The A-type spectra had flat power trends representative of pure shot noise distributions, The Ktype PSD was flat apart from frequency enhancements across the range of ≈ 1−10 mHz, indicative of p-mode wave signatures, as we would expect from a Solar-like K-type star . By contrast, the M-type spectra revealed spectral slopes and frequency enhancements across the range of ≈ 1 − 10 mHz. As a result, it was unclear whether the frequency enhancements were a result of nanoflare activity governed by a power-law relationship, or the capture of p-mode wave signatures. To investigate this further, we employed Monte Carlo models of nanoflare activity to examine whether pure flare signatures have the ability to manifest as spectral power enhancements in their corresponding PSDs.
A grid of 6100 Monte Carlo models were constructed that replicate the exposure time, cadence, and duration of the NGTS observations, but with each resulting time series generated from different combinations of power-law indices, α, and flare decay timescales, τ . Each of the time series were added to synthetic shot noise distributions to simulate realistic NGTS lightcurves. These were examined using identical statistical and Fourier-based techniques, with the results cross-correlated to the observational findings. Importantly, we found evidence that time series comprised of nothing but impulsive nanoflare signatures and Poisson-based shot noise are able to demonstrate spectral power peaks across the frequency range ≈ 1 − 10 mHz, suggesting that previously detected p-mode signatures may actually arise from nanoflare activity in the host star. Combining both the statistical and PSD benchmarks, we find evidence for stellar nanoflare activity across the sampled M-type stars for a power-law index α = 3.25 ± 0.20 and a decay timescale τ = 200 ± 100 s.
Looking to the future, higher cadence observations from instruments such as HiPERCAM (Dhillon et al. 2016) may allow for the more rapid accumulation of suitable number statistics, plus the ability to investigate potential nanoflare signals across a number of different color photometry bands.As the nanoflaring parameters observed in the Sun by Jess et al. (2019) varied with wavelength, we would expect a similar result for stellar nanoflares. Such multi-wavelength observations could allow for a limited analysis of how the nanoflare signals differ throughout the stellar atmosphere. Investigating young Sun-like stars with HiPERCAM is also an area of interest, as the high cadence and multi-color observation can overcome the difficulties of increased flare contrast on these young and active stars. Due to their highly active x-ray emission and coronal temperatures (Johnstone & Güdel 2015), we expect a very high degree of nanoflare activity, possibly leading to stellar coronal heating via nanoflaring. Furthermore, a follow up observational campaign could leverage the large sky sampling of the NGTS to examine the presence of nanoflare signatures on other spectral classifications, particularly M1 and M5 spectral types, and investigate how the convective boundary affects the nanoflare power-law indices. This larger star sample could also investigate the source of the dip below the idealized Gaussian at approximately −0.90 σ N in the statistical distribution of the M-stars. To further improve the Fourier-based PSD analyses, we propose a more continuous observational platform that will further increase the frequency resolution possible, e.g. the Tran-siting Exoplanet Survey Satellite (TESS; Ricker et al. 2014), which can operate in both 240 s and 20 s cadences. The 20 s cadence data is part of the extended mission program that will begin operations in July 2020. The obvious advantages of space-based observations would allow us to minimize any high frequency (scintillation) noise present in the stellar lightcurves, while also allowing for a much higher frequency resolution in the subsequent PSD analyses. Longerduration observations have been proposed to study stellar oscillations in greater detail (Ball et al. 2018), and this capability would extend the same advantages to our nanoflare PSD analyses.  Table A4.    Table A4. The Spectral type, NGTS identifier, Gaia source ID, Tess Input Catalog (TIC) ID, RA, Dec, Stellar Mass (in Solar mass units), Stellar Radius (in Solar radi units), Stellar Luminosity (in Solar luminosity units), Distance (in parsecs), Macroscopic Flare Rate (per hour) and the ratio log Lx L Bol for the stars used in the analysis. The Stellar masses, radi. and luminosity data is from the Tess Input Catalog release V8. (Stassun et al. 2018) The ratio log Lx L Bol is calculated by comparing the log ratio of the x-ray luminosity as observed by 4XMM XMM-Newton Serendipitous Source Catalog (Webb et al. 2020), with the stars luminosity. A ratio of log Lx L Bol ∼ −3 is expected for x-ray saturated, young and active M stars (Kastner et al. 2003) . Figure B8 shows the trendlines (calculated over ±6 frequency elements or ±0.427 mHz) of Fourier power spectral densities (PSDs) for the example A K, and M stars, as well as for a modeled time series with a power-law index α = 3.25 and a flare decay timescale τ = 245 ± 24.5 s. This plot highlights the agreement in the observational M-type and modeled time series PSDs, with comparable spectral slopes of approximately β = −0.30 ± 0.05, and peaks around ≈0.8 mHz . This is in contrast with the Aand K-type PSD, which are relatively flat and featureless by comparison. Figure B8. The Fourier power spectral density (PSD) trendlines calculated over ±6 frequency elements (±0.427 mHz) for example A-K-and M-type stellar sources NGTS J025840.5-120246 (red line), NGTS J030000.7-105633 (orange line) and NGTS J030047.1-113651 (blue line), and a modeled time series corresponding to a power-law index α = 3.25 and a flare decay timescale τ = 245 ± 24.5 s (black line), displayed in normalized units of σ 2 N /mHz. It can be seen that the A-type and K-type spectra are relatively flat across all frequencies with no evidence of distinct peak frequencies. The K-type does show some slight power enhancement between ≈ 1 − 10 mHz, consistent with stellar p-mode oscillations, as have been previously observed in K-type Solar-like stars (e.g. Chaplin et al. 2009) . Contrarily, the M-type PSD has a primary power peak at ≈0.8 mHz, followed by decreasing spectral power exhibiting a spectral slope of β = −0.30 ± 0.05, followed by numerous power peaks in the range of 3 − 10 mHz, consistent with previous links to stellar p-mode spectra. The synthetic PSD is remarkably similar to the NGTS M-type stellar source, with peaks and spectral slopes in the same range and magnitude (see Figure 7 for the full range of peak frequencies and spectral slopes in modeled time series PSDs.

C. OBSERVATIONAL CONSIDERATIONS
The detectability of nanoflare signals via statistical and periodic analyses are dependent on the underlying observational parameters, including, • Time Series Length: The statistical analysis is dependent on the number of frames, N . The error in statistical analyses scales with √ N , while the signal scales with N . Periodic analysis is also dependent on the time series length, but crucially on the length of successive uninterrupted frames. Increasing the duration of the observations will provide increased fre-quency resolution. The periodic signal also benefits from increased number statistics, as the number of nanoflares captured increases with longer observing sequences, hence providing more accurate quantification of any associated periodicities. We have investigated modeled lightcurves (which are not subject to day/night cycles), and found increasing the number of successive frames had the effect of increasing the ratio of nanoflare power above the noise floor in the range 1 − 5 mHz (i.e., the lower panel of Figure 7). As a result, the nanoflare periodic signatures became more prominent over the noise. We expect space-based (e.g., TESS) observation to allow us to uncover more of the underlying spectral slopes, particularly for the highest power-law values.
• Cadence: Shorter cadences will allow for increased Nyquist frequencies to better resolve rapid and short-lived periodic signatures. Sub-second cadences (e.g., HiPERCAM, with exposures on the order of milliseconds; Dhillon et al. 2016) could allow for a very large frequency range and excellent number statistics to be achieved in a very short observation window.
• Apparent Magnitude: As the observed magnitude increases, the scintillation noise begins to increase also (this is not an issue with space-based observations). This would affect the frequency distribution of the noise, since the scintillation introduces a frequency-dependent noise component that needs to be considered. Searching for nanoflares embedded within this more complex noise distribution would require the seeding of a scintillation model into the numerical simulations. A future study could explore high-magnitude stars, to determine whether the increased scintillation is balanced by the increased nanoflare signal, or future space observations (e.g., TESS) could mitigate this entirely. However, the star itself should still be of low intrinsic stellar brightness; see below.
• Intrinsic Stellar Brightness: Brighter stars have increased quiescent flux, and therefore a more pronounced noise floor that must be combated when searching for nanoflare signals on top of this brighter background. This means the contrast between the nanoflare signals and the background becomes a challenging issue. Even at solar-like luminosities, the detection of microflare energies becomes difficult, let alone nanoflares on stellar sources that cannot be spatially resolved.