A study of flares in the ultra-cool regime from SPECULOOS-South Free

ABSTRACT

1 INTRODUCTION

In the search for planets capable of supporting life, ultra-cool dwarfs (UCDs) make compelling hosts. In our local stellar neighbourhood, UCDs are plentiful (Chabrier 2003; Henry 2004) and predicted to host large numbers of planets (Cantrell, Henry & White 2013; Dressing & Charbonneau 2015; Ballard & Johnson 2016). Moreover, their small sizes and low temperatures make it easier to detect Earth-sized, habitable-zone planets around these objects, and to probe those planets’ atmospheres for biosignatures (Kaltenegger & Traub 2009; Seager, Deming & Valenti 2009; de Wit & Seager 2013), than for any other type of star.

Despite their promise, serious questions remain about the habitability of planets around red dwarfs. Several authors find it unlikely that these extremely cool, red stars would provide their planets with enough UV photons to initiate specific pre-biotic chemistry pathways (Rimmer et al. 2018) or enough visible light for photosynthesis to occur (Lehmer et al. 2018; Mullan & Bais 2018; Covone et al. 2021). UCDs are also especially active objects (West et al. 2015; Williams et al. 2015; Gizis et al. 2017; Paudel et al. 2018; Günther et al. 2020), producing energetic stellar flares and large-scale photometric variability, which results in treacherously variable conditions for their planets. The dynamic relationship between a planet and its host is a major factor in evaluating how conducive a planet’s environment is for life. For UCDs, however, this relationship is poorly understood.

Stellar flares may be a major determinant in whether a planet can initiate and sustain life. Flares are explosive events caused by magnetic recombination in the upper atmosphere of a star (Benz & Güdel 2010). This sudden eruption of magnetic energy usually occurs around active regions, such as stellar spots, and causes bursts of particles and electromagnetic radiation. The spectra of this radiation resembles a blackbody with an effective temperature of 9000 K (Shibayama et al. 2013). Often, though not always, powerful flares will come accompanied with a coronal mass ejection (CME) event, where clouds of charged particles are directionally ejected from the star. While flares can be destructive through atmospheric erosion (Lammer et al. 2007), ozone depletion (Segura et al. 2010; Tilley et al. 2019), and even extinction events, they can also be an essential power source for life. It is possible that flares could provide the missing energy at the bluer end of the spectrum needed for cool, red dwarfs to initiate pre-biotic chemistry (Buccino, Lemarchand & Mauas 2007; Ranjan, Wordsworth & Sasselov 2017; Rimmer et al. 2018) and for photosynthesis (Mullan & Bais 2018; Lingam & Loeb 2019a). The additional UV energy may also affect the evolution of a planet’s atmospheric chemistry (Segura et al. 2010; Vida et al. 2017). In summary, stellar flares have extreme and far-reaching consequences on the planets hosted by cool stars. Therefore, it is essential to study these stars’ flaring activity, and assess how applicable our current understanding of stellar activity is to the ultra-cool regime.

Stellar activity and rotation are also closely connected. Magnetized stellar winds dictate the loss of angular momentum, which slows a star’s rotation over time. These stellar winds are dependent on the structure and properties of the magnetic field. As the rotation slows, this decreases the magnetic activity in a process known as spin-down (Skumanich 1972; Noyes et al. 1984). Due to this effect, rapidly rotating stars flare much more frequently than slow rotators (Skumanich 1986; Davenport et al. 2019; Mondrik et al. 2018; Medina et al. 2020). Spots and faculae on the surface of stars and brown dwarfs cause periodic photometric variations as they come in and out of view, on the same time-scale as the object’s rotation. This allows the rotation period to be deduced directly from the photometry. Therefore, the rotation of a star can provide valuable insights into the magnetic dynamo, the mechanism that generates a star’s magnetic field.

This magnetic dynamo is poorly constrained for fully convective low-mass objects [with masses ≤0.35 M_⊙; Chabrier & Baraffe (1997)], where it is believed to differ significantly from the solar model. Despite this predicted difference, recent work has shown that relationships between activity and rotation remain consistent from partially to fully convective stars (Wright & Drake 2016; Newton et al. 2017; Wright et al. 2018). Spin down is, however, believed to occur on slower time-scales for fully convective stars, which accounts for the enhanced activity of mid-to-late M dwarfs (West et al. 2008; Newton et al. 2017; Jackman et al. 2021). Therefore, rotation, activity, and the relationship between the two, are extremely useful probes of the underlying magnetic activity of UCDs.

Due to the promising nature of M dwarfs as planetary hosts, there have been several detailed studies of their flaring activity within the past decade. Space telescopes, such as the Kepler/K2 missions (Borucki et al. 2010; Howell et al. 2014), allowed the first insights into the flares of bright M dwarfs (Davenport et al. 2014; Hawley et al. 2014; Lurie et al. 2015; Silverberg et al. 2016) and of (small numbers of) UCDs (Gizis et al. 2013, 2017; Paudel et al. 2018). Additionally, the recently launched Transiting Exoplanet Survey Satellite [TESS; Ricker et al. (2015)] has facilitated studies of the flaring and rotating activity of cool stars (Günther et al. 2020; Medina et al. 2020; Seli et al. 2021). Several ground-based photometric surveys, such as MEarth (Nutzman & Charbonneau 2008), the Next Generation Transit Survey (Wheatley et al. 2018), the All-Sky Automated Survey for Supernovae [ASAS-SN; Shappee et al. (2014)], and Evryscope (Law et al. 2015), have also carried out detailed flare studies that include M dwarfs (West et al. 2015; Howard et al. 2019, 2020a, b; Mondrik et al. 2018; Martínez et al. 2019; Schmidt et al. 2019; Jackman et al. 2021).

The most frequent stellar flares are small, fast, and difficult to detect above photometric scatter (Lacy, Moffett & Evans 1976). Due to the intrinsic faintness of UCDs in the visible, it is difficult to achieve the high photometric precisions necessary to constrain their flaring activity. However, it is possible to perform small, dedicated studies of the much rarer, high-energy flares on UCDs (Gizis et al. 2013; Paudel et al. 2018; Jackman et al. 2019). Therefore, to obtain a sufficient sample size, previous large flare studies have focused on hotter stars, up to mid-M dwarfs, where photometric precisions are much higher. This has resulted in limited flare statistics for the coolest stars, for which we would require a large, high-cadence photometric survey optimized for UCDs.

The Search for habitable Planets EClipsing ULtra-cOOl Stars (SPECULOOS) project (Gillon 2018; Burdanov et al. 2018; Delrez et al. 2018; Jehin et al. 2018; Sebastian et al. 2021) aims to search for transiting planets around the nearest (within 40 pc) UCD stars. SPECULOOS’s motivation is to provide temperate, terrestrial planets for detailed atmospheric characterization with the James Webb Space Telescope (Gardner et al. 2006) and future extremely large telescopes. The SPECULOOS target catalogue of ultra-cool objects is defined in Sebastian et al. (2021). SPECULOOS comprises a network of 1-m class telescopes spread across the Northern and Southern Hemispheres. The largest facility, the SPECULOOS-Southern Observatory (SSO) in Cerro Paranal, Chile, consists of four identical, robotic telescopes. These telescopes are named Io, Europa, Callisto, and Ganymede. Each telescope operates independently and in a robotic mode, following plans written by SPECULOOS’s automatic scheduler, SPeculoos Observatory sChedule maKer (Sebastian et al. 2021). Additionally, the SPECULOOS-Northern Observatory (Niraula et al. 2020) in Tenerife and the Search And characterIsatioN of Transiting EXoplanets (SAINT-EX) facility (Demory et al. 2020) in San Pedro Mártir each have one 1-m telescope. The SSO began official scientific operations in 2019 January, followed by SAINT-EX in 2019 March and the SNO in 2019 June.

In this paper, we present a study of flares and rotation of SSO targets, observed over a span of almost 2 yr. Section 2 defines the SSO data sample. Section 3 outlines how we generate and clean our global light curves, and how we combine an automated flare detection algorithm with manual vetting to obtain our final flare sample. We describe modelling the flares, calculating the flare energies, and estimating flare rates. This section also includes an assessment of the flare sample’s completeness. In Section 4, we measure the rotation periods in our sample. The results of our flare and rotation analyses are presented in Section 5 and discussed in Section 6, which includes a contextualization of the impact of flares on the potential for life on planets around UCDs.

2 SPECULOOS-SOUTH DATA SAMPLE

For this study, we defined a data set spanning from 2018 June 1 to 2020 March 23, using observations from all four telescopes in the SSO. While this start date is before official scientific operations began, it marks a point of stability in the commissioning phase. After this date, we performed no major maintenance to the observatory; donuts (McCormac et al. 2013), our auto-guiding software for precise pointing, was in use; and our operating strategy had been finalized [see Sebastian et al. (2021) for details]. The end of the data package was defined as the date that ESO Paranal Observatory shutdown due to the coronavirus pandemic.

In the time period between 2018 June 1 and 2020 March 23, there are 661 potential nights of observation. Ganymede, the final telescope installed in the SSO, started commissioning on 2018 September 30, therefore, between the four telescopes there is a cumulative total of 2523 nights. However, primarily due to weather loss, the SSO has lost 23–24 per cent of observing time, resulting in 1931 nights of observation. Over the 1931 (combined) nights in this sample, we have observed 176 unique photometric targets for at least one night. These observations have typical exposure times of 20–60 s. The majority, though not all, of this data sample are in the SPECULOOS target list (86 per cent), as defined in Sebastian et al. (2021).

This target list is divided into three main scientific programs, which are described in more detail in Sebastian et al. (2021). The objects in the SSO sample that are not in the target list are exclusively objects with spectral types around or earlier than M6. It is likely that these objects were part of the commissioning phase of the telescopes or were removed from the target list (due to reclassifying their spectral type) before official scientific operations began. For the objects in the target list, we extract radii, masses, effective temperatures, and spectral type classifications from Sebastian et al. (2021). The stellar parameters for non-target-list objects are calculated in the same way.

We remove any objects that are not part of SPECULOOS’s usual ‘survey mode’, such as targets observed for follow-up and monitoring of TRAnsiting Planets and PlanetesImals Small Telescope (TRAPPIST)-1’s transits. We note that, in SPECULOOS’s ‘survey mode’, there are no simultaneous observations with multiple telescopes during the course of a night. During this time period, the operational temperature of the CCD was increased from −70 to −60°C in 2018 October. The choice to raise the temperature of the CCD was due to the effect on the quantum efficiency of the detector, improving our sensitivity at the red limit, while the increase in dark current was found to be negligible. This temperature change introduces an offset in the differential flux between the nights before and after. While this offset has no consequence on our flare program, it affects long-term photometric trends. Therefore, when recovering rotation periods (in Section 4), we split the light curves that straddle this temperature change and analyse the before and after sections independently. As this temperature change happens during the first few months of our data set, the majority of targets were either observed entirely before or after it. Therefore, we only need to split a handful of our light curves. If a target has been observed by multiple telescopes, we do not combine their light curves, as each will experience slightly different instrumental systematics. However, we do check that the rotation periods estimated in Section 4 are present in observations from all telescopes.

If we choose only the objects that have been observed more than 20 h with one telescope, this provides us instead with 154 targets (of which 134, or ∼87 per cent, are in the SPECULOOS target list). We define the observation time as the sum of the span of all observation nights (start of night to end of night), where we exclude any gap longer than 15 min. These 154 objects define the SPECULOOS-South data sample. This sample covers a range of M dwarfs in spectral type, extending from M4 into the early L dwarf regime (with masses of 0.07–0.2 M_⊙), as shown in Fig. 1. All objects in this sample will therefore be fully convective, as the convection limit occurs around 0.35 M_⊙ (Chabrier & Baraffe 1997). While M4 and M5 objects are not considered UCDs, we include them in this sample to explore any differences between mid-M, late-M, and L dwarfs.

Figure 1.

Spectral-type distribution of the objects in the SPECULOOS-South data sample from 2018 June 1 to 2020 March 23. We indicate which objects are in the SPECULOOS target list (Sebastian et al. 2021).

Open in new tabDownload slide

In this sample, 80 per cent of the objects have been observed for less than 77 h (Fig. 2), and 50 per cent have been observed for less than 51 h. The short observation times in this sample will inevitably bias us towards detecting flares on targets with high flare rates in Section 3.3 and detecting the clearest rotation periods for fast rotators in Section 4.

Figure 2.

A cumulative histogram of the number of days observed for the objects in the SSO data sample. Of the 176 objects, 14 per cent (24) have been observed for less than 20 h, 50 per cent have been observed for less than 51 h, and 80 per cent have been observed for less than 77 h.

Open in new tabDownload slide

From Sebastian et al. (2021), we confirm that none of the target list objects in the SSO data sample are known binaries. We could not find any record in the literature of binarity for any of the 20 non-target list objects, except from the M4 star, Gaia DR2 4450376396936878336 (Lépine & Bongiorno 2007). We flag this object, however, we find no evidence that this object behaves anomalously. This star is flaring, however, it only flares twice, therefore, it is not included in our analysis in Section 5.3. We also do not find a rotation period for this object, therefore it does not impact our results in Section 4.

3 GENERATING THE SPECULOOS-SOUTH FLARE SAMPLE

3.1 Global light curves

To search for flares, we extract global light curves for all targets in the data sample, using the SSO Pipeline described in Murray et al. (2020). The SSO Pipeline performs automated differential photometry by generating an ‘artificial light curve’ from a normalized, weighted average of reference stars in the field. Only objects above a brightness threshold (optimized for each field) are chosen as reference stars. Reference stars whose differential light curves exhibit photometric variability, or that are far in spatial distance on-sky from the target, are then weighted down in an iterative algorithm, based on Broeg, Fernández & Neuhåuser (2005). It is worth noting that the target star is removed from this process. This means that, except for the distance weighting, this ‘artificial light curve’ is entirely independent of the target star’s light curve. This pipeline then corrects for the second-order effects of atmospheric water absorption on our differential photometry, using measurements of precipitable water vapour from the ground-based Low Humidity and Temperature PROfiling radiometer instrument (Kerber et al. 2012).

To generate global light curves, we perform our differential photometry process on the entire time series of a target at once. We obtain differential light curves for 13 different aperture sizes. For each aperture, the aperture size, comparison stars, and weightings do not change over the time span of observations. Using global light curves allows us to study their long-term photometric variability, but restricts our ability to optimize the photometry night by night.

To mitigate any systematics caused by changing atmospheric conditions, we carefully select the ‘best’ aperture for each target. As we only use one aperture size across all nights, it is important to take time at this step. Choosing the best aperture for a target involves a fine balance between an aperture that is large enough to avoid losing flux for nights with sub-optimal observing conditions (e.g. a larger seeing), and an aperture that is small enough to prevent blending from nearby stars and extra white-noise contamination from the background. We select the aperture by eye that minimizes the correlation between the light-curve flux and the seeing for the whole light curve, and has no clear contamination from neighbouring stars. In addition, we implement a bad weather flag and thresholds for the background sky level, airmass, and seeing (above which the data are severely impacted), described in more detail below.

3.2 Cleaning the light curves

3.2.1 Removing ‘bad’ observations

In the context of this paper, ‘bad’ observations are defined as those that are significantly affected by the observing conditions of the night, where distinguishing real stellar variability from ground-based systematics would be extremely difficult.

To reduce the impact of these observations on our flare study, we implement a bad weather flag, as defined in Murray et al. (2020). As previously mentioned, the artificial light curve has no relation to the target’s light curve, therefore it is not impacted by flares on the target star. This distinction means that we can use the artificial light curve as an independent reflection of the photometric conditions at that time to identify where the quality of our light curves deteriorates. Data points in our light curves are flagged based on the RMS of the section of the artificial light curve surrounding that data point (local RMS). We define the length of the local section to be ±0.01 d. We make the assumption here that there is a threshold for photometric scatter in the artificial light curve (which therefore must be present in the highly weighted reference stars’ light curves) above which we cannot obtain good photometric precision. We define this threshold as a local RMS of 8 per cent, as defined in Murray et al. (2020).

We also include strict cuts on specific observation parameters. We exclude data points in our light curves where the sky background level is greater than 4000 counts per pixel, the airmass is greater than 2.5, or the full-width at half maximum (FWHM) of the point spread function is greater than 2.6 arcsec (the SSO cameras have a pixel scale of 0.35 arcsec pixel⁻¹).

From 2019 May 1 to June 19, we experienced an issue with moving dust on the CCD window of Ganymede. While stationary dust can be easily corrected by flat images taken at the start of each night, if dust moves across the CCD window during the night, this leads to residuals in the flat correction and structures in the final differential light curves. These structures can mimic planetary transits, though the frames affected by dust are easy to identify from the raw images. However, as we do not currently have a robust correction technique for moving dust, we removed all observations taken by Ganymede during this period. This amounts to less than 2 per cent of total observations, therefore the impact is marginal.

3.3 Flare detection

As the SSO data sample is constrained to only 154 targets, it is possible to identify flares manually. Therefore, we did not find it necessary to implement a fully automated, complex flare-detection algorithm. Instead, we decided to perform flare detection in two parts: a simple, automated flare-detection algorithm to extract all the flare candidates, followed by a manual vetting process to confirm them. Both stages of this flare detection process are demonstrated in Fig. 3 for target Gaia DR2 3200303384927512960. We note that this target is one of the most frequently flaring objects studied by Seli et al. (2021). This object is a particularly challenging case due to its rapid variability, where flare decays are difficult to separate from photometric modulations.

Figure 3.

Demonstration of the flare detection process on six nights of observations of an M7 object (J = 10.7 mag), which exhibits both frequent flaring and short-period rotation. The four rows show different views of the same global light curve, zoomed in to individual nights on the second row, and only the five manually verified flares on the bottom two rows. The grey points are the unbinned data, while the black points are binned every 5 min. We do a simple least-squares fit of a sine wave with a period of 0.33 d, shown in light blue. It is clear from this fit that the rotation pattern is not perfectly sinusoidal; however, it provides a way to visualize the periodicity. The medium-blue data points are those initial flare regions flagged by the automatic flare detection algorithm described in Section 3.3.1. The vertical, medium-blue lines are the flare candidates that have also been confirmed manually. During visual inspection of this light curve, we remove two flares (one in the first plot on the second row and the first flare in the second plot on the second row), due to their small amplitude, which were found too difficult to detach from the photometric scatter. The best-fitting Davenport et al. (2014) models for each flare are shown in orange. The case where we have a poor fit to the template is shown in the bottom-left plot.

Open in new tabDownload slide

3.3.1 Automatically searching for flares using light-curve gradients

It is possible to have significant time gaps within the light-curve sections of 60 points due to bad weather, however, this scenario is rare. Including data points that are far apart in time would lower the reliability of the RMS in this section of the light curve, however, when weather conditions are adverse enough to close the telescope, or seriously affect photometric precision, the RMS directly before or after is likely to be enhanced (due to poor observing conditions that do not quite meet the threshold to be removed). We note that the ‘nearest 60 points’ are only used for this flare detection method, following Lienhard et al. (2020), whereas time boxes are used for calculating the running median and RMS in Sections 3.3.2 and 3.3.3 and the rest of this work.

In their paper, Lienhard et al. (2020) determine A_thresh = 12 by inspection of the smallest flares that they intended to remove; however, this value is specific for the photometric precision and typical exposure times of TRAPPIST light curves. We would expect the precisions obtained by TRAPPIST’s 60-cm telescopes to differ significantly from those achieved by SPECULOOS’s 1-m telescopes. Instead, for the SSO light curves, we derived a value of A_thresh = 5. We find a lower value for A_thresh partly due to the increased photometric precision of the SSO and partly because we are choosing to manually verify the flare candidates that this algorithm detects, so we can afford to overdetect at this stage.

By choosing a flare detection method based on the shape of a flare, and not outlier detection, this may limit the diversity of flare morphologies that we will detect. Since there have been few large-scale studies done to date on UCD flares, such as those in our SSO sample, any difference in flare structure between earlier- and later-spectral-type objects is still unknown. Therefore, we decided to focus on the flares most resembling the standard flare shape (Moffett 1974), while flagging more unusual light-curve behaviour in the manual vetting stage (see Section 3.3.4).

3.3.2 Small versus large flares

To both maximize the detection of small flares and optimize the modelling of high-amplitude flares with slow recovery times, we split the flare candidate sample into two categories: small and large flares. We classify large flare candidates as those with at least two data points in the flare region more than seven times the running local standard deviation (standard deviation of the surrounding 80 data points) above the running local median (defined similarly as the median of the surrounding 80 data points); otherwise, we classified it as a small flare candidate. The flare region starts at least two points before the peak. We approximate the end of the flare region by using least-squares to fit the flare decline with a sum of a fast and a slow-decaying exponential, as in Davenport et al. (2014). The decay time of the slower exponential decay is used to estimate the time of the end of the flare region.

3.3.3 Validating and vetting flares

Once a flare candidate is classified as large, then we use the surrounding 320 points (defined as global), instead of 80 points (local), to calculate the running median and standard deviation for the next step of validation.

If the criteria in Section 3.3.2 is met, we then run several additional automatic quality checks to separate flare candidates from flare-like signals replicated by noise. By ensuring that there is more than one data point in the flare region after the peak, we remove cosmic rays from the flare candidates. Additionally, at least two points in the flare region must be more than two times the running local standard deviation above the running local median of the light curve. The local standard deviation and median are as defined in Section 3.3.2. When calculating the running local median and standard deviation, we mask all potential flare candidates in the light curve, and then interpolate over the flare region.

This automated algorithm provides us with a collection of flare candidates to follow up with manual vetting to then obtain our final flare sample.

3.3.4 Manual vetting

Once we have a collection of flare candidates from the method detailed in Sections 3.3.1 and 3.3.3, then the second part of our flare detection process is to manually inspect these candidates and obtain our final flare sample. We consider each flare candidate and only confirm those which can be clearly identified as matching the standard flare shape of a sharp flux increase followed by a slower, exponential flux decay (Moffett 1974; Davenport et al. 2014). We also ensure that the flares cannot be attributed to rapid changes in atmospheric conditions.

As a final validation step, we check through the global light curves in the SSO sample. The only flare-like structures remaining are due to the following events: a flare occurs before the start of the night (where we catch the tail end of the flare) or at the very end of a night, the structures are low-amplitude and so cannot be distinguished from noise, the structure diverges significantly from the standard flare model, or they are correlated with systematics. Examples of structures identified in the SSO data set that deviate from the standard flare model, and are not clearly overlapping flares, are symmetric flares (with similar rise and decay times) and flares with a sudden increase and linear decline. During this manual vetting stage, we do not add any of these potential missed flares into the flare sample, as they do not meet the criteria to be clearly identified as a flare with the standard flare shape.

3.4 Modelling flares

In order to model each flare, we coarsely flatten the light-curve surrounding that flare by dividing the light-curve by a median filtre, using the local or global criteria as in Section 3.3.2.

We chose to model our flares by fitting the empirical flare template described in Davenport et al. (2014). In this paper, the authors generate a median flare template from Kepler observations of 885 ‘classical’ flares on the M4 star GJ 1243. They found a sharp flux increase, modelled with a fourth-order polynomial, followed by an initial fast exponential decay and subsequent slower exponential decay. This model requires a time for the peak of the flare, a relative amplitude on the normalized light curve, and an FWHM_flare, which corresponds to the flare’s decay time. We demonstrate fitting this flare model to flares of differing amplitudes and decay times in Fig. 3.

While this model works well for the more classically shaped flares, it struggles to represent complex flares that do not fit the standard flare morphology (Davenport et al. 2014). This template is yet to be tested on a statistically large number of low-mass objects, such as UCDs, whose flare profiles may differ significantly. Complex flares include those with multiple peaks (Davenport et al. 2014), oscillations (Anfinogentov et al. 2013), and those that are closely entangled with variability. The overlapping of multiple flares is likely the culprit of more unusual light-curve structure. During this manual vetting, we divide flare regions containing clear multiple peaks into separate flares and each of the new separated flares are added to the flare list. We also flag (but do not remove) 27 flares that either do not fit the standard flare model, or are not easily separable into multiple flares. This corresponds to 11 per cent of the total flare sample. An example of a possible flare overlap is shown in the bottom-left plot in Fig. 3.

3.5 Extracting flare energies

For this integral, we decided to use the normalized light curve (corrected with the same local or global median filter described in Section 3.4) directly for the energy calculation, rather than using the flare model. In doing so, we note that the energies we calculate provide a lower bound on the flare energy, as we may miss the flare’s peak. It is also possible that if there was a large flare on a rapidly rotating star, the underlying F_mean may ‘smooth out’ photometric structure, leading to the flare energy being overestimated. However, to be classed a ‘large’ flare, a flare needs to be more than seven times the running global standard deviation (as defined in Section 3.3.2). This standard deviation would be highly inflated from the periodic photometric modulations, and therefore this flare would have to be very high energy for this to become an issue. As high energy flares are rare (Gershberg 1972; Lacy et al. 1976), this scenario is not a major concern. We chose not to use the flare model because the flares in our light curves are well-sampled, and when there is increased photometric noise in the light curve, or more complex flare shapes from overlapping flares (see Section 3.4), the fit to the flare model can be unreliable. Integrating the light curve directly means that bursts of flares occurring in short succession are counted as one flare of increased energy. While this simplification may affect the calculation of our flare rates and energies, it should not significantly affect our results as we only detect 27 unusually shaped flares that may be a blend of multiple flares. We note that while flares are typically approximated as blackbodies with temperatures from 9000 to 10 000 K, the flare temperature has been observed to vary outside this range, both between and within flare events (Howard et al. 2020b).

3.6 Calculating flare rates

As flares are stochastic events, the gaps in observations due to the day-night cycle or bad weather loss do not hinder our ability to calculate flaring rates, which are only dependent on the total time on-sky. It is likely that we underestimate flaring rates for the several reasons flares are not found by our automatic detection algorithm, outlined in Section 3.3.4.

We calculate our flare rates as the number of flares detected per target divided by the observation time (summed across all telescopes). Almost all of our targets have been observed for less than 200 h, therefore we have a lower limit on individual stars’ flaring rates of ∼0.12 d⁻¹.

3.7 Completeness of the flare sample

To evaluate the completeness of our flare sample (the minimum energy flares recoverable by our flare detection process for different spectral types), we perform artificial flare injection-recovery tests. Using the global light curves from the 154 different targets with at least 20 h of observation, we mask all of the flares detected in Section 3.3. We use the Davenport et al. (2014) flare model to generate 100 000 artificial flares with amplitudes drawn from a lognormal distribution between relative fluxes of 0.001 and 5, and FWHM_flare drawn from a uniform distribution between 30 s (typical exposure time for SPECULOOS) and 1 h. We calculate the energy of the resulting flares for each star, and divide into six energy bins in log₁₀ space from 10²⁸ to 10³⁴ erg.

We note that here we have followed the work of Davenport et al. (2014), Davenport (2016), Günther et al. (2020), among others, and decided not to consider any relationship between amplitude and duration. While we may be injecting ‘unphysical’ flares, this approach allows us to explore a larger parameter space of different flare morphologies. The amplitude–duration relationship for flare has also not been well-studied for the case of UCDs, where it may differ from that seen in earlier M dwarfs. We do not see a clear relationship between fitted amplitude and duration in our sample (and so cannot confirm any other observed amplitude–duration relationships, e.g. Hawley et al. 2014), though our estimates of amplitude and duration rely on a good fit to the Davenport template, which is significantly more challenging with messy ground-based data. Our purpose with these tests is to examine the decay time, amplitude, energy, and spectral type limits of our flare detection method, not to reproduce a realistic flare frequency distribution (FFD). To model a realistic flare distribution, we would need to consider that low energy flares are more common than high energy flares, as discussed in Section 5.3. However, then we would have significantly more low energy flares injected and recovered, making our sensitivity results more robust for small amplitudes/decay times, and less robust for high energy flares.

For each energy bin, treated independently, we then randomly select five flares and inject them into our differential light curves at a random observed time, taking care to ensure the ±5 data points surrounding the time of the flare’s peak are within 0.01 d (14.4 min). This prevents those five flares from occurring too close to the start or end of night or during gaps in our observations, which we would always miss. We allow flares to overlap to reflect the scenario we see often in real light curves, however as there are only ever a maximum of five flares in a light curve at once, the impact on our recovery results is limited.

We then run the automatic part of our flare detection method, described in Sections 3.3.1 and 3.3.3, on all injected light curves and record the recovered flares. We consider a flare to be recovered if the recovered time for its peak flux is within 1 FWHM_flare (of that flare) of the injection time. We allow for a more flexible recovery time because of the decision to allow flares to overlap, which can lead to more complex structures with multiple peaks. If two or more flares overlap within 1 FWHM_flare, they may be detected as only one flare. We assume if they occur further apart in time then they should be easily separable. This process is then repeated for each energy bin five times, making sure to use different flares from our artificial sample in each iteration.

The manual vetting step of our flare detection technique was not feasible for checking the approximately 23 000 flares we injected. This may introduce a bias into our results. It is more likely during manual vetting that we would conflate low-amplitude flares or flares on high photometric scatter targets with noise, therefore by removing this stage, we have potentially inflated our recovery rates for the lowest energy flares, or for the faintest target stars.

Fig. 4 shows the results of the injection-recovery tests, carried out for approximately 23 000 artificial flares. We are able to detect a significant proportion of flares with amplitudes above 1 per cent. Due to the short exposure times of SPECULOOS, we have the advantage of being very sensitive to flares with a short duration (with FWHM_flare < 5 min). However, we are more limited in detecting longer duration flares (with FWHM_flare ≥ 60 min) due to the day-night cycle, with typical uninterrupted observation windows of 4–8 h. There will also be some stars which have nights of observation less than 4 h, due to weather, or limited visibility, and these will be the most difficult targets on which to detect slowly decaying flares. However, we note that during the manual vetting (Section 3.3.4), we inspect the light curves with flares removed and we do not find any undetected high energy flares remaining in our sample. We also assess our limitations with the photometric RMS of our light curves. When the RMS of a global light curve exceeds ∼0.7 per cent (for 5-min binning), we begin to see a drop in detection efficiency, resulting in an increase in the minimum detectable amplitude. Only three targets in our data set exceed this RMS.

Figure 4.

Fraction of injected flares recovered as a function of flare amplitude and FWHM_flare using the automatic flare detection algorithm described in Section 3.3.

Open in new tabDownload slide

For a given flare amplitude and FWHM_flare, the flare energy will vary depending on the effective temperature and radius of the star (equation 9). The recovery fraction for different flare energies across all stars is shown in Fig. 5. However, the minimum detectable energy of recovered flares varies with spectral type, as shown in Fig. 6, which facilitates the detection of lower energy flares on our coolest dwarfs. For low-energy flares, we have two competing biases. The light curves for the later, and fainter, M and L dwarfs have higher photometric scatter, which makes it difficult to detect small-amplitude (and therefore lower energy) flares. However, as flares have a strong white light component (Namekata et al. 2017), there is an increased contrast between flares and the stellar spectra of red dwarfs, which should make lower energy flares easier to identify in our coolest stars (Allred, Kowalski & Carlsson 2015; Schmidt et al. 2019). SPECULOOS is also optimized to observe ultra-cool objects, and due to this specificity, we only see a very minimal increase in RMS with later spectral types.

Figure 5.

Fraction of injected flares recovered as a function of flare energy.

Open in new tabDownload slide

Figure 6.

Fraction of injected flares recovered as a function of flare energy and the object’s spectral type. We find that as we move to cooler stars we are able to detect lower energy flares. The upper and lower flare energy limits, determined by the amplitude and FWHM_flare limits of the artificial flare sample and the radius and T_eff of each star, are shown in orange. We are only impacted by the sample limits for spectral types earlier than M5 and later than L0.

Open in new tabDownload slide

The results from our artificial flare injection-recovery tests allow us to calculate the recovery fraction of our flare detection method, R(E), which reflects the low recovery rates of low energy flares. The average R(E) for all stars is shown in Fig. 5, as well as its spectral dependence in Fig. 6. As our sample is already small, we are unable to reduce detection-sensitivity effects by simply limiting our sample to flares with energies above a minimum recoverable energy (Paudel et al. 2018), or energies where a sufficient fraction of flares are recovered (Davenport 2016). Instead, we decided to calculate flare occurrence rates in Section 5.4 using the method in Jackman et al. (2021), which models the decline in sensitivity for low energy flares.

4 IDENTIFYING ROTATION PERIODS

With the previously identified flares and bad observations masked, we search for rotation periods in the 154 low-mass objects with more than 20 h of observation. We apply the Lomb–Scargle periodogram analysis (Lomb 1976; Scargle 1982) to our global light curves, binned every 20 min. The SSO light curves are not uniformly sampled, as we have gaps in our data, not only from the day/night cycle, but also from bad weather, masked flares, and changes to our observation strategy. Therefore, careful treatment of the resulting periodograms is essential to remove aliases.

We searched for periods between the Nyquist limit (twice the bin size) and the entire observation window. While we cannot be fully confident in periods greater than half their observation window, by expanding the period range we search, we can estimate a lower bound on long-period rotation. We then remove all peaks in the periodogram with a false alarm probability below 3σ (0.0027). Any target which has at least one significant Lomb–Scargle peak undergoes visual inspection of their periodogram and the corresponding phase-folded light curves for all possible rotation periods (all significant peaks in the periodogram). In addition, by comparing with periodograms of the time stamps, airmass, and FWHM, we can eliminate signals which arise from the non-uniform sampling or ground-based systematics.

We decided to apply a similar classification system as Newton et al. (2016) to our rotating objects. Classifying our rotators helps account for the difficulties in observing from the ground, both through atmospheric systematics, and through irregular, non-uniform sampling. When examining each light curve, we ask ourselves several questions:

• Is the period clearly visible in the phase-folded light curve?

• Was the object observed for long enough to span multiple periods?

• Is the frequency an alias of the ‘day signal’, seen as integer multiples in frequency space (periods of 1, 0.5, 0.33 d, etc.)?

• Is there a correlation with systematics (as seen in the airmass or FWHM periodograms)?

• Can the period be seen by eye in the un-phased light curve?

• If we observed this object with more than one telescope, does this period fit them all?

• Can we easily disentangle the ‘real’ period from its one day aliases?

• Is the amplitude of the periodic signal above the level of noise in the light curve?

If a rotator passes all the above criteria, then we class it as ‘A’. If it fails any of the above, but the rotation still seems likely, then we class it as a ‘B’ grade rotator. Most commonly, the ‘B’ rotators are convincing, but we do not observe multiple cycles, or we cannot easily choose between the period and its 1-d aliases. Any light curves for which we see some periodic structure, but we cannot easily determine a period, we class as ‘U’. This can result from broad peaks in the Lomb–Scargle periodogram, or multiple possible periods due to lack of observations, or a very low amplitude periodic signal that is comparable to the noise level. Because of these ambiguities in the period measurements, we do not attempt to place errors on our period estimates. If we do not detect any periodic signal or we cannot remove correlations with systematics (such as for very crowded fields), then we class the object as ‘N’. Finally, in addition to the classes defined in Newton et al. (2016), we add an extra ‘L’ grade. This is for the light curves where we see clear long-period rotation, but the period is similar or longer than the time window observed. For these objects the best we can do is to estimate the lowest possible period we measure, and acknowledge that there are large uncertainties on these period values.

5 RESULTS

The SSO data sample, along with the results of the flare and rotation analyses, are presented as supplementary material.

5.1 SPECULOOS-South flare sample

From the SSO data set described in Section 2, we identify 234 flares. We find that of our 154 unique targets, 78 are flaring (50 per cent). These flaring stars span the spectral type range from M4 (T_eff = 3160 K) to L0 (T_eff = 2313 K). Fig. 7 shows the spectral type distribution and proportions of the flaring objects identified in the SSO sample. From this figure, we see the proportion of flaring stars stays consistently above 60 per cent for objects of spectral type M5–M7. We see that the rate of flaring stars begins to decline around M8 (∼30 per cent), with no detected flares for any object beyond L0. The coolest flaring star we detect is a 2313 K, M9.6V object (which is rounded to L0 in Fig. 7). However, for the L dwarfs, we are limited by the small sample size, and therefore cannot make any conclusions about whether the fraction of flaring objects continues to decrease beyond late M dwarfs. Likewise, we do not have any objects earlier than M4; however, these earlier M dwarfs have been well-studied by TESS, Kepler, and MEarth. Despite our small sample of M4 dwarfs, we see an apparent reduction in activity for stars earlier than M5 that agrees with the previous results [∼30 per cent for M4 with TESS, Günther et al. (2020), 25 per cent for M5 with ASAS-SN, Martínez et al. (2019)]. Several authors find a steep rise in flaring fractions around spectral type M4 (West et al. 2004; Yang et al. 2017; Martínez et al. 2019; Günther et al. 2020). Specifically, West et al. (2004) analysed the H α emission of cool stars (which roughly correlates with flaring activity, Yang et al. 2017; Martínez et al. 2019) and determined that the fraction of active stars rose monotonically from spectral type M0 to M8, peaked at M8 (∼80 per cent active) and then declined to L4. Our results show an increase from M4 to a broad peak from M5–M7, followed by a decrease from M8 to L2.

Figure 7.

Top: histogram of objects in the SSO data sample as a function of spectral type. The flaring objects are marked in light blue. Bottom: fraction of flaring objects as a function of spectral type.

Open in new tabDownload slide

We probe the parameter space of high-flare-rate stars with low-to-mid energy flares. Due to the relatively short baselines of SPECULOOS observations, if a target does not flare frequently, then it is unlikely we would identify it as flaring. This also means we are less likely to detect the rarer, high-energy flares (Gershberg 1972; Lacy et al. 1976). In conjunction with the difficulties arising from the day-night cycle (typical night observations are 4–8 h), which complicates our detection of very slowly decaying flares (see Fig. 4), it is unlikely that superflares with energies >10³³ erg (Shibayama et al. 2013) will be detected by our survey. From studying the completeness of our flare sample, we also conclude it is unlikely that we will detect any flares below 10^29.5 erg for any spectral type.

5.2 Flare energies

In this flare sample, we recover flares ranging in energy from 10^29.2 to 10^32.7 erg, with a median of 10^30.6 erg. We do not detect any superflares, with energies between 10³³ and 10³⁸ erg (Shibayama et al. 2013).

From our flare injection-recovery analysis, we see that we have a spectral dependence in the flares we can retrieve. This dependence means that we will struggle to recover flares on the earliest M dwarfs below an energy of ∼10³¹ erg. As we move towards later spectral type objects, we are able to recover significantly higher fractions of lower energy flares, as shown in Fig. 6.

5.3 Flare frequency distribution

To generate FFDs, we extract only the objects for which we have at least three flares to obtain a good linear fit while including as many stars as possible. There are 31 objects in the SSO flaring sample (of 78 objects) that have at least three flares (the 31 objects have between 3 and 19 flares, with an average of 5 flares). We bin the energies per 0.1 erg and calculate the cumulative frequency of flares greater than that energy. In Fig. 8, we present the best-fitting power laws from equation (11). We also present the individual FFDs and line fits for each of the 31 objects in Appendix  A (Fig. A1), to demonstrate the diversity of power-law fits between stars of the same spectral type. We also performed the same FFD fitting using only the actual flare energies (unbinned) which did not change our results. We note that these results are not corrected for the detection efficiency as in Section 5.4 due to the complexity of applying it to the cumulative flare rates for low numbers of flares. Instead, we correct for the recovery rate in the next section, where we combine flares from multiple stars to yield a much larger data set.

Figure 8.

FFDs for every star with at least three flares. We find the best power-law fit to the cumulative flare rates against flare energy for every star, shown by straight lines. The spectral type of each star is indicated by its colour. The abiogenesis zones (Rimmer et al. 2018) are also shown for each star as the regions shaded in green (similarly to Günther et al. 2020), and the photosynthesis thresholds (Lingam & Loeb 2019a) are shown in orange.

Open in new tabDownload slide

By extrapolating the relationship between the flare energy and cumulative frequency for each star, we can predict the amount of energy that flares would deliver to the planets orbiting those stars. However, extrapolating the power–law relationship from our parameter space into the high-energy-flare regime is dangerously unreliable, as it is common to see breaks in these power–law relationships (Silverberg et al. 2016; Paudel et al. 2018). If we instead treat these linear fits as upper limits, then we can give estimates for the maximum frequency of high-energy flares. We note that the power laws also vary with the energy binning and ranges chosen for fitting.

We see a tentative trend that as the stars get cooler, the flares we detect become less energetic and less frequent. However, we see a large diversity of FFD profiles even within the same spectral type. This could be due to variation in stellar age or metallicity. Since we are biased towards detecting lower energy flares for cooler stars, it is possible we more significantly underestimate the rates of lower energy flares for our earlier M stars. This effect may result in steeper power-law slopes, further separating the mid-M dwarfs from the late-M and L dwarfs. However, we only have one M4, M9, and L0 object with at least three flares, and within the larger sample of M5–M7 stars, we see little distinction. The decline in flare rates for cooler spectral types could be partially explained by our decreased sensitivity for early-M stars. If only the infrequent, high energy flares are within reach, we only pick up the objects that flare more often. Conversely, since we are more sensitive to detecting lower energy flares on late-M and L dwarfs, we are able to identify stars with lower flaring frequencies. This effect explains the lack of low frequency flaring stars of early spectral types, but not the lack of high frequency flaring stars of later spectral types. As we have a small sample of flaring ultra-cool stars (see Fig. 7), we cannot confirm whether they, as a whole, flare less frequently. Paudel et al. (2018) similarly find that the flare rates for low-energy flares decrease as they move towards later spectral types, with L0 and L1 dwarfs having the lowest flaring rates. Critically, however, they find shallower slopes for cool stars [also in agreement with Gizis et al. (2013) and Mullan & Bais (2018)], implying that they have higher occurrence rates for the high-energy flares that are inaccessible to us. High frequency of these flares could be sufficient for a star to enter the abiogenesis zone. However, we are unable to confirm this trend, due to our small parameter space and the large uncertainties on our individual power-law slopes.

If we take a more conservative threshold of at least five flares to fit the power law, then we extract FFDs for only M5-M7 stars (for which we have a much larger flare sample). For this sub-sample, we find values for α in the range 1.2–2, roughly in agreement with Paudel et al. (2018), who measure 1.3–2 with their sample of 10 UCDs (for a smaller sample, but a similar spectral type range). Paudel et al. (2018) discuss various reasons for the variations in FFD slopes that are not dependent on age or spectral type, such as rotation, stellar spot coverage, and magnetic field topology.

5.3.1 Converting bolometric flare energies to U-band

5.3.2 Prebiotic chemistry

Here, we consider the laboratory work of Rimmer et al. (2018) in defining ‘abiogenesis zones’ around each of our potential planet hosts, outside of which it is unlikely a planet would receive enough energy for the following pre-biotic photochemistry scenario to occur. Specifically, the pathway considered in this work is the synthesis of ribonucteotides, as a precursor to ribonucleic acid (RNA; Patel et al. 2015; Sutherland 2017; Xu et al. 2018). By considering stellar flares as a mechanism for providing this UV energy, we can calculate the abiogenesis zones around the stars in our sample from their FFDs. Whether or not this chemistry is possible on a planet does not tell us if life has originated there, but instead, whether this mechanism can allow that planet to generate this first building block for RNA. Conversely, if a planet does not receive the necessary energy for this reaction, it does not rule out alternate pre-biotic pathways for the origins of life.

If we extrapolate the FFD power laws to E = 10³⁴ ergs, we have only one star that provides the necessary UV flux to reach the abiogenesis zone (see Fig. 8). This star is an M6 object. If we extrapolate the FFD power laws instead to E = 10³⁶ ergs, then we have 13 objects that reach the abiogenesis zone (with spectral types consisting of 1 M4, 1M8, 8 M6, and 3 M7). However, as our sample is confined to probing the low energy flare regime, we have to be careful to not over-extrapolate to high energies. Therefore, we do not extend the power-law predictions further than E = 10³⁴ ergs. Almost all the FFDs will eventually reach the abiogenesis zone, but at very high energy the uncertainties in our power-law fit would be too large to draw any conclusions.

5.3.3 Photosynthesis

5.4 Applying our sensitivity to the flare sample

To compare the average flare occurrence rates for different spectral types, we must account for the incompleteness of our sample. Due to the difficulty of detecting small amplitude flares above photometric scatter, we likely underestimate the frequency of low energy flares. This is often seen as a non-linear ‘tail-off’ from the expected power law in log-log space (equation 11).

Due to the small population of flares on the earliest and latest stars in our sample, we build up flare numbers by combining spectral types into the following five bins: M4–5, M6, M7, M8, and M9–L0. We calculate every star’s unique recovery fraction, R(E), based on the results of the injection-recovery tests. We bin flare energies into 12 logarithmically spaced bins from 10²⁸ to 10³⁴ ergs. For every star, we find a value for the recovery fraction in each energy bin by calculating the fraction of flares with energies in that bin that were recovered. We smooth the recovery fraction using a Wiener filter of three bins. We average the recovery fractions within each spectral type bin to get the recovery fraction of an ‘average’ star.

We then extract the observed flare occurrence rates from our flare sample. First, we isolate only the stars in each spectral type bin. We include all flaring and non-flaring stars in the SSO sample to avoid overestimating our flaring rates. For every energy bin, E, we sum the total number of flares with energies ≥E. Then we divide the total number of flares by the sum of the stars’ total observation times to produce flaring rates. In doing so, we make the assumption that 10 stars observed for 10 h is equivalent to 1 star observed for 100 h.

We fit equation (16) to our observed flare occurrence rates and energies using a Markov Chain Monte Carlo (MCMC) parameter optimization procedure. To implement MCMC, we use the emcee Python package (Foreman-Mackey et al. 2013). We used 32 walkers for 10 000 steps and the last 2000 steps to sample the posterior distribution. As in Ilin et al. (2019) and Jackman et al. (2021), we multiplied our errors by R(E)^−0.5 to account for larger uncertainties in recovering the smallest energy flares. The observed flare occurrence rates for each spectral type bin are presented in Table 1 and the results of the best-fitting power laws for each spectral type bin are shown in Fig. 9 and Table 2. We obtain similar values of α = 1.88 ± 0.05, 1.72 ± 0.02, 1.82 ± 0.02, 1.89 ± 0.07, and 1.81 ± 0.08 for M4–5, M6, M7, M8, and M9–L0, respectively. This work demonstrates that, on average, the power–law relationship does not change with spectral type in the mid- to late-M regime, despite the large variations within each spectral type. FFDs with similar gradients, but offset (different y-intercepts), implies that these spectral types produce similar relative proportions of high and low energy flares, but the coolest stars have lower rates of flares of all energies.

Figure 9.

The computed flare occurrence rates (equation 10), and their 1σ uncertainties, against bolometric flare energy for an average star of each spectral type bin. These flare rates have been corrected for our incompleteness at low flare energies and represent the ‘intrinsic’ flaring rate (equivalent to equation (16) at very high energies). The solid lines mark up to the highest energy flare for each spectral type bin, after that the ‘extrapolated’ region is marked by the dotted line.

Open in new tabDownload slide

Our results are consistent with three similar studies of the relation between flaring rate and energy for cool stars. As previously mentioned, Paudel et al. (2018) measured a range of α from 1.3 to 2.0 for 10 UCDs from K2 light curves. Raetz et al. (2020) compiled K2 light curves and derived a value for α of 1.83 ± 0.05 for FFDs of the coolest stars in their sample (with spectral types M5–M6). Gizis et al. (2017) found α = 1.8 for three UCDs. However, our results have noticeably lower values of α than several other works on mid-to-late M-dwarfs. From Kepler light curves, Yang et al. (2017) and Yang & Liu (2019) generated flare catalogues and found on average M dwarfs had power-law indices of α = 2.09 ± 0.10 (for fully convective stars) and α = 2.07 ± 0.35, respectively. For a sample of cool dwarfs observed by TESS, Medina et al. (2020) found α = 1.98 ± 0.02. Seli et al. (2021) specifically isolated TESS’s observations of stars similar in spectral type to TRAPPIST-1, for which they found α = 2.11. This work also modified the FFD for TRAPPIST-1 generated by Vida et al. (2017) to include recovery rate, updating their value of α from 1.59 to 2.03 ± 0.02. Lin et al. (2021) used EDEN and K2 light curves to analyse the flaring activity of the nearby active M dwarf Wolf 359, for which they found α = 2.13 ± 0.14. We could find a lower value of α if the incompleteness of the sample is not fully described by equation (16). We see a slightly earlier tail-off at low energies compared to the Jackman et al. (2021) model for all spectral types, which may be evidence of this. Alternatively, previous works may find higher values for α than this work as they have studied higher energy flares on UCD targets using space-based observations. We are limited in this study to only frequent, low energy flares due to our observing strategy. Therefore, it is possible that power laws may steepen at higher flaring energies.

This analysis agrees with the results of Section 5.3, which show that there is a decline in flaring rates as we move to the coolest stars. The injection-recovery tests demonstrate that if flares were present on the lowest-mass stars we would be able to detect them; therefore, these flares must occur too infrequently to be detected with our survey strategy.

5.5 Flares and rotation

From the SSO data sample, we recover 69 (24 A, 22 B, and 23 L) rotators, with periods ranging from 2.2 h to 65 d, as shown in Fig. 10. We also find 29 U and 60 N class objects. From here on, we refer to ‘rotators’ as only grades A, B, and L.

We find 41 targets for which we can detect both flares and rotation. Therefore, of the 78 flaring objects described above, 53 per cent have clear rotation. Alternatively, we detect flares on 59 per cent of our rotators. It appears that while the fraction of rotators across all spectral types stays consistently between 20–50 per cent, there are increasing proportions of slow rotators for the later M9 and L0 stars.

We see clearly that the very fast rotators are much more likely to flare than slow rotators. By comparing the rotation period in steps of 0.25 d, P, with the proportion of stars with rotation periods ≤P that flare (see Fig. 11), we demonstrate that the likelihood that we detect a star as flaring decreases as the rotation slows. In Fig. 11 we only include our 46 ‘A’ and ‘B’ grade rotators, as the periods of the ‘L’ (long period) grade rotators have very large uncertainties. In this sample, at least 74 per cent of fast-rotating stars, with P ≤ 2 d, flare. Comparatively, 59 per cent of all rotators flare, 63 per cent of ‘A’ and ‘B’ rotators flare, 42 per cent of stars with no detected rotation flare, and 50 per cent of all stars flare. We note that the majority (33/46) of ‘A’ and ‘B’ rotators have periods ≤2 d.

Figure 10.

Top: histogram of objects in the SSO data sample as a function of spectral type, divided by rotation class. The ‘rotators’ are shown in varying shades of blue, and the ‘non-rotators’ are shown in shades of orange. Bottom: fraction of different rotator classes for each spectral type.

Open in new tabDownload slide

Figure 11.

The proportion of stars with rotation periods ≤P that have detected flares, against rotation period P, in steps of 0.25 d. We only include our 46 ‘A’ and ‘B’ rotators, which range in period from 0.09 to 25.7 d. The colour scale shows the number of stars with rotation period ≤P.

Open in new tabDownload slide

5.5.1 Comparison with MEarth

Within the SSO Sample, 20 targets also appear in the Newton et al. (2018) MEarth-South rotation sample. Of these 20 objects, our results agree for eight, and they suggest periods not found by MEarth for another seven. For the eight objects in agreement: two have long period rotation measured by MEarth, with periods too long to be measured by SPECULOOS (for which we find no rotation); four have no detected rotation with either survey; and two have a similar measured period by both surveys (though one has only a tentative detection with the SSO). For the remaining 12 objects, five have clear, short periods detected by SPECULOOS but not by MEarth, and two have long-period, low-amplitude estimates given by SPECULOOS with no detection in MEarth. We discount the remaining five objects classified as U/N in both surveys, with missing or disagreeing periods. All of the rotation periods measured by the two surveys are shown in Fig. 12, and their overlapping observations are reported in Table 3.

Figure 12.

A comparison of the rotation periods measured in this SSO data sample with the rotation periods from MEarth (Newton et al. 2018). The objects observed only by SPECULOOS are in light blue and the objects observed only by MEarth are in light orange. The dark blue points are the measured period in SPECULOOS for objects that have also been observed by MEarth. The different rotation classes are indicated by the marker shapes, with the uncertain (U) rotation periods made transparent. Note that the L rotators, a class unique to SPECULOOS, provide only lower limit estimates of the periods, due to our shorter observing baselines.

Open in new tabDownload slide

The difference in rotation periods measured by the two surveys is likely a result of their different observing strategies. While SPECULOOS performs continuous monitoring of every target for 4–8 h each night over several weeks, MEarth cycles through multiple targets during a night, returning to each at 20–30 min cadences. This advantages SPECULOOS to observe very short, <5 h period rotators, and MEarth to measure very long, >50 d rotation periods. We also derive a long-period rotation estimate with SPECULOOS for two objects with small amplitudes that may not be evident in the MEarth data, due to limitations on the lowest mass stars. MEarth experiences a drop-off in recovery rate for stars with M < 0.2 M_⊙ (which includes all of the objects in our sample) likely due to systematics and the precipitable water effect (Newton et al. 2018).

The Newton et al. (2018) sample shows a clear dichotomy of fast rotators, with periods less than 10 d, and slow rotators, with periods greater than 70 d. We are unable to confirm this, as the majority of our sample have been observed for less than a 10-d span. For the clear slow rotators, which we observe for less than a full phase, we can only assign a long-period estimate (L). The apparent gap in rotation periods from our sample between 1 and 2 d is likely a bias from the difficulty in disentangling real rotation from 1-d aliases.

5.5.2 Comparison with TESS

Seli et al. (2021) study the relationships between age, activity, and rotation for a sample of 248 ‘TRAPPIST-1 analogues’ from 30-min TESS Full Frame Image (FFI) observations. They detect a total of 94 flare events on 21 stars. We compare their target catalogue to the SSO data sample and find 35 objects that appear in both. Of these 35 objects, we detect 19 as flaring, including five of the seven targets that Seli et al. (2021) identify as flaring. The five objects that both surveys identify as flaring stars are Gaia DR2 2331849006126794880, 2349207644734247808, 3200303384927512960 (as shown in Fig. 3), 4825880783419986432, and 5055805741577757824. The two flaring stars detected from TESS light curves, on which SPECULOOS does not detect flares are Gaia DR2 4967628688601251200 and 5637175400984142336. For 14 of the targets in this overlapping sample of 35 objects neither survey detects any flares. The 21 overlapping objects (excluding the 14 that have no flares detected by either study) are presented in Table 4.

For 42 stars in their sample, Seli et al. (2021) detect rotation periods. Out of the stars that appear in both UCD samples, Seli et al. (2021) find rotation periods for four objects that appear in our catalogue. We find rotation periods for 16 of the 35 objects (with rotation classes A, B, or L). Therefore, we can also compare both our flaring and rotation results for this handful of targets. For three of these objects, we detect an alias of the period extracted from TESS 30-min FFIs. For each of these three objects, we inspected the SSO light curves and periodograms. In each case, the rotation period measured by Seli et al. (2021) is the second most promising period, however, the peaks in the periodogram are very close in power, and the chosen period usually gives a better phase coverage. This is the reason for the ‘B’ grade we assign to rotators Gaia DR2 2349207644734247808 and 4967628688601251200. Gaia DR2 3200303384927512960 is assigned a grade of ‘A’, however, it is a unique case. For this target, we detect a rotation period of 0.334 d, while Seli et al. (2021) measure a period of 0.500 62 d. Rotation periods of 0.33, 0.5, or 1 d are special scenarios for ground-based observations where, over short baselines, we only see the same part of the phase on every observation night. Therefore, due to the limit imposed by the day-night cycle, we are unable to distinguish between these periods. The target for which we do not detect any rotation period, Gaia DR2 4971892010576979840, has an amplitude for periodic photometric variation of 2.1 mmag [measured by Seli et al. (2021)]. Upon inspection of this object’s SSO light curve, this periodicity is too small to detect above the level of photometric scatter, and residual systematics.

Since TESS is a space-based mission that does not contend with the day-night cycle or ground-based systematics, and they perform 27-d continuous monitoring, this provides a large number of phases for fast-rotating stars, which can counteract lower photometric precisions for red dwarfs. Therefore, we would expect rotation periods derived from TESS light curves to be more accurate than those from ground-based observations.

For their catalogue of stars, Seli et al. (2021) find that ∼8 per cent are flaring. This flaring fraction contrasts this work, where we detect flares on 28 of the 57 M7–M9 objects within the SSO sample (49 per cent). There may be several reasons for this discrepancy. First, Seli et al. (2021) detect flares from light curves with a low time resolution of 30 min. TESS is also not optimized for very red objects such as UCDs, therefore the photometric precision is lower than for SPECULOOS. Both of these factors inhibit a flare survey on TESS data from detecting low energy flares (which have small amplitudes and short durations) on UCDs. Instead, continuous 27-d monitoring of each TESS field allows Seli et al. (2021) to probe the rarer, high energy (even superflare) flare regime. This is demonstrated by their recovery fraction, which is less than 20 per cent for flares with energies below 10^31.5 ergs, compared to SSO recovery fractions of 50–70 per cent (depending on spectral type) for the same energy. Additionally, in this paper, we demonstrate that cooler stars flare less frequently at all energies, therefore it is likely that a survey focused on high energy flares on cool stars would detect significantly fewer flares. Finally, space-based data are less impacted than ground-based data by time-varying systematics due to the Earth’s rapidly changing atmosphere, therefore the risk of false positives is presumably lower.

6 DISCUSSION

6.1 Extension into the ultra-cool regime

The breadth of the SSO flare sample allows us to compare the flaring and rotating behaviour of mid-M stars with UCDs. We predominantly detect flares on M5–M7 stars (∼60–70 per cent flare), and decreasing proportions of flaring stars for both earlier and later M dwarfs. Since there is a spectral dependency to our flare detection (Fig. 6), it is likely that we are missing the lower energy flares on mid-M dwarfs, and therefore underestimating the proportion of flaring stars. On the other hand, these higher mass objects do produce rarer, high-energy flares more frequently than late-M dwarfs (Lacy et al. 1976). These high-energy flares should be easier to detect. Our results also agree with other surveys, with around ∼30 per cent of M4 stars flaring. It is possible that these surveys are limited in different ways for mid-M stars, which exist at the redder end of their samples, rather than at the bluer end of ours.

Lacy et al. (1976) showed that the flares on mid-M dwarfs were more frequent and lower energy than on early-M dwarfs. Hilton et al. (2010) found that this trend continued into the late-M dwarf regime (M6–M8). From our findings, we see little difference between spectral types up to M7. However, our results extend beyond M8, where we find that flares of all energies appear to become less frequent. From our injection-recovery tests, we conclude that if higher energy flares were present on the reddest stars in our sample, then we would detect them; therefore, the lack of high-energy flares for late-M and L dwarfs on the sampled time-scales of 0.8 to 8.7 d is likely real.

We see a relationship between flaring activity and rotation, especially for our very fastest rotators. The stars with rapid rotation periods of <5 d, are much more likely to flare than the rest of the SSO sample. At least 74 per cent of the stars with rotation periods faster than 2 d flare, compared to around 59 per cent of all rotating stars and 42 per cent of non-rotating stars. Therefore, we are more likely to detect flares on stars where we can determine a rotation period than stars we cannot. Possible explanations are that fast rotators flare much more frequently, more energetically, or are more likely to flare than stars with no detectable rotation period. Being unable to detect a rotation period may be a symptom that a star has a period too long to be measured by our survey, or the photometric variations on the rotation time-scale are too small to be seen, possibly due to small spot covering fractions, a low spot-photosphere contrast, or a particular orientation of spots (e.g. distributed at the poles). This contrast between rotators and non-rotators is further enhanced for the very fastest periods (P ≤ 2 d). Assessing the spectral distribution of our fastest rotators shows that almost all of the M4–M7 stars with periods less than 2 d flare. Comparatively, very few objects with spectral types later than M8 flare, fast-rotating or otherwise. Stars with fast rotation will have more magnetic energy available; therefore, these results are consistent with faster rotators exhibiting more active behaviours (Allred et al. 2015; Davenport 2016; Newton et al. 2017; Martínez et al. 2019; Yang & Liu 2019; Günther et al. 2020). Specific activity-rotation studies on cool stars, such as those carried out by West et al. (2015), Stelzer et al. (2016), Medina et al. (2020), and Raetz et al. (2020), have also found that fast rotating stars flare more frequently than slow rotators. Seli et al. (2021) found that ∼50 per cent of their fast-rotating (P < 5 d) late-M stars flared, compared to ∼70 per cent of our sample with the same periods.

In this work, we have compared results for SPECULOOS targets that also appear in the rotation study performed by Newton et al. (2018) on MEarth data, and in the activity-rotation study from Seli et al. (2021) on TESS observations. SPECULOOS and MEarth are two ground-based surveys focused on cool stars, but with very different observing strategies. SPECULOOS is a fast-cadence survey that targets UCDs for relatively short baselines of 4–8 d. MEarth, however, revisits its M dwarf targets only every 20–30 min, over a period of several months. With SPECULOOS’s light curves, therefore, we are able to study the short-term photometric variability of an object, such as its flaring and fast rotation (<10 d), whereas with MEarth’s light curves, we can monitor an object’s long-term variability, such as any slow rotation (>50 d) and stellar surface evolution. SPECULOOS and MEarth, however, are both ground-based surveys that have to contest with the challenges of observing through the Earth’s atmosphere and the loss of data from the day-night cycle. TESS, a space telescope, is able to perform continuous 27 d monitoring of each field, however UCDs are not its ideal targets. Therefore, TESS light curves of UCDs are mostly only available from the 30-min FFIs, and have generally lower precisions. SPECULOOS has the short exposure times and high precisions necessary to identify frequent low energy flares, however for rotation, the irregular data sampling and large gaps (limited by weather conditions and the observing cycle) make it difficult to identify the correct period above its aliases. For this reason, we recommend caution when interpreting the rotation results from this survey, and provide ‘rotation classes’ to help identify targets with rotation period uncertainties. TESS is only able to detect the rarer high-energy flares (with large amplitudes and slow decay times), however, its ability to provide continuous phase coverage of fast-rotating objects, and to observe up to hundreds of cycles, make it well suited for rotation studies. Comparisons such as these, between ground and space-based surveys, and between fast and slow-cadence observing strategies highlight that they are highly complementary. Leveraging the differences and the overlaps between different surveys allows us to fully explore the photometric activity parameter space of UCDs, and unearth more information about the subsequent impact on any potential planets they might host.

6.2 Planetary impact

Often used to put stellar flares into context, the largest solar flare ever directly observed on Earth was the ‘Carrington Event’ in 1859 (Carrington 1859; Hodgson 1859). This flare released an energy of 10³² erg and the associated CME hit the Earth’s magnetosphere, resulting in the largest geomagnetic storm on record. This storm disrupted telegraph systems and produced auroras reaching almost as far as the Equator. While our flare sample is relatively low-energy compared to the wealth of recorded flares on other stellar types, in our sample we detect 14 flares with energies greater than 10³² erg on 10 different stars. Characterizing these interactions between planets and their host stars, especially those which have the capability to initiate or destroy life, is vital to our understanding of planetary habitability.

6.2.1 Habitability of Earth-sized planets around UCDs

From calculating FFDs for the stars in our sample, we can place energetic bounds on whether flares can drive photosynthesis (Section 5.3.3) or provide their planets with enough UV flux for the pre-biotic chemistry described in Section 5.3 to occur. If we extrapolate the FFD power laws to E = 10³⁴ erg, we have only one M6 star that provides this necessary UV flux. However, there are several caveats to these results. As we see a non-linear ‘tail-off’ effect for low energy flares due to our detection limits, we likely underestimate low energy flares rates, resulting in shallower FFD power laws in Section 5.3. Therefore, without accounting for the completeness of our flare sample, this work only places upper bounds on the stars whose planets can reach the abiogenesis zone, or which can sustain photosynthesis. As we move towards cooler stars, we see that the frequency of flaring decreases (see Fig. 8). However, there are only one of each of the M4, M9, and L0 stars, with at least three flares detected. Due to our low detection efficiency for mid-M dwarfs, we are biased to detect only the most frequently flaring stars; therefore, it is likely that we would see an artificial increase in flare rate for M4 objects. In contrast, we should be more sensitive to flares of all energies on our late-M and L dwarf stars, so this bias does not explain the decrease in flare frequency for these stars.

By including our detection sensitivity in Section 5.4 and considering an ‘average’ star for each spectral type, we can conclude more general relationships between flare frequency and spectral type. While there is a wide diversity of flare frequencies for individual stars’ within a given spectral type, on average cooler stars have lower flaring rates. If this trend holds for a larger sample of UCDs, it may have severe consequences on the search for life around the very lowest mass objects. Similar conclusions have been reached by work on the flaring activity of TRAPPIST-1 (Ducrot et al. 2020; Glazier et al. 2020; Seli et al. 2021) and ‘TRAPPIST-1 analogues’ (Seli et al. 2021).

If the majority of late-M and L dwarfs have flaring rates too low to initiate the synthesis of ribonucleotides, then for life to originate it must happen through a different mechanism. Possible alternatives are surface hydrothermal vents (Rimmer & Shorttle 2019), aerial biospheres (Lingam & Loeb 2019b), impact-shock synthesis (Furukawa et al. 2015), or extraterrestrial delivery (Rimmer & Shorttle 2019; Krijt et al. 2017). We also find that none of our sample meets the requirement set out by Lingam & Loeb (2019a) to sustain an Earth-like biosphere, due to the significantly lower photosynthetically active radiation (with wavelengths of 400–750 nm) incident within UCD habitable zones. These results are consistent with the work of Mullan & Bais (2018) and Lingam & Loeb (2019a). However, we are considering an Earth-centric viewpoint, and it has been suggested that photosynthesis could proceed using deep-sea hydrothermal vents (Beatty et al. 2005) or at infrared wavelengths, where UCDs are brightest (Scalo et al. 2007; Takizawa et al. 2017; Claudi et al. 2021).

The activity of M dwarfs is predicted to decrease as they age (West et al. 2008, 2015; Paudel et al. 2019). Therefore, it is possible that earlier in these stars’ life cycles they would produce more energetic and frequent flares, with enough UV flux to trigger pre-biotic chemistry on their planets, or even for photosynthesis. The results in this paper only represent a snapshot of each object at its current point in time. We do not consider age in this work due to the huge uncertainties on stellar age estimates for these objects.

Modelling of the impact of stellar flares on modern Earth-like analogues by Tilley et al. (2019) shows that stars with at least 0.1 flare per day with energies above 10³⁴ erg will strip the ozone layers from any terrestrial planets they host. We do not extend our FFDs to energies above 10³⁴ erg as we are probing the low-energy flare regime; therefore, we do not consider ozone depletion in our analysis.

Here, we have only discussed two aspects of how flares could impact the habitability of planets around UCDs. In reality, we would need to apply a more holistic approach. There are a multitude of other factors (such as orbit dynamics, stellar variability, tidal locking, atmospheric composition, stellar age, etc.) that influence the environments on planets hosted by UCDs, and other possible pathways for life that have not been considered here.

7 CONCLUSIONS

By analysing the high-precision photometry produced by SPECULOOS, we can monitor the activity of late-M and UCDs. This allows us to probe further into the ultra-cool regime, expanding our understanding of the relationships between activity and rotation into a new parameter space. This work can be used in conjunction with previous and ongoing flare studies on earlier M dwarfs, such as those carried out on TESS, K2, and MEarth observations, to assess the diversity of M dwarfs and how well we can extrapolate the results of these studies to lower mass objects. We demonstrate the benefit of a large-scale flare study of UCDs, such as this, to extending existing flare catalogues.

Using our flare sample to make more general predictions about flaring activity as we move to very low-mass stars, we address two different scenarios in which flares may assist the origin and sustenance of life on planets they host. First, the quiescent UV flux of cool stars is too low for the synthesis of ribonucleotides, a major step in pre-biotic chemistry. Secondly, the visible light incident on planets around cool stars is also found to be too low for photosynthesis to occur. In both cases, while we find it is possible that flares on active M dwarfs could provide the necessary extra radiation, the frequency of flares of all energies appears to decrease with spectral type. These two considerations alone suggest UCD systems may not be favourable sites for abiogenesis, though other factors – known and unknown – certainly impact the likelihood of finding life in these systems as well.

SUPPORTING INFORMATION

ssodata.dat

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

The research leading to these results has received funding from the European Research Council (ERC) under the FP/2007–2013 ERC grant agreement no. 336480, and under the European Union’s Horizon 2020 research and innovation programme (grants agreements no. 679030 & 803193/BEBOP); from an Actions de Recherche Concertée (ARC) grant, financed by the Wallonia–Brussels Federation, from the Balzan Prize Foundation, from the BELSPO/BRAIN2.0 research program (PORTAL project), from the Science and Technology Facilities Council (STFC; grant no. ST/S00193X/1), and from F.R.S-FNRS (Research Project ID T010920F). This work was also partially supported by a grant from the Simons Foundation (PI Queloz, grant number 327127), as well as by the MERAC foundation (PI Triaud). PI Gillon is F.R.S.-FNRS Senior Research Associate. MNG acknowledges support from MIT’s Kavli Institute as a Juan Carlos Torres Fellow. BOD acknowledges support from the Swiss National Science Foundation (PP00P2-163967 and PP00P2-190080). This research is also supported by funding from STFC grants no. ST/S00193X/1 and ST/S00305/1. Finally, we would like to thank the anonymous referee for their helpful, detailed comments.

DATA AVAILABILITY

The data underlying this article is not currently public. The SPECULOOS-South Consortium will make all SPECULOOS-South Facility reduced data products available to the ESO Science Archive Facility following the regular Phase 3 process as described at http://www.eso.org/sci/observing/phase3.html.

REFERENCES

APPENDIX A: FLARE FREQUENCY DISTRIBUTIONS

In this appendix, we present the FFDs for each object with at least three flares. This includes the binned cumulative frequency against bolometric flare energy data and the best-fitting lines.

Figure A1.

The FFDs for each of the 31 objects with at least three flares. The cumulative flare rates against bolometric flare energies (binned every 10^0.1 ergs in log space) are shown by the points as well as the best-fitting power law for each star. The plots are ordered in increasing spectral type and the colour of the points and lines correspond to that object’s spectral type. The thresholds to reach the abiogenesis zone and perform oxygenic photosynthesis are shown for each object by the green and orange lines, respectively.

Open in new tabDownload slide

Author notes