Statistical Signatures of Nanoflare Activity. III. Evidence of Enhanced Nanoflaring Rates in Fully Convective stars as Observed by the NGTS

Previous examinations of fully convective M-dwarf stars have highlighted enhanced rates of nanoflare activity on these distant stellar sources. However, the specific role the convective boundary, which is believed to be present for spectral types earlier than M2.5V, plays on the observed nanoflare rates is not yet known. Here, we utilize a combination of statistical and Fourier techniques to examine M-dwarf stellar lightcurves that lie on either side of the convective boundary. We find that fully convective M2.5V (and later subtypes) stars have greatly enhanced nanoflare rates compared with their pre-dynamo mode-transition counterparts. Specifically, we derive a flaring power-law index in the region of 3.00 ± 0.20, alongside a decay timescale of 200 ± 100 s for M2.5V and M3V stars, matching those seen in prior observations of similar stellar subtypes. Interestingly, M4V stars exhibit longer decay timescales of 450 ± 50 s, along with an increased power-law index of 3.10 ± 0.18, suggesting an interplay between the rate of nanoflare occurrence and the intrinsic plasma parameters, e.g., the underlying Lundquist number. In contrast, partially convective (i.e., earlier subtypes from M0V to M2V) M-dwarf stars exhibit very weak nanoflare activity, which is not easily identifiable using statistical or Fourier techniques. This suggests that fully convective stellar atmospheres favor small-scale magnetic reconnection, leading to implications for the flare-energy budgets of these stars. Understanding why small-scale reconnection is enhanced in fully convective atmospheres may help solve questions relating to the dynamo behavior of these stellar sources.


INTRODUCTION
Magnetic reconnection is a fundamental physical process in conducting plasmas that allows for the conversion of magnetic energy through the rearrangement of magnetic fields (e.g., Priest 1986;Reale 2007;Cargill et al. 2015).Magnetic reconnection came to initial prominence in astronomy Corresponding author: S. D. T. Grant samuel.grant@qub.ac.uk as a proposed mechanism for observed solar flare activity.The first derived reconnection model was presented by Sweet (1958) and Parker (1957) as a two-dimensional magnetohydrodynamic (MHD) configuration, where there exists a long, thin diffusion region enabling magnetic fields to reconnect.The reconnection rate of this Sweet-Parker model was found to be inversely proportional to the square root of the dimensionless Lundquist number, S, which is defined as; where L is the length of the diffusion region, v A is the Alfvén speed, and η is the plasma resistivity (a full derivation is presented by Priest & Forbes 2007).Predicted Sweet-Parker reconnection rates for solar conditions could not recreate actual observations and the derived energetics from Sweet-Parker events were orders-of-magnitude smaller than what was observed in the solar corona (e.g., Crosby et al. 1993).An alternative model was presented by Petschek (1964), which allowed for faster rates by permitting reconnection across far shorter length scales of the diffusion region (Priest & Forbes 2007).This remedies an issue with the Sweet-Parker model (reconnection rate ∝ 1/ √ S) due to the Petschek reconnection rate being inversely proportional to the logarithm of the Lundquist number, which limits the influence of plasma conductivity and provides more robust similarities to the characteristics of large flare events (e.g., Aschwanden 2020).
Observationally, magnetic reconnection is characterised by an impulsive brightening as magnetic energy is converted into localised plasma heating, and is classically seen as stochastic, macroscopic events (see the reviews of Cargill & Klimchuk 2004;Benz & Güdel 2010;Fletcher et al. 2011;Benz 2017).Subsequently, the plasma gradually cools over an extended period, which manifests as an exponentially decaying intensity from the time of maximum brightness (e.g., Moffett 1974;Moffett & Bopp 1976;Kowalski et al. 2013;Pitkin et al. 2014).The quantification of this decay process is measured through the e-folding time, τ , which is the time taken for the flare luminosity to decrease by a factor 1/e.The magnitude of this value is dependent on the underlying local plasma conditions, such as the efficiencies of evaporative, non-evaporative, conductive, and radiative cooling processes (Antiochos & Sturrock 1978).Solar and stellar flare energies are governed by a power-law relationship (Aschwanden et al. 2000).Here, the power-law exponent governs the frequency, dN/dE, of flaring events with an associated energy, E, through the relationship, where α represents the power-law index.The nature of a power-law relation dictates that low-energy flares will be many times more frequent than larger events, and that smallscale events become more energetically important as the power-law index, α, increases.A range of power-law indices have been documented across varying solar and stellar flare energy windows, from 1.35 ≤ α ≤ 2.90 (Berghmans et al. 1998;Krucker & Benz 1998;Aschwanden 1999;Parnell & Jupp 2000;Benz & Krucker 2002;Winebarger et al. 2002;Aschwanden & Freeland 2012;Aschwanden et al. 2014Aschwanden et al. , 2015)).
The stochastic nature of large solar flaring eventsindicates that they are too infrequent to be viable heating mechanism for the extra-ordinary temperatures of the solar corona, known as the coronal heating paradox.Instead, nanoflares, with individual energies around 10 9 times less than their large-scale counterparts, were proposed as an alternative due to their higher occurrence rates (Parker 1988).In order to be considered as consequential to the flare energy budget, and thus atmospheric heating, it has been established that the minimum requirement is α ≥ 2 (Parker 1988;Hudson 1991).Due to their individual low energies, nanoflares are typically embedded within the noise floor of the measured intensity signals, leading to difficulties identifying individual nanoflare events.However, their higher occurrence frequency means they can be recovered from time series data using statistical techniques that do not rely on the individual identification of macroscopic intensity signals.
Building on the work of Terzo et al. (2011) and Jess et al. (2014), Jess et al. (2019, henceforth referred to as Paper I) developed a robust method for nanoflare investigation.Through Monte-Carlo simulations, realistic nanoflare lightcurves were generated for a wide range of α and τ values, coupled with precise modeling of the noise characteristics of solar observables.Through this, Paper I was able to uncover nanoflare signatures in solar coronal observations, manifesting as asymmetric contributions to the intensity fluctuation distributions of coronal images, consistent with power-law distributions on the order of 1.82 ≤ α ≤ 1.90.Despite showing that the solar active region under study did not appear to contain the necessary nanoflare activity to influence coronal heating, Paper I provided a comprehensive method for analysing small-scale flare activity in intensity time series, and suggested the same techniques could be applied to stellar observations.Subsequently, Dillon et al. (2020, henceforth referred to as Paper II) utilised the techniques of Paper I on stellar lightcurves from A, K and M-type stars to investigate whether signals previously interpreted as p-mode oscillations in dMe flare stars could in-fact be caused by nanoflares.Fourier analysis of each spectral class showed no enhanced power around p-mode frequencies (i.e., 1 − 1000 s) in Atypes, power enhancement at p-mode frequencies in the Ktype, and power enhancements across the entire frequency spectrum in M-dwarfs.These enhancements were classically seen as evidence of global wave activity generated in both the K and M-type stars, however Monte-Carlo simulations revealed that the M-type stars produced the asymmetric intensity fluctuation distribution effects consistent with nanoflares, as opposed to the symmetrical effects observed in the K-type distributions that are consistent with dominant oscillatory behaviour.The M-dwarf flare activity produced a power-law index of α = 3.25 ± 0.20, greater than previously reported values, and enabled Paper II to show that the flaring rate was high enough for nanoflare signals to appear quasi-periodic when an entire stellar disk is integrated into a single lightcurve, thus explaining their influence on the resulting Fourier power spectra.The reason for the enhanced nanoflare power-law indices for M dwarfs was not known, but it was theorized that the fully convective nature of these stars may be responsible.
While solar-like stars have a combination of convective and radiative zones bridging their core and visible surface, some stars operate in a fully convective manner.The change from partially-to fully-convective interiors has been related to the 'convective boundary', distinguished by a lack of tachocline in later stars that are fully convective.The tachocline is a thin region of the stellar interior at the boundary between the radiative and convective zones that contains large radial shears due to the imbalance between the rigid radiative zone and the differentially rotating outer convective zone (Spiegel & Zahn 1992;Browning 2008).Wright & Drake (2016) estimated that this transition occurs in M-dwarf stars around M3V and later, with recent studies suggesting a more precise transition at approximately M2.1 − 2.3V (Mullan & Houdebine 2020).
Convection is a primary driver of magnetic reconnection in stars (Pedersen et al. 2017).As magnetic reconnection is the driving force behind flares, changes to the convective nature of a star have important implications for the resulting flare dynamics.The tachocline is thought to play a role in strengthening magnetic fields in partially-convective stars such as the Sun, as the shear forces across the region can convert poloidal fields into stronger toroidal configurations (Parfrey & Menou 2007).However, it is important to note that the tachocline is not necessary for magnetic field generation, since fully convective stars also exhibit magnetism, where the dynamo is theorized to be driven by helical turbulence (Durney et al. 1993;Browning 2008;Pipin & Seehafer 2009), but this change in dynamo is still under debate.Indeed, Wright & Drake (2016) and Wright et al. (2018) investigated the relationship between stellar rotation and activity levels for fully-convective late M-type dwarf stars.They found that the rotation/activity relationship for fully convective stars was almost indistinguishable from partially convective stars, suggesting the solar-type dynamo is independent of the presence of a tachocline.
Returning to the results of Durney et al. (1993), Browning (2008) and Pipin & Seehafer (2009), it is possible to hypothesize that the enhancement of nanoflaring rates is linked to the conversion of late M-type interiors to fully-convective, and hence to the consequent changes induced in the helical dynamo processes.Previous examinations of stellar flares on late-type MV stars have found a range of power-law indices.Early space-based observations produced values of α ∼ 1.5 (Collura et al. 1988), although a range of partiallyand fully-convective MV stars were included in the sample studied.Subsequent studies of late-type stars produced in-dices α > 2 and showed that there were no discrepancies between ground-and space-based observatories in these calculations (Robinson et al. 1995(Robinson et al. , 1999)).A trend developed where reported power-law indices increased as the complexity of the methods for isolating marginal flare signals developed (e.g., Güdel et al. 2003;Güdel 2004;Welsh et al. 2006;Hawley et al. 2014a), with values as high as α = 2.7.However, this remains below the nanoflare power-law index (α ≈ 3.25) presented in Paper II.This may be a result of the novel detection techniques of Paper I providing unprecedented access to the lower energy spectrum of flares, as it has been seen in solar observations that the power-law index, α, can change depending on the size and energetics of the flares under consideration in the sample (e.g., Wang & Dai 2013;Ryan et al. 2016;Milligan et al. 2020).Therefore, there may exist a similar noticable shift in the gradient of the power-law index around the transition between large and small-scale flares in MV stars.Hence, it is important to examine the power-law indices associated with nanoflare activity across a wide range of MV spectral sub-types with the techniques described in Paper I and Paper II to better understand the influence of fully-convective interiors in the generation of nanoflares.
As discussed in Paper II, an alternative source for the enhanced rate of small-scale reconnection may be the fullyconvective stars having plasma with a higher resistivity value (Mohanty et al. 2002), which lowers the associated plasma Lundquist numbers.Small-scale flaring has been shown to occur more favorably via Sweet-Parker reconnection (Tsuneta & Katsukawa 2004), thus enhanced nanoflaring can be expected in stars with low plasma Lundquist numbers.If nanoflare rates are enhanced in fully convective stars, then investigating whether this is due to the change in dynamo, or down to the plasma resistivity, could answer important questions regarding the dynamo physics in operation in these stars.

BACKGROUND TO PREVIOUS STATISTICAL STELLAR NANOFLARE ANALYSIS
Paper II used a combination of statistical and Fourier analyses to investigate nanoflare populations in a variety of stellar sources.Through comparisons between observational lightcurves and synthetic time series with simulated nanoflare signals, enhanced nanoflare activity in late-type MV stars was observed, through a larger power-law index when the nanoflare occurrence rates were plotted as a function of their underlying energy.In order to determine why the power-law indices in these stars were significantly larger than in other stellar and solar studies, we re-apply the anal- the change in the statistical and Fourier parameters associated with the embedded nanoflare conditions) either side of the fully-convective boundary, we aim to identify the role a fully-convective atmosphere plays in nanoflare occurrence rates.This work constitutes the third contribution to a series of studies on nanoflare behavior in solar and stellar atmospheres.Given the consistency of techniques applied across these studies, it is prudent to provide contextual information on previous works.For detailed discussions of the analysis procedures and modeling set-ups, it is advised to consult Paper I and Paper II directly.
The initial detection method for nanoflare signatures involves the statistical analyses of quiescent intensity fluctuations following a traditional Z-scores approach (Sprinthall 1990).A histogram of the fluctuations is generated, with two distinct statistical deviations away from a standardized Gaussian distribution providing evidence of nanoflares.These signatures, identified by Terzo et al. (2011) and Jess et al. (2014), were diagnosed through the modeling of the observed lightcurves by embedding the rapid rises and exponential decays of intensities around the noise limit that are associated with the energetics of nanoflares.The first detection characteristic is a negative median offset of the intensity fluctuation distribution, whereby the median value of the histogram is < 0 σ N , i.e., offset from the mean of the distribution that is equal to 0 σ N following the application of polynomial detrending.This was shown to be associated with the exponentially decaying nature of the nanoflare lightcurve.The decay phase produces more significant negative fluctuations below the mean than the relatively brief elevated signal of the flare event, thus providing an offset between the median and mean that is directly measurable.The second signature of nanoflare activity is an excess of fluctuations at ∼ 2 σ N , which is caused by the associated energetics of nanoflares producing consistent peak brightenings around ∼ 2 σ N (see Jess et al. 2014 and Paper I for further details on flare energy modeling).This gives rise to an asymmetric distribution with a slight excess of fluctuations visible at ∼ 2 σ N in the corresponding histogram, which is best categorized through Fisher skewness coefficients.In addition to the primary nanoflare indicators described above, benchmarks on the shapes and widths of the intensity fluctuation distributions are provided through calculation of the kurtosis and ζ values, where ζ is the ratio of the full-width at eighth-maximum to that of the full-width at half-maximum (i.e., FW 1 8 M-to-FWHM ratio) of the resulting distribution.These parameters allow for additional information on the nature of the distribution, with kurtosis linked to the prevalence of outliers at high σ N values, and deviations from the established ζ value for a normal distribution (ζ = 1.73) revealing non-linear distributions of measurements around the mean.
As established in Paper I & Paper II, a quiescent lightcurve exhibiting both of these statistical signals, identified through the median offset and Fisher skewness, contains embedded nanoflare signatures.Terzo et al. (2011) presented an observation of nanoflare activity that only induced a median offset in the distribution, due to the weak nanoflare signal not producing sufficient peak amplitudes to influence the positive wing of the distribution.However, both effects presenting in observations, and corroborated by the four chief diagnostics, cannot be explained by any alternative mechanism.In an observation of purely ambient solar or stellar plasma, noise would be the only artifact.However, noise fluctuations follow a standard Gaussian distribution as a result of Poisson statistics tending to a Gaussian in the limit of large number statistics (Terzo et al. 2011).Solar and stellar observations contain a variety of oscillatory phenomena, ranging from the signatures of MHD waves in their atmospheres (see the review of Jess et al. 2023) to modulation caused by the rotation of the star and starspots.In the context of this study, these signals have no influence over the asymmetric effects under consideration, as the sinusoidal nature of linear oscillations produces a symmetric distribution around the mean.Starspots, characterized by periodic reductions in intensity for times related to rotation, will decrease both the median and mean of the sample, providing a minimal effect in the offset between the two.Their effect on the timeseries is also mitigated through the application of polynomial detrending, as the periodicity of their effect is on the order of multiple days.MHD waves are also capable of steepening into non-linear shocks that produce notable intensity enhancements over a short period (e.g., Carlsson & Stein 1997;Grant et al. 2018).However, these shocks display a 'saw-tooth' pattern to their intensity morphology, in contrast to the exponential decay of flares, and thus would not reproduce the statistical effects described above.Therefore, with these effects discounted, and any macroscopic flare events robustly removed from the data (see Section 3), there is a confidence that the statistical effects derived from the intensity fluctuation distributions are positive signatures of nanoflare activity.
Through inspection of nanoflare observations with insight from simulations, Terzo et al. (2011) and Jess et al. (2014) found the interval between nanoflares in a single lightcurve to be ∼ 360 s, a similar frequency to ubiquitous p-mode signatures (∼ 1 − 10 mHz; Andrews 1989; Rodríguez-López et al. 2014;Rodríguez et al. 2016).Subsequently, Paper I were able to derive the frequency of flaring events given their power-law index, showing that the greater the power-law index, the higher the frequency of nanoflares in a given sample (see Figure 4 of Paper I).With the low flaring intervals and high frequencies found, it was suggested that nanoflare signals, when integrated across a field-of-view, can no longer be viewed as stochastic events like their macroscopic counterparts.Instead, they can be considered to be a quasi-periodic phenomena, particularly given the exceptionally large power law index of α ∼ 3.25 reported by Paper II.
The quasi-periodic nature of nanoflares was successfully utilized in Paper II to further verify their detection.When inspecting the statistical distributions of K-and M-type stars, there was the potential, but marginal, signature of nanoflares in K-types, as opposed to clear nanoflare signals in the M-types.Paper II subsequently employed Fourier techniques to distinguish the quasi-periodicities of stellar nanoflares.Power spectral densities (PSDs) were computed for the longest continuous time series common across all stars in order to maximize the resulting frequency resolution.These PSDs revealed power at p-mode frequencies in both types, however the nature of that power differed, with M-types showing power enhancement across a wide range of frequencies, synonymous with nanoflare activity from a range of interval times/frequencies, as opposed to the strictly oscillatory nature of K-type atmospheres (see Figure 3 of Paper II).Additionally, the PSDs in M-types displayed a prominent spectral slope following the peak energy value, which was also found to be due to the underlying nanoflare signal.Fourier analysis is therefore a valuable tool in both confirming nanoflare signal in the source, and differentiating between those atmospheres with wave and nanoflare interplay (i.e., partially-convective interiors) and where nanoflares dominate the energy landscape (i.e., fullyconvective interiors).
Throughout these works, the observational findings have been substantiated and interpreted through simulated nanoflare time series.Paper I devised the methodology and analysis techniques, based on the Monte-Carlo modeling of flare intensities with additionally added camera-specific noise signatures, coupled with a range of typical nanoflare amplitudes and decay rates, and provides a detailed description of the set-up.Paper II took these techniques and tailored them for a strictly stellar scenario, including the remodeling of the noise profiles and resizing the observed area to mimic a full stellar disk, as opposed to a sub-section of the solar atmosphere used in Paper I. Flare energies between 10 22 − 10 25 erg, typical of nanoflares, were used alongside a linear scaling relationship to reproduce the observed counts.Each simulated lightcurve has two variables controlling the nanoflare input: the power-law index, α, and the e-folding time, τ .Lightcurves were then generated for a dense grid of parameter ranges, 1 ≤ α ≤ 4 (in steps of 0.05) & 5 ≤ τ ≤ 500 (in steps of 5 s), consistent with previous viable observation ranges (Terzo et al. 2011;Jess et al. 2014).This synthesis produced 6100 lightcurves embedded with characteristic noise and unique nanoflare configurations.Subsequent comparison between synthetic and observed lightcurves is achieved through forming distributions of each synthetic time-series, and identifying the configuration with matching median offset, Fisher skewness, kurtosis and ζ values (see Figure 5 of Paper I).Given the stringent requirements of this matching criteria, unique solutions are found for the observables in both Paper I and Paper II.Given the computational demands, the simulations presented in Paper II are used in this study, thus further details on their reproducability can be found in the original paper.
In this paper, we apply proven statistical nanoflare analysis techniques to a wide range of M-type stars that lie on either side of the predicted fully-convective boundary.We compare their statistical and Fourier properties to the simulations generated in Paper II in order to determine the probable underlying nanoflare conditions and the effect of the convective boundary on the uncovered nanoflare properties.

OBSERVATIONS WITH NGTS
To remain consistent with Paper II, the Next Generation Transit Survey (NGTS; Wheatley et al. 2018) was utilized to obtain the observations.The long time series (each in excess of 10 5 frames) and short cadence (∼ 12 s) available for thousands of M-type stars allow for the accumulation of suitable number statistics necessary for nanoflare analyses.The initial spectral classification generated by the NGTS pipeline (which utilizes Spectral Energy Distribution fitting, see section 5.1.1 in Wheatley et al. 2018) was combined with stellar parameters from the TESS Input Catalog Version 8 (TIC V8; Stassun et al. 2018), to ensure robust spectral sub-type identification.
To ensure the ideal sample of objects for study, a number of selection criteria were applied to the catalog extracted from NGTS to remove unwanted artifacts.Initially, the magnitude of the stars were constrained to ensure similar photon noise characteristics for each object.Thus, only stars with magnitudes matching the range of the previous study of Paper II (spanning NGTS magnitudes of ∼ 12 − 14) were progressed.This ensured that the magnitude of the fluctuations in each spectral type were approximately equivalent, with M0V and M4V stars exhibiting standard deviations of 1.6 ± 0.1% and 2.0 ± 0.6%, respectively.The average standard deviation of the time-series across the full sample was 2.0 ± 0.8% of the mean, thus an equivalent fluctuation profile can be applied to the modeling.Next, in a similar manner to Jackson et al. (2023), complementary photometric and astrometric data from Gaia DR3 (Gaia Collaboration et al. 2016, 2022) were utilized to exclude unwanted candidates from the sample.In particular, astrometric excess noise analysis was applied to exclude any binary systems (Evans 2018), and the photometric filtering processes of Arenou et al. (2018) identified any blended sources.A final step was to consider the rotation rate of the stars in the sample.Given the length of the NGTS time series, it is not possible to extract the periods associated with slow rotating M-dwarfs (i.e., above 30 days).Therefore, generalized Lomb-Scargle techniques (Lomb 1976;Scargle 1982) were applied to the NGTS data to identify fast rotators and exclude them from the sample.This is done to ensure the sample contains similarly 'slow' rotators, as defined in previous studies (e.g., Mondrik et al. 2019), thus removing the influence rotation rate has on increasing stellar activity (West et al. 2008;Candelaresi et al. 2014).Only two candidates exhibited definable periodicities and were thus excluded, an M2V & an M4V with periods of 21.1 days and 19.8 days, respectively.
After accounting for magnitude considerations, avoiding blended sources, ensuring TIC matching and excluding fast rotators, we were able to find 5 stars for each spectral subtype, consisting of M0V, M1V, M2V, M2.5V, M3V, and M4V.The stellar properties (NGTS identifier, RA/Dec, magnitude, etc.) of these candidates are provided in Table A1.Only one suitable M5V star with TIC-derived stellar parameters could be identified, and no sub-types later than this were found.The intrinsic brightness of M-dwarfs decreases with increasing sub-type (Yang et al. 2017), leading to difficulty in identifying suitable candidate stars with the desired brightness properties.Future investigations of post-M4V stars may be fruitful, but identifying a suitable number of candidates may prove difficult with existing instrumentation.Hence, we limit our current study within the range of M0V -M4V, where we have multiple candidates available for compari- A standardized Gaussian profile is overplotted in each panel using a dashed red line for reference.The M4V-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 M0Vtype intensity fluctuations provide effectively zero negative median offset, and no elevated occurrences at ∼ 2 σN .This is inconsistent with clear statistical signatures of nanoflare activity, with the resulting distribution remaining more 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 negative median offsets and occurrence excesses, respectively, are clearly visible for the M4V stellar source.For improved clarity, the blue and gold lines display the corresponding distributions in each zoomed panel.
son.This range also overlaps well with the predicted dynamo mode transition to fully convective(M2.1 − 2.3V; Mullan & Houdebine 2020), making it suitable for the study of the role that fully-convective starsnplay in the resulting nanoflare activity.
The lightcurves were background corrected and flat-fielded via the NGTS data reduction pipeline described and visualized in Wheatley et al. (2018).This pipeline calculates a relative error in the flux at each data point in the time series.This error correlates with cloudy weather and/or high airmass values.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.
To prepare the data for statistical analysis, each lightcurve was detrended by a low-order polynomial so the mean value is zero.Then, the time series is subsequently renormalized by its own standard deviation, σ N .Next, the lightcurves extracted for each observing sequence were examined for the presence of non-Gaussian intensity enhancements such as macroscopic flare signatures following the methodology described by Paper II. Emission signatures exceeding 3σ N above the mean value, lasting continually for a minimum of 1 minute (5 datapoints), were identified in each lightcurve.Based on a normal distribution, the probability of these event presenting through Gaussian-Poisson noise is ≲ 2 × 10 −13 , hence allowing for robust detection of macroscopic flaring activity.Every star, apart from the M2V candidate NGTS J062005.7-372555,demonstrated macroscopic flare signatures, resulting in the removal of a further ∼ 0.2 − 2.5% of the remaining M-type time series.The degree of macroscopic flare emission varied with the spectral sub-type, with M4V stars exhibiting approximately five times more detected flares than the M0V stellar types, consistent with previous studies (Hawley et al. 2014a;Yang et al. 2017).Once the larger-scale flare signatures had been identified, they were subsequently removed from the time series using an interval of ±5 minutes (±25 datapoints) from the first and last detection above the 3σ N threshold (see Figure 1).Removing these signatures allows for the assumption of normality in the intensity distribution, i.e. shot and readout noise combined with ambient stellar intensity (Terrell 1977;Delouille et al. 2008).This has been shown as valid in Terzo et al.A1.To ensure consistency with previous stellar nanoflare investigations, the filtering steps employed were identical to those used in Paper II, with the filtered lightcurves subsequently cropped to 97 060 datapoints each to match the number statistics from the previous study.This allows a direct comparison to be made with the work of Paper II, since the previously published nanoflare simulations can be re-used due to identical number statistics, filtering techniques, desired α (power-law index) and τ (efolding time) ranges, in addition to specific NGTS-modeled noise characteristics.

ANALYSIS AND DISCUSSION
To investigate the possible changing nanoflare properties with spectral type, we utilized the statistical and Fourier analysis techniques outlined in Section 2. As outlined, nanoflares give rise to 2 distinct statistical signatures, which can be used to diagnose stellar nanoflare activity.We present two example histograms of intensity fluctuations in Figure 2 for stars NGTS J052346.3-361114(M0V spectral type; top panel) and NGTS J050423.8-373021(M4V spectral type; lower panel).From Figure 2 it is clear that opposite ends of the included spectral types, which lie on either side of the predicted fullyconvective boundary, demonstrate distinctly different statistical signatures.The M0V star exhibits weak nanoflare signatures, with a marginal negative median offset and no elevated intensity fluctuations at ∼ 2 σ N .On the contrary, the M4V star has a clear excess of ∼ 2 σ N intensity fluctuations in addition to a prominent negative median offset.The signatures of the M4V star shown in Figure 2 are consistent with previous positive stellar nanoflare identifications in Paper II.The distinct increase of visible nanoflare signatures within the expected regime of full convection indicates that the enhanced nanoflare rates are related to the underlying convective nature of the star.
These examples illustrated in Figure 2 clearly identify the vastly different nanoflare signatures present at either end of the investigated range of spectral sub-types.To better examine the change in nanoflare activity across the given spectral range (M0V -M4V), the derived properties were averaged according to their specific spectral type following the bootstrap method documented by Efron et al. (1979).Straightforward averaging of features that are dependent on the underlying stellar plasma conditions from multiple stars is challenging due to the uncertain behavior of the standard errors of the given parameters.Hence, bootstrapping techniques are used extensively throughout the physical sciences to better calculate confidence intervals for data following non-standard or unknown distributions (Simpson & Mayer-Hasselwander 1986;Desmars et al. 2009;Yao et al. 2017).
Figure 3 shows the change in the median offset, kurtosis, Fisher skewness, and ζ values, respectively, as a function of spectral sub-type.The results are also tabulated in Table 1.From Figure 3, we find a distinct change in the nanoflare statistical signatures as a function of spectral sub-type, suggesting the convective boundary may play an important role in the generation of efficient nanoflare conditions.We find that M2.5V (and beyond) stars exhibit distinct nanoflare statistical signatures that are consistent with those put forward by Paper II.Specifically, the average median offset for the pre-M2.5Vstars exhibits a large spread around a weakly offset value (upper panel of Figure 3), while the post-M2.5Vstars The Fisher skewness value is effectively zero for pre-M2.5Vstars (second panel from bottom in Figure 3), suggesting no, or very weak, nanoflare activity.From M2.5V onward, there is a clear increasing trend in the Fisher skewness value of the fluctuation distribution, with the M4V subtype displaying a Fisher skewness equal to 0.051 ± 0.014, providing strong evidence for the presence of nanoflares.
In the additional distribution diagnostics, the relationship is less clear.Regarding the kurtosis (second panel from top in Figure 3) there appears to be a trend in that statistical kurtosis increases across the full sample, between M0V and M4V.However, the exact nature of this relationship is obscured by the large uncertainty in the kurtosis for spectral types around the dynamo mode transition.In particular, the M3V sub-type appears to exhibit a decrease in kurtosis, though this may be due to the large uncertainty in M2V & M2.5V producing abnormally large values.
There is no clear trend visible in the corresponding ζ values (lower panel of Figure 3 As described in Section 2 (and in detail in Paper I), the products of the four distribution diagnostics can be used to derive the power-law index, α, and the nanoflare decay timescale, τ , for each observation.Through the calculation of median offset, Fisher skewness, kurtosis, and ζ for each of the 6100 synthesized lightcurves, these can be compared to the values seen derived in Table 1.This was achieved by considering each diagnostic individually.The values of the statistical diagnostic within a range of ±1σ N are directly N /mHz.The crosses in each panel depict the individual power values as a function of frequency, while the solid red line reveals a trendline calculated over ±6 frequency elements (±0.478 mHz).It can be seen that the PSD for the M0V star is relatively flat, with small-amplitude power enhancements in the range 3 − 10 mHz, which is consistent with typical p-mode oscillations.On the contrary, the PSD for the M4V star exhibits a clear enhancement of spectral energy at lower frequencies, resulting in a spectral slope of β = −0.57± 0.05 that begins at 0.32 ± 0.04 mHz, followed by numerous power peaks in the range of 1 − 10 mHz, which is consistent with the presence of both nanoflare activity and p-mode oscillations.compared to the corresponding simulated signatures, to determine which α and τ values match.The values of these nanoflare paramaters which equal all four of the observed diagnostics are the derived α and τ values.For pre-M2.5V stars, it was not possible to establish values for the powerlaw index and e-folding time that were self-consistent with the Monte Carlo models provided by Paper II.For example, it was possible to find self-similarity between the observational and model power-law indices, but this resulted in decay timescales that were incompatible and inconsistent.As a result, we are unable to define nanoflare characteristics for pre-M2.5Vstars, suggesting that nanoflare activity may be very weak (or not present) on these specific stellar sub-types.
The statistical parameters for the M2.5V, M3V and M4V stars, which are believed to be beyond the convective boundary and therefore best described as 'fully convective', exhibit values consistent with the power-law indices of α = 2.25 ± 0.25 or α = 3.00 ± 0.25, α = 2.25 ± 0.20 or α = 3.00 ± 0.20, and α = 2.30 ± 0.20 or α = 3.10 ± 0.20, alongside the e-folding timescales of τ = 200 ± 100 s, τ = 200 ± 100 s, and τ = 450 ± 50 s, respectively (see Table 2).As highlighted by Paper II, the approximate symmetry of the statistical distributions about their peak values leads to ambiguity in the derived power-law indices (see, e.g., the bands of similar values shown in each panel in Figure 5 of Paper II ).As a result, it is possible to map each subtype onto two distinct solutions for the power-law index.Irrespective of this ambiguity, both sets of possible nanoflare conditions are highly active (i.e., α > 2), in stark contrast to the effectively zero statistical nanoflare signals observed in the pre-M2.5Vspectral sub-types.The larger uncertainty in the M2.5V power-law indices are due to the larger uncertainty associated with the kurtosis value for these spectral sub-types.M2.5V stars are at the boundary of predicted full convection, so a larger spread in their nanoflare properties would be expected if full convection is the cause of the spectral 'break' in associated power-law indices.Interestingly, the M4V stars display evidence for longer efolding timescales when compared to their M2.5V and M3V counterparts.This may imply that the power-law index is marginally greater than for the earlier spectral classes.As previously discussed, constant ζ values are seen throughout the spectral sub-type range, and are thought to be due to the statistical effects of larger power-law indices being negated by the slower decay timescales associated with those stars (see Paper I for a more thorough discussion of this interplay).The specific values for the e-folding timescales for the M2.5V and M3V stars of τ = 200 ± 100 s, are consistent with the previous work of Paper II, who studied similar stellar types.
Overall, the changes in the statistical parameters indicate that post-dynamo mode transition M-dwarf stars (i.e., M2.5V and later and fully convective) exhibit greatly enhanced stellar nanoflare activity when compared to the partially convective pre-dynamo mode transition M-dwarfs that show littleto-no evidence for nanoflare activity.
As highlighted in Paper II, the examination of Fourier signatures, which are derived directly from the stellar lightcurves, can help disambiguate any derived nanoflare characteristics and further substantiate the evidence for specific activity levels.Following the methods documented by Welch (1961) and Vaughan (2012), power spectral densities (PSDs) were derived from the stellar time series.The longest continuous time series (i.e., the longest uninterrupted series of frames) common to all stars was 2095 datapoints, slightly shorter than the 2316 consecutive frames employed by Paper II.This resulted in the frequency resolution being slightly reduced from ∆f = 0.0356 mHz to ∆f = 0.0398 mHz in the present study.In order to readily compare the observational PSDs to those calculated from the Monte Carlo nanoflare models of Paper II, the Fourier signatures needed to be recalculated adhering to the new frequency resolution.Hence, utilizing the new frequency resolution, we re-computed the PSDs and corresponding 'heat map' of the simulated Fourier properties (c.f., Figure 7 of Dillon et al. 2020) as a function of both the nanoflare power-law and e-folding time.The recalculated heat map is displayed in Figure 6.Due to the change in frequency resolution being a relatively small value (0.0042 mHz), no noticeable deviations from Figure 6 and the original distribution (Figure 7 of Dillon et al. 2020) can be As with the statistical signatures shown in Figure 3, there are dramatic differences in the Fourier properties between M0V and M4V stars.As seen in Figure 4, the M0V has an effectively flat power spectrum (suggesting no nanoflare signal is present; Dillon et al. 2020), which is contrasted by the M4V star that demonstrates a spectral slope of β = −0.57± 0.05 between the frequencies ∼ 0.3 − 6.0 mHz.In Figure 4, the black crosses represent the individual frequency-dependent power measurements, while the solid red line depicts a trendline established over ±6 frequency elements (±0.478 mHz).In the lower panel of Figure 4, a PSD slope is consistent with enhanced rates of stellar nanoflare activity, which begins at the 'turning point' of 0.32 ± 0.04 mHz.As defined by Paper II, the turning point is defined as the initial peak before the gradual reduction in Fourier power with increasing frequency.It must be noted that both PSD plots shown in Figure 4 (i.e., for M0V and M4V spectral types) exhibit numerous power peaks in the range of 1 − 10 mHz, consistent with both stellar nanoflare signatures (Dillon et al. 2020) and the presence of p-mode oscillations generated in the convective layers of M-dwarf stellar sources (M-dwarf stars are believed to exhibit solar-like oscillations, hence p-modes synonymous with the typical solar frequency range; Rodríguez-López et al. 2014;Rodríguez et al. 2016).As the entire range of spectral types included in this study (M0V -M4V) are expected to exhibit p-mode oscillations, the peak frequencies within this interval are not conclusive evidence alone of nanoflare activity.
The averaged (following bootstrap procedures) Fourier properties per spectral type are shown in Figure 5, and tabulated in Table 3.As with the averaged statistical signatures shown in Figure 3, there is a marked change in Fourier features consistent with nanoflare activity for spectral classifications M2.5V and later.Evidence for this is shown in the averaged PSD spectral gradient (lower panel of Figure 5), where pre-M2.5Vstars have relatively flat spectral slopes (β ∼ 0), yet stellar sources past the convective boundary at M2.5V and later demonstrate increased magnitude spectral slopes in the range of −0.6 ≤ β ≤ −0.3.Note that the peak frequency values (upper panel of Figure 5) are relatively consistent across all M-dwarf stellar sources, approximately in the range of 2 − 4 mHz.As discussed above, this alone does not constitute evidence of nanoflare activity since all of these sources are expected to demonstrate p-mode oscillations spanning that particular frequency interval (Guenther et al. 2008;Rodríguez-López et al. 2014).
The corresponding 'turning point', where the spectral slopes are observed to begin, is, of course, equal to zero for the pre-M2.5Vstars since they do not exhibit any associated spectral slopes (middle panel of Figure 5).However, for spectral classifications beyond M2.5V,where the stars are believed to be fully convective, a relatively constant value (when errors are included) in the range of 0.3 ≤ f ≤ 0.9 mHz is found, which is consistent with the previous work of Paper II.In simulated nanoflare lightcurves documented by Paper II, an increased flare decay rate (i.e., a longer τ value) gave rise to a decreased frequency of the Fourier turning point.Examination of the middle panel of Figure 5 shows that while the turning point frequencies are distinctly different from the pre-M2.5Vstars, there does seem to be tentative evidence that the average turning point frequency decreases across the M2.5V, M3V, and M4V spectral types.This is further evidenced in Table 3, where the turning points of the M2.5V, M3V, and M4V stars are computed as 0.762 ± 0.105 mHz, 0.684 ± 0.063 mHz, and 0.467 ± 0.103 mHz, respectively.The evidence suggests that the e-folding timescales associated with the M4V stars are longer than their M2.5Vcounterparts, which is consistent with the intensity fluctuation statistical signatures discussed above.
Comparing the derived Fourier properties to the heat maps shown in Figure 6, it is possible to estimate the power-law indices and decay timescales for each of the M2.5V, M3V, and M4V stellar types that shown clear evidence for nanoflare activity.We find power-law indices of α = 3.00 ± 0.15, α = 3.00 ± 0.15, and α = 3.10 ± 0.15, alongside nanoflare e-folding timescales of τ = 200 ± 100 s, τ = 250 ± 100 s, and τ = 450 ± 50 s, for the M2.5V, M3V, and M4V spectral types, respectively (see Table 4).Importantly, these values are consistent with the statistical analyses, with the Fourier techniques providing additional benchmarks to validate the nanoflare properties extracted from the observational time series and resolve the ambiguity in power-law index arising from the statistical analysis.In contrast to the statistical mapping, the derived Fourier parameters of the M3V stars are consistent with a marginal e-folding time enhancement compared to the M2.5V classifications.This is likely related to the same physical processes that caused enhanced e-folding timescales in the M4V star.However, this is difficult to ascertain due to the relatively large errors in determining the plasma decay rate over the entire stellar surface.
Combining the Fourier and statistical analyses (see Table 5), we find that the fully convective M2.5V and M3V subtypes exhibit nanoflare power-law indices of α = 3.00±0.20 and α = 3.00 ± 0.18, respectively.The M2.5V sub-types are consistent with a decay timescale of τ = 200 ± 100 s, whereas the M3V stars display tentative evidence for a slightly enhanced e-folding timescale of τ = 225 ± 100 s.These e-folding timescales and power-law indices are values consistent with similar M-dwarf spectral types studied by Paper II, whereas M4V stars exhibit elevated power-law indices of α = 3.10 ± 0.18, with an increased decay timescale of τ = 450 ± 50 s.With these properties confirmed, the behavior of these flares in comparison to M-dwarf flare samples as-a-whole can be inferred.It has been established that a general relationship between the flare duration, t, and energy holds for observable flare populations, namely that t ∝ E x , where x ≈ 0.33 for solar and G-type stellar flares (Veronig et al. 2002;Maehara et al. 2015), whereas it drops to x ≈ 0.2 for solar microflares (Christe et al. 2008).In Chang et al. (2015), a directly comparable relationship between the e-folding time (in minutes) and the flare energy was defined from a sample of 420 energetic M-dwarf flares (E ≃ 10 31 − 10 34 erg).The log-log fit of the data was established at a high statistical significance as, log τ = (0.57± 0.05) log E − (15.61 ± 1.57) .(3) Taking the peak energies for nanoflares as E = 10 25 erg, Equation 3 produces e-folding times of 1 ≤ τ ≤ 258 s.From inspection, the derived e-folding times of the M2.5V and M3V populations are consistent with the upper boundary of predicted values, whereas the M4V values lie outwith the derived relationship.This is not necessarily unexpected, since a complementary study by Howard et al. (2019) found a broken power-law index relationship between the e-folding time and flare energy, where at E ≤ 10 33 erg, τ remained approximately constant instead of following the trend associated with Equation 3. The authors attributed this to the limitations of flare characterizations around the detection limit, but the flares in this study suggest that the effect may be physical.Equation 3 was derived from flare energies orders of magnitude above those under consideration in this study, thus the disparity between the predicted e-folding time and that seen in M4V stars may be indicative of a transition from largescale Petschek reconnection to the Sweet-Parker process.As was proposed by Tsuneta & Katsukawa (2004), small-scale pico/nano-scale flares occur more favorably via Sweet-Parker than Petschek reconnection.As Paper II suggest, this would explain a discontinuity in the power-law relationship between nanoflares and their larger scale counterparts, which remain driven by Petschek-like reconnection (Loureiro & Uzdensky 2016).The Sweet-Parker reconnection process is inversely proportional to the square root of the plasma Lundquist number, which is itself inversely proportional to the plasma resistivity.As such, Sweet-Parker reconnection is more favorable in poorly conducting plasmas.The increased decay timescale of τ = 450 ± 50 s, alongside the associated increased powerlaw index of α = 3.10 ± 0.18, found for the M4V sub-type may be related to increased plasma resistivity, which matches expectations for mid-to-late M-dwarfs (Mohanty et al. 2002).Caution is required however ; these increased α values are within 1 σ N , and the τ values within 3 σ N of the uncertainties of the less enhanced M2.5 and M3V stars, so this trend cannot yet be considered statistically significant.Future investigation of M5V and later sub-types is required to determine if there is a statistically significant trend exceeding 3 σ N confidence in the observed properties.This could be complimented by multi-color observations that would allow for lower uncertainty in the τ value at each color band due to the reliance on underlying plasma properties, which are naturally more separated across color bands due to their associated temperature sensitivities.
In contrast to the fully convective sub-types, pre-dynamo mode transition M0V -M2V stars exhibited weak (if any) nanoflare signals, suggesting that fully convective stellar atmospheres lead to a large enhancement of nanoflare activity.
While the observed trend of fully convective stars exhibiting enhanced nanoflare activity is clear, the exact mechanism leading to this is still a matter of debate.While the Sweet-Parker hypothesis is plausible, there is also a potential issue.If enhanced nanoflare activity occurs in the corona, it would lead to enhanced heating of that plasma.Consequently, this would lower the resistivity and hence lower the rate of Sweet-Parker reconnection.This 'feedback loop' behavior may reach some natural and stable equilibrium, but it may be necessary to incorporate additional theory to ensure the stability of this mechanism.Referring to the original nanoflare mechanism theorized by Parker (1988) may provide this.In that paper, Parker suggested random convective motion in the photosphere causes 'shuffling' and subsequent deformation and braiding of the photospheric footpoints of the coronal magnetic fields and consequently the generation of free energy.The coupling of the magnetic field lines between the photosphere and corona provides the framework to allow this free energy to flow into the corona.This energy is then dissipated in coronal current sheets, leading to small-scale reconnection.As such, enhanced heating leading to decreased resistivity would improve the magnetic coupling between these footpoints and the corona, consequently enhancing the flow of free energy available for nanoflare activity.One can imagine a combination of these scenarios, wherein the sympathetic transfer of hot plasma and free energy through these coupled fields regulates the resistivity and drives a stable rate of Sweet-Parker reconnection.
To uncover the source of this enhanced activity, it's vital to obtain two sets of observations; multi-band photometry, and observations of later MV star types.The multi-color observation of these stars will allow us to make a limited analysis of the change in nanoflaring properties across different wavelengths and consequently the contribution at different atmospheric heights.Comparing relative photospheric and coronal signatures could diagnose the underlying mechanism powering this enhanced nanoflare activity.The multi-color analysis should also provide a lower uncertainty in the τ values.Secondly, sourcing M5V and beyond stars would allow the continuation of the trend in flare decay rate (if any) to be investigated.If later MV stars continue to exhibit enhanced activity it would support the Sweet-Parker reconnection theory, as it would suggest the enhanced resistivity is key.Ultimately, observations of later MV stars, and across multiple photometry bands will need to be coupled with detailed physical modeling to try and uncover what changes in these stars are driving their nanoflare behavior.
Regardless of the specific physical mechanism causing this enhancement across the convective boundary, it is there.The observational evidence points to nanoflare contributions increasing significantly in the fully convective M2.5V and later stars.This novel result is independent of the modeled nanoflare lightcurves, which serve only to diagnose the parameters of the nanoflare signatures within observed lightcurves.It is also independent of the range of stellar luminosities present in the sample.It is established that greater macroscopic flare rates in later-type stars, as seen in Table A1, are influenced by the reduced luminosity threshold of these stars.It is therefore of interest to investigate whether the reduced flare detection threshold may influence the nanoflare study presented here.Average luminosities for the M0 and M4 stars, representing the largest range of spectral classes under consideration, are L M 0 = 0.068L ⊙ and L M 4 = 0.014L ⊙ .Therefore, the energy rates associated with these luminosities can be estimated as E M 0 = 2.6 × 10 32 erg s −1 and E M 4 = 5.5 × 10 31 erg s −1 , respectively.This order-of-magnitude drop in the energy rate associated with the fundamental stellar brightness agrees with previously detected flare energies (e.g., Rodríguez Martínez et al. 2020).In the present study, a 1σ N deviation is modeled in the simulation as 5 × 10 24 erg (Jess et al. 2019;Dillon et al. 2020), and by projecting an order-of-magnitude energy threshold decrease between M0 and M4 classifications to nanoflare conditions, 1σ N deviations in M4V stars equates to flares with energies ∼ 10 23 erg, still well within the typical nanoflare regime.This differential in flare energy sampling is also not sufficient to explain the lack of nanoflare enhancement in early-type MV stars.The flare frequency, dN/dE from Equation 2, associated with the M2.5V nanoflares calculated from the α and τ values seen in Table 5 reveals a two order-of-magnitude increase in nanoflare frequency between M0 and M4 (10 −46 − 10 −44 erg −1 cm −2 s −1 ), consistent with previous solar studies (Purkhart & Veronig 2022).Therefore, in a scenario where early-type M-dwarfs had the same flaring profile as their later-type counterparts, there would still be an ample frequency of 1σ N signatures present in their histograms and Fourier spectra.Given that there is no such behavior seen, we posit this as evidence that the nanoflare frequency spectrum detected in this work exists only in the fully-convective sample.Finding the source of this convective divide should be a key focus of future studies.
The enhanced small-scale flare rates in fully convective stars holds profound implications for the energy budgets of those stellar sources.The energy output of rapid and continuous nanoflares may be a major component of the overall stellar energy budget, yet are hidden within the noise envelope of the observations and can only be extracted through use of large-scale statistical and Fourier analyses.The question of whether the enhanced flaring visible in post-dynamo mode transition M2.5V -M4V stars is due to the helical dynamo or altered plasma Lundquist conditions in these stars  (2018), that solar and stellar dynamos operate independent of a tachocline.As a result, it is of paramount importance to source sufficient late M-type stellar time series for follow-up analyses.

CONCLUSIONS
Evidence for stellar nanoflares has been observed on a further 15 post-dynamo mode transition (M2.5V,M3V, and M4V classification) stars, with nanoflare power-law indices and e-folding times consistent with the enhanced rates of nanoflare activity put forward by Paper II.The marked increase in nanoflare activity is coincident with M2.5V and later sub-types, suggesting that the change from partial to fully convective atmospheres may be responsible.The postdynamo mode transition stars exhibit nanoflare rates that are enhanced from those seen at larger energies in other stars and the Sun, with power-law indices found to be in the region of α = 3.00 ± 0.20 for M2.5V and M3V sub-types, with slightly larger values of α = 3.10 ± 0.18 for M4V subtypes.Given the relation between power-law index and lowenergy flare frequency, it is clear that the atmospheres of late MV stars have an energy budget dominated by small-scale flaring.Whereas observational evidence of nanoflares being the dominant heating mechanism in the solar atmosphere remains elusive, their energy output in fully-convective stars may well be sufficient to produce bulk heating.The decay timescales for M2.5V and M3V stars were found to be on the order of τ = 200 ± 100 s, while evidence was presented for increased plasma e-folding times of τ = 450 ± 50 s in the M4V stars, suggesting the presence of Sweet-Parker reconnection processes.It must be noted these enhanced values for the M4V star remain within 1 − 3σ N of the M2.5 and M3V stars, so we cannot yet consider these to be fully distinct.
On the contrary, pre-dynamo mode transition M-dwarf (M0V, M1V, and M2V classification) stars exhibit marginal statistical or Fourier-based nanoflare signals, indicating that the large power-law index for MV stars reported in Paper II is not uniform across all spectral types.Instead, it is implied that a fully-convective interior is necessary to exhibit the α ≥ 3 that distinguish them from other stellar candidates.Additionally, the underlying reason why fully convective atmospheres lead to enhanced nanoflare activity should be explored, i.e., is this due to an altered dynamo, or due to other plasma changes such as modification of the corresponding Lundquist number?One avenue of exploration would be examining M5V (and later) stellar types, to investigate if there is a continuing and more statistically significant trend in the flare decay rate and associated power-law index, which could be linked to increasing plasma resistivity, and thus increased Sweet-Parker reconnection rates.It is likely such observations would need to be coupled to detailed theoretical and modeling efforts using well-developed numerical simulations (e.g., Takahashi et al. 2011;Tenerani et al. 2015;Shi et al. 2018;Papini et al. 2019).
Additionally, sampling late-type M dwarfs may reveal the traditional observational signatures of nanoflares.Since the advent of M-dwarf flare studies following the seminal works of Gershberg (1972) and Lacy et al. (1976), it has been clear that the higher luminosity of early-type MV stars can skew flare population studies in favor of less luminous later-type stars due to a lower intensity threshold for flare detection.This has led to subsequent studies constraining the population of stars under consideration, or studying a single star in detail (e.g., Hawley et al. 2014b;Davenport et al. 2014).This is the case in our present study, where we limited the population of MV stars to those with comparable magnitudes.However, utilizing the reduction in surface temperature in late-type MV stars would allow for the intensity threshold of a macroscopic flare detection to be reduced.A sample of MV stars with varying magnitudes would allow for the 2σ N nanoflare intensity excursions found in this study to potentially lie above the minimum observable detection threshold in cooler late-type MV stars.This could confirm the temporal morphology and occurrence frequencies of flares at these energies and provide further insight into the validity of previously defined relationships, such as Equation 3. Additionally, the statistical techniques employed here could provide the first signatures of even smaller flares, such as the proposed pico-flare energy regime (Katsukawa & Tsuneta 2001;Katsukawa 2003).
It goes without saying that enhanced small-scale reconnection in fully convective stars may mean that nanoflare activity could be a significant component of their overall energy budget.Large-scale multi-year studies of stellar nanoflare rates in fully convective M-dwarfs would further our understanding of nanoflare behavior across different activity cycles, which would further shine light on the ubiquity and role nanoflares play in these dynamic host stars.This can be achieved through further use of large-scale sky surveys (like the NGTS) and space-based observations from the likes of the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014), alongside targeted campaigns using high-cadence observational platforms, such as HiPERCAM (Dhillon et al. 2016), or multi-band photometry such as the Rapid Eye Mount (REM) telescope (Antonelli 2005)

Figure 1 .
Figure1.Sample lightcurves of NGTS J052346.3-361114(M0V-type; top panel), and NGTS J050423.8-373021(M4V-type; lower panel) spanning 24 000 s, each normalized by their respective standard deviations, σN .The region highlighted in red in the lower panel denotes the intensity values that are removed from consideration due to an excursion above 3σN , as is the formality for macroscopic flare signatures.

Figure 2 .
Figure 2. Histograms of intensity fluctuations, each normalized by their respective standard deviations, σN , for the NGTS J052346.3-361114(M0V-type; top panel), and NGTS J050423.8-373021(M4V-type; lower panel) lightcurves.A standardized Gaussian profile is overplotted in each panel using a dashed red line for reference.The M4V-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 M0Vtype intensity fluctuations provide effectively zero negative median offset, and no elevated occurrences at ∼ 2 σN .This is inconsistent with clear statistical signatures of nanoflare activity, with the resulting distribution remaining more 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 negative median offsets and occurrence excesses, respectively, are clearly visible for the M4V stellar source.For improved clarity, the blue and gold lines display the corresponding distributions in each zoomed panel.
(2011)  andJess et al. (2014) alongside Paper I and Paper II, and allows for sensitivity in the detection of any residual intensity deviations within the expected intensity ranges of a normal distribution.The number of macroscopic flares removed were used to calculate approximate flare rates for the M stars, which are displayed in Table

Figure 3 .
Figure 3.The bootstrap-averaged statistical properties of the intensity fluctuation histograms for each stellar classification.Beyond the convective boundary, at approximately M2.5V and later, sub-types begin to exhibit statistical signatures that are consistent with the presence of nanoflare activity, including larger median offsets (top panel), increasing levels of kurtosis (second panel from top), and higher Fisher skewness values (second panel from bottom).The ζ (FW 1 8 M-to-FWHM ratio) values do not vary significantly as a function of stellar classification.However, this is likely due to the interplay between the power-law index of the nanoflares and the duration of the e-folding timescales, which are able to counteract the statistical effects of one another.demonstrate a larger consistent offset magnitude (with less uncertainty) of approximately −0.05σ N .The Fisher skewness value is effectively zero for pre-M2.5Vstars (second panel from bottom in Figure3), suggesting no, or very weak, nanoflare activity.From M2.5V onward, there is a clear increasing trend in the Fisher skewness value of the fluctuation distribution, with the M4V subtype displaying a Fisher skewness equal to 0.051 ± 0.014, providing strong evidence for the presence of nanoflares.In the additional distribution diagnostics, the relationship is less clear.Regarding the kurtosis (second panel from top in Figure3) there appears to be a trend in that statistical kurtosis increases across the full sample, between M0V and M4V.However, the exact nature of this relationship is obscured by the large uncertainty in the kurtosis for spectral types around the dynamo mode transition.In particular, the M3V sub-type appears to exhibit a decrease in kurtosis, though this may be due to the large uncertainty in M2V & M2.5V producing abnormally large values.There is no clear trend visible in the corresponding ζ values (lower panel of Figure3) as a function of spectral subtype.It must be remembered that the ζ value is a measure of the deviation away from a standard Gaussian distribution, which has a value of ζ = 1.73.As discussed in Paper I, in- ) as a function of spectral subtype.It must be remembered that the ζ value is a measure of the deviation away from a standard Gaussian distribution, which has a value of ζ = 1.73.As discussed in Paper I, in-creased nanoflare decay timescales (i.e., larger τ values) result in broader tails of the intensity fluctuation distributions, hence giving rise to ζ > 1.73.On the contrary, large powerlaw indices help reduce the widths of the tails in the intensity fluctuation distributions due to the superposition of positive intensity fluctuations (e.g., new nanoflares) superimposed on top of decaying (i.e., negative) intensity fluctuations, which result in ζ < 1.73.As such, the interplay between the powerlaw index and the nanoflare e-folding time produces the specific value of ζ measured.As such, the relatively consistent values of ζ found across the spectral range M0V -M4V may result from an increased nanoflare rate expected for M4V stars being negated by an increase in the associated decay timescales of the resulting nanoflares, i.e., a larger α term being coupled with longer τ values.

Figure 4 .
Figure 4.The Fourier power spectral densities (PSDs) for example M0V (upper panel) and M4V (lower panel) stellar sources, displayed in normalized units of σ 2N /mHz.The crosses in each panel depict the individual power values as a function of frequency, while the solid red line reveals a trendline calculated over ±6 frequency elements (±0.478 mHz).It can be seen that the PSD for the M0V star is relatively flat, with small-amplitude power enhancements in the range 3 − 10 mHz, which is consistent with typical p-mode oscillations.On the contrary, the PSD for the M4V star exhibits a clear enhancement of spectral energy at lower frequencies, resulting in a spectral slope of β = −0.57± 0.05 that begins at 0.32 ± 0.04 mHz, followed by numerous power peaks in the range of 1 − 10 mHz, which is consistent with the presence of both nanoflare activity and p-mode oscillations.

Figure 5 .
Figure5.The bootstrap-averaged properties of the Fourier power spectral densities (PSDs) across each spectral type.The upper panel displays the peak frequency values (in mHz), which are found to reside within the range of approximately 1 − 4 mHz, which is consistent with both nanoflare activity and p-mode oscillations, and therefore cannot be used as an indicator of nanoflare activity by itself.The middle and lower panels display the turning point frequencies (in mHz) and subsequent spectral slopes, respectively, as a function of stellar classification.When compared to the Monte Carlo nanoflare simulation outputs depicted in Figure6, the distinct jump in turning point frequency and spectral gradient at the convective boundary (M2.5V)provides clear evidence of prominent nanoflare activity in M2.5V -M4V stellar sources.
Figure 6 shows a 'heat map' of the simulated PSDs (c.f., Figure 7 of Dillon et al. 2020), which has been recalculated for the 2095 element long time series employed in the present study.While the frequency resolution is slightly reduced (∆f = 0.0398 mHz) than that utilized by Dillon et al. (2020, ∆f = 0.0356 mHz), the overall trends and evolution remain consistent across the power-law index and e-folding timescale values.

Figure 6 .
Figure 6.A reproduction of Figure 7 from Paper II, with the constituent PSDs re-calculated for 2095 datapoints to match the longest continuous time series used in the present study.The primary peak frequencies (lower-left), spectral slopes (upper-left), dominant frequencies following detrending (upper-right), and the percentage of nanoflare power above the noise floor in the range of 1 − 5 mHz (lower-right), is displayed as a function of the power-law index, α, and the decay timescale, τ , used to generate the synthetic time series.While a few individual values differ, the overall trends and the magnitude of the derived signals are consistent with the PSD properties generated from 2316 datapoints and reported by Paper II.

Table 1 .
Averaged characteristics of the statistical properties by each spectral type.

Table 2 .
Nanoflare parameters per spectral type, derived from statistical properties of Monte-Carlo modeled nanoflare timeseries.The approximately symmetrical distribution of statistical properties leads to an ambiguity in the derived powerlaw indices, hence α1 and α2.

Table 3 .
Average characteristics of the Fourier PSD properties by each spectral type.

Table 4 .
Nanoflare parameters per spectral type, derived from Fourier properties of Monte-Carlo modeled nanoflare timeseries.There is no ambiguity in the derived powerlaw indices.

Table 5 .
Wright & Drake (2016)per spectral type, derived from combined statistical and Fourier properties of Monte-Carlo modeled nanoflare timeseries.Mohanty et al. 2002).If the nanoflare e-folding times continue to increase with increasing M-dwarf sub-type, it would support the scenario of increased plasma resistivity leading to increased small-scale flaring via Sweet-Parker reconnection.This would appear to support the findings ofWright & Drake (2016);Wright et al.
to investigate the nanoflare signature across layers of the Stellar atmosphere.