The population of M dwarfs observed at low radio frequencies

Coherent low-frequency ($\lesssim 200$ MHz) radio emission from stars encodes the conditions of the outer corona, mass-ejection events, and space weather. Previous low-frequency searches for radio emitting stellar systems have lacked the sensitivity to detect the general population, instead largely focusing on targeted studies of anomalously active stars. Here we present 19 detections of coherent radio emission associated with known M~dwarfs from a blind flux-limited low-frequency survey. Our detections show that coherent radio emission is ubiquitous across the M~dwarf main sequence, and that the radio luminosity is independent of known coronal and chromospheric activity indicators. While plasma emission can generate the low-frequency emission from the most chromospherically active stars of our sample, the origin of the radio emission from the most quiescent sources is yet to be ascertained. Large-scale analogues of the magnetospheric processes seen in gas-giant planets likely drive the radio emission associated with these quiescent stars. The slowest-rotating stars of this sample are candidate systems to search for star-planet interaction signatures.

. Radio and optical properties of the 144 MHz detected M dwarfs. Sp. Type, d, S V , S V /S I , M , N , P rot , and T B correspond to spectral type, distance, Stokes V flux density, percentage of circularly polarised radio light, number of detections, number of exposures, stellar rotation, and 144 MHz brightness temperature, respectively. Known close (< 3 AU separation) binaries are identified by an asterisk, and the reported spectral type is for the total system. T B is a lower limit as it was calculated assuming that the emission emanates from the total photospheric surface. Potential ECMI sites are likely confined to much smaller areas, as observed on Jupiter 11,19 . Stellar rotation periods are not reported if unknown. Uncertainties on rotational periods are not listed if they are smaller than the last significant figure.

Name
Sp. Type from most main-sequence stars, which is often either hours-long with low or zero polarisation, or consists of short-duration ( 30 min), highly-polarised bursts 2,26,37 . The low-frequency emission from stars with low X-ray luminosities and long ( 2 d) rotation periods is more likely to be generated by ECMI, resembling the emission processes seen in the magnetospheres of sub-stellar objects such as ultracool dwarfs and gas-giant planets 11,19,26 . For the stars in our sample that are coronally active and rapidly rotating (e.g. AD Leo, DO Cep), plasma emission could plausibly generate the observed radio emission. With the first blindly formed sample of coherently-emitting low-frequency stellar systems, we can investigate the potential origins of the emission by identifying any trends with stellar properties. In particular, it is unclear how the coherent emission is related to the coronae of the stars. Radio emission associated with known radio stars is powered by high-energy supra-thermal electrons with energies of ∼ 10 keV−1 MeV 1,10 .  Figure 2. 144 MHz radio luminosity L ν,rad as a function of rotation period P rot (left panel), and radio surface flux density (L ν,rad /R 2 * , where R * is stellar radius) as a function of Rossby number R o (right panel), for ten of our M dwarf stellar systems. Chromospherically active, quiescent, and unknown chromospheric activity stellar systems are plotted as triangles, diamonds, and circles, respectively. We find no correlation between L ν,rad and P rot , or with L ν,rad /R 2 * and R o , with 95% confidence. The uncertainties in P rot and R o are smaller than the symbol size unless otherwise indicated. The error bars represent 1-σ uncertainty.
In stellar coronae, the electrons are accelerated by the chromospheric and coronal dissipation of magnetic fields 1,10 . This is evidenced by a tendency of chromospherically-active stars to be radio detected, a fact that has been used to probe the stellar dynamo 10 . In sub-stellar objects that lack a dominant corona, like brown dwarfs and gas giants, the electrons are thought to be accelerated at large distances from the central object due to either a breakdown of co-rotation between the object's magnetic field and its magnetospheric plasma, or an interaction with an orbiting satellite 12,26 .
To determine whether the coherent radio emission we detect is powered by processes in the corona or chromosphere, we compare the radio luminosity of our sample with known coronal and chromospheric activity indicators that have been used to understand the origin of gigahertz-frequency stellar emission. We find a conspicuous lack of correlation between the radio luminosity of our sample with well-known chromospheric activity indicators such as Hα luminosity, rotation period, and Rossby number as shown in Figure 2 and the Supplementary Information. The Kendall-τ rank correlation test rejects a relationship between the radio luminosity and rotation period, as well as between radio surface flux density and the Rossby number, at ≈ 95% significance.
The lack of correlation with the Rossby number is interesting because previous gigahertz-frequency observations have found a single trend of radio surface flux density with Rossby number that ranges from Sun-like G2 stars to L4 brown dwarfs 27,28,29 . Our low-frequency detected M dwarf sample spans over three orders of magnitude of Rossby number, with the largest Rossby numbers ever observed for radio-bright M dwarfs. Only substantially larger and slower-rotating G dwarfs have been detected in the radio with such large Rossby numbers 27 . We also do not predominantly detect rapid rotators, unlike gigahertz-frequency observations, with a constant detection rate for stars with rotation periods from 0.1 to 100 days (see Supplementary Information Section "Rotation and detection rate"). This is in contrast to the canonical expectation that radio bright stars are chromospherically active and rapidly rotating 27,29 . While Figure 2 likely has stars producing low-frequency emission via plasma and ECMI processes, the combination of the lack of correlation of radio surface flux density with Rossby number, the largest Rossby numbers for radio-bright M dwarfs, and the independence of our detection rate with stellar rotation period supports the idea that the low-frequency radio emission we detect is different than that commonly detected for M and ultracool dwarfs at gigahertz frequencies 27 .Additionally, our sample also shows no clear dependence on the coronal temperature and density traced by soft X-ray luminosity, as has also been shown for ultracool dwarfs 30 (see Supplementary Information Section "Güdel-Benz relationship").
It is also important to note that a high coronal temperature does not immediately imply that plasma emission is producing the low-frequency emission. For example, all 20 detections of CR Dra have an identical circularly polarised handedness, implying that the emission site is situated in a stable magnetic field arrangement 31 . Plasma emission would be expected to flip in circularly polarised handedness as it emerges from areas of complex magnetic geometry 2 . Therefore, even though CR Dra does have a large coronal base density, we suggest that we are observing ECMI from this star. In fact, for all of the stars of which we have multiple detections (Table 1), no star has flipped circular polarity. Consistent circularly-polarised handedness of persistent radio emission has also previously been observed on some active stars at gigahertzfrequencies, and has been used to support ECMI as the emission mechanism 2 .
The lack of a radio luminosity -stellar activity relationship, and the likelihood that the low-frequency emission is generated by ECMI in our quiescent stars, suggests that the radiation from the quiescent members of the population has similarities to auroral radio emission seen in sub-stellar objects. In particular, a coherent emission mechanism appears ubiquitous across the M dwarf main sequence and independent of whether the core of a M dwarf is partially or fully convective. Taken together, these two facts suggest that the electron acceleration site for some of our stars could be extrinsically sourced and that the emission could be a stellar analogue of the radio aurorae seen in brown dwarfs and gas-giant planets.
This sub-stellar auroral behaviour in M dwarfs is not unexpected. M dwarfs are known to retain strong magnetic fields over gigayear timescales 32 . As a result, their Alfvén surfaces likely extend to tens of stellar radii. Unlike the Sun, whose ECMI radiation is powered by processes within the corona, auroral processes in M dwarfs can transport energy from within the Alfvén surface back to the star to power intense ECMI emission. We note that one of our detections, WX UMa, provides direct evidence for a distant source site. WX UMa has a large-scale dipolar magnetic field with a strength of ∼ 3.5 kG at the stellar surface 32 (cyclotron frequency of ≈ 10 GHz) which means that the 120 MHz emission from WX UMa must originate at a distance of at least 5.5 stellar radii from the star 33 . Such a conclusion has the important implication that M dwarf coronae can be structurally similar to planetary magnetosphere in that they are dominated by a global dipole with a closed field topology.
Our sample likely consists of systems that are producing low-frequency emission via either plasma or ECMI processes. For the systems emitting via ECMI, there are two known mechanisms that power hours-long radio emission: co-rotation breakdown and interaction with a conducting and/or magnetised obstacle 26,34 . ECMI emission from co-rotation breakdown is powered by stellar rotation and has been used to explain the auroral emission from Jupiter to brown dwarfs and some ultra-cool dwarfs (M7 and later type) 26 . If some of the stars in our sample have magnetospheric parameters comparable to ultra-cool dwarfs, a simple empirical scaling of the radiated power with rotation period shows that this mechanism is only feasible for a small sub-sample of our stars which have rapid rotation (period 2 day) 34 (see also Supplementary Information Sections "Radio emission mechanisms" and "Rotation period and detection rate"). Therefore, the ECMI mechanism for stars in our sample with long rotation periods could be driven by star-exoplanet interactions 3,4 or by a hitherto unknown energy release mechanism, possibly from sites in the equatorial current sheets that extend into interplanetary space 35 . Exoplanet-driven ECMI has been theoretically modelled based on our understanding of Jupiter-Io interactions 4,8 . The radio flux density of our sample is consistent with theoretical expectations for interaction with the most common type of terrestrial planets found around M dwarfs 5 (see Supplementary Information Section "Radio emission mechanisms"). Although such planets are ubiquitous around M dwarfs, the small detected fraction could be explained by the angular beaming of ECMI emission combined with the requirement that only M dwarfs with surface magnetic field strengths larger than about 50 Gauss can generate emission at our observation frequency 36 (see Supplementary Information Section "Variability"). However, only follow-up studies to test if periodicity is evident in the radio emission from these stars will be able to validate such a suggestion.
Our sample can be interpreted along with previous gigahertz-frequency studies 2,26 , to arrive at a unified phenomenology of radio emission from M dwarfs. We suggest that M dwarfs, unlike the Sun, can harbour both star-like and gas-giant-like behaviour in their radio emission. Their coronae likely consist of short Xray bright coronal loops that flare occasionally, and a more tenuous and extended coronal gas supported by a large-scale magnetic field. The polarized short-duration bursts commonly detected at gigahertz frequencies likely originate in flaring coronal loops and therefore trace coronal and chromospheric flaring activity seen at other wavelengths 10 . On the other hand, the long-duration polarised emission we are detecting originates in the large-scale corona (outside the loops), from a quasi-continuous acceleration mechanism similar to that seen on Jupiter.
This unified model does not restrict the presence of ECMI emission to low frequencies. In highly-active M dwarfs with strong large-scale magnetic fields, the emission can extend to the gigahertz regime. Such emission has been recently detected in dMe-type flare stars with kilogauss-level magnetic fields 2,37 . Most previous efforts were restricted to highly active stars in which the distinction between continuous acceleration due to frequent coronal flaring and an extra-coronal process as proposed here is difficult to discern. The untargeted character of our detection strategy has allowed us to infer that such coherent emission processes are not isolated to chromospherically-active, rapidly-rotating M dwarfs. Our model can be tested via long-term monitoring of the sample in different spectral bands. Co-rotation breakdown in our fast-rotating stars will generate emission that is modulated at the rotation period of the star. On the other hand, exoplanet-induced emission must show similar modulation but over the orbital period of the putative planet, which can be confirmed by radial-velocity observations. In such a situation, it will be possible to assess the magnetic field strength of both the star and exoplanet, the latter being possible if the large-scale magnetic topology of the star is known and if the star-planet interaction is dipolar 4,38 . Therefore, some of our stars in this sample represent excellent candidates for searching for star-exoplanet interaction signatures.

ACKNOWLEDGEMENTS
The LOFAR data in this manuscript were (partly) processed by the LOFAR Two-Metre Sky Survey (LoTSS) team. This team made use of the LOFAR direction independent calibration pipeline (https: //github.com/lofar-astron/prefactor), which was deployed by the LOFAR e-infragroup on the Dutch National Grid infrastructure with support of the SURF Co-operative through grants e-infra 170194 e-infra 180169 39 43 ; ASTROPY, a community-developed core PYTHON package for astronomy 44 ; and NUMPY 45 .
Author Contributions: JRC initiated the LOFAR project that led to the discovery of the sources, conducted the crossmatching analysis, processed the extracted source visibilities, and wrote the manuscript. JRC and HKV developed the detection strategy. HKV led the theoretical interpretation of the detections and contributed substantially to the manuscript. TWS processed the radio data with LoTSS DR2 data reduction software developed by TWS, MJH, CT, WLW, FdG, JRC and other members of the LoTSS survey collaboration. BJSP helped develop the project and contributed to the manuscript. IED helped process individual stellar system datasets. JS and PNB processed the ELAIS-N1 deep-field data. HJAR is the principal investigator of the broader LOFAR Two-Metre Sky Survey. RJvW developed the extraction and re-calibration pipelines which optimized the calibration towards a target of interest. PZ provided comprehensive feedback on the manuscript. All authors commented on the manuscript.
Competing interests: The authors have no competing interests with respect to this manuscript. Data availability: LOFAR visibilities taken before 2020 are publicly available via the LOFAR Long Term Archive. The XMM-Newton data on LP 212-62, G 240-45, LP 169-22, and GJ 1151 are available through the XMM-Newton Science Archive (XSA). All other data used in the manuscript have been sourced from the public domain.

LoTSS circular polarisation data processing and source detection
Moderate resolution (20 ), wide-band (120 to 168 MHz) circular polarisation dirty images are standard products from the LoTSS-DR2 pipeline (described briefly in Sec. 5.1 of Shimwell et al. 21 and by Tasse et al. 46 ) which, to date, have been used to image approximately 1000 LoTSS pointings covering ≈20% of the Northern sky. The 8 hour LoTSS pointings are observed with 1 s integration and 12.5 kHz spectral resolution. The flux density scale of the images is corrected by matching sources in each field with NVSS and then assuming an average spectral index of sources. The median sensitivity of the Stokes V images is 80µJy beam −1 .
Each individual 20 Stokes V image was searched for ≥4σ sources using the Aegean source finder 47,48 , where σ is local-rms noise. To ensure a detection was reliable and not the product of a noise spike, any significant Stokes V source was required to be co-incident with a Stokes I source in a 6 image, as the corresponding Stokes I image has independent noise properties.
To determine the amount of Stokes I to V leakage, we measured the flux density at the brightest Stokes V pixel co-incident with the position of all bright (> 125 mJy), compact (< 6 largest angular size) Stokes I sources. The corresponding Rician distribution peaks at a Stokes I to V ratio of |V /I| = (0.14 ± 0.07)%. Based on this, we consider any ≥ 4σ Stokes V source with |V /I| > 1% as a reliable Stokes V detection. While we are confident that we have accurately characterised the leakage into Stokes V, to be conservative we exclude 2 regions around sources with a Stokes I flux density >125 mJy. Such a cut does not mean we are missing the bright end of the stellar system population, as evidenced by the resulting source counts distribution (see Supplementary Information Section "144 MHz source counts of M dwarf stellar systems").
Once significant and reliable Stokes V emission was identified, a post processing step was carried out to enable us to study the emission in more detail. After running the LoTSS-DR2 pipeline, all sources from a given LoTSS pointing were subtracted with the direction dependent calibration solutions, apart from sources in small region around the target of interest. The datasets were then phase shifted towards the region of interest and corrected for the LOFAR station beam attenuation. These steps are followed by three cycles of "TEC and phase" self-calibration on 10-20 sec timescales, which helped correct for any ionospheric distortions by solving for variations in total electron content (TEC). That step was followed by seven rounds of diagonal gain calibration on timescales of 10-60 min, depending on the amount of source flux density available. The "TEC and phase" solutions were pre-applied before solving for the diagonal gains. DPPP was employed for the calibration 49 and the images were made using WSClean 50 . The final products are small (<5 GB) datasets for each pointing that are calibrated in the direction of the target and with sources away from the target region subtracted.
All images were made with a robust parameter of −0.5 and the bandwidth/time-resolution was set to ensure that the source was 3σ in each image. The Stokes I and V images for all our detections are available in Supplementary Information Figure 1.
All LoTSS fields neighbouring to a confirmed detection were also inspected to test if the source was additionally detected in a different epoch. The source was considered undetected in a neighbouring field if the local rms was comparable to or better than the local rms of the field in which the source was detected, as expressed in Table 1. On average, this involved inspecting 2 to 3 fields per detection. CR Dra and GJ 625 are part of the ELAIS-N1 deep field which has 21 independent 8 hour LOFAR observations taken at approximately the same local sidereal time.

ROSAT and other X-ray flux measurements
The majority of our low-frequency detections have a counterpart in the Second ROSAT all-sky survey (2RXS) source catalogue 51 . We converted the ROSAT count rate to a 0.1 to 2.4 keV flux using the conversion factor 52,53 CF = (5.30HR1 + 8.31) × 10 −12 erg cm −2 counts −1 , where HR1 is the hardness ratio in the soft ROSAT band.
For the sources that did not have a counterpart in 2RXS, we searched if they had been observed by XMM-Newton or Chandra. Any flux measurements or upper limits were converted to the 0.1 to 2.4 keV energy band. 2MASS J14333139+3417472 was detected by Chandra several times in a survey of the Bootes field. We report the flux listed in the Chandra Source Catalog 54 .
XMM-Newton observed LP 212-62 serendipitously for ≈ 30 ks at ≈ 2.5 away from the pointing centre on 2007-11-11 (ObsID: 0503630201). We reduced the observation following standard routines 58,59 . The observation data were ingested using the XMM-Newton Science Analysis System (version 16.0.0), and background flares were filtered from the European Photon Imaging Camera data. A circular source region of radius ≈ 16 was used for spectral and timing extraction, with background spectra and timing extracted on the same CCD from a region with twice the area. Photon redistribution matrices and ancillary response files were constructed for the spectra and lightcurve for the pn, MOS1, and MOS2 cameras using standard calibration files.
The average 0.2 to 10.0 keV count rate for LP 212-62 during the observation was ≈ 0.08 count s −1 , but for ≈30 min the source flared to a count rate of ≈ 0.25 count s −1 . This flare was filtered before constructing the quiescent spectrum of LP 212-62. The quiescent spectrum of LP 212-62 was best fit by a single photoelectrically absorbed optically-thin plasma model (XSPEC 60 model apec) with a temperature of 0.41 ± 0.04 keV. Such a model provides a 0.2 to 2.4 keV flux of (9.5 ± 0.7)×10 −14 erg s −1 cm −2 , consistent with the 2RXS upper-limit at the position of LP 212-62.
GJ 1151 was also observed by XMM-Newton for ≈ 10 ks on 2018-11-01 (ObsID: 0820911301). The data reduction performed was identical to the process outlined for LP 212-62, with the exception that we discarded the pn data as GJ 1151's position on a chip edge meant the total effective area is not well calibrated by the standard pipeline. The average 0.2 to 10.0 keV count rate for GJ 1151 was ≈ 0.015 count s −1 , and we did not detect any X-ray flares. The spectrum was reasonably fit by an apec model with a temperature of 0.25 +0.09 −0.06 keV. The estimated 0.2 to 2.4 kev flux from such a model is (2.0±1.1) ×10 −14 erg s −1 cm −2 , consistent with the non-detection of GJ 1151 in an earlier 3 ks Chandra pointing 57 .
Finally, we observed G 240-45, LP 169-22, and 2MASS J09481615+5114518 with XMM-Newton for ≈ 15, 13, and 8 ks, respectively (ObsIDs: 0864480301, 086448101, 0864480201; PI: Callingham). Following identical data reduction procedures to those described above, we detect G 240-45 at ≈4 σ in the pn data. It is not detected in the individual MOS cameras. We do not detect any flares in the pn data. The spectrum is best fit by an apec model with a temperature of 0.3 +0.3 −0.1 keV. The estimated 0.2 to 2.4 kev flux from such a model is (1.1±0.9) ×10 −15 erg s −1 cm −2 . Therefore, even though G 240-45 is ≈3.5 times further away than GJ 1151, it has a similar X-ray luminosity.
2MASS J09481615+5114518 was detected in both pn and MOS data. There was no clear flares evident in the pn lightcurve. The spectrum is best fit by an apec model with a temperature of 0.14 +0.04 −0.05 keV. The estimated 0.2 to 2.4 kev flux from such a model is (2.1±0.8) ×10 −14 erg s −1 cm −2 .
LP 169-22 was not detected in our the XMM-Newton observation. While close to the X-ray detected star TYC 3451-1027-1, we derive a conservative three-sigma 0.2 to 2.4 kev flux upper-limit of 1.8 ×10 −14 erg s −1 cm −2 . This corresponds to an X-ray luminosity <2.3 ×10 26 erg s −1 , making LP 169-22 a quiescent candidate of our sample.

Radio source association with Gaia DR2
A benefit of wide-field radio surveys is that implicit biases inherent to targeted observations are bypassed, ensuring we observe a more complete sample of the stellar population outside of active stars and close binaries 10 . However, even with the LoTSS astrometric precision of ≈ 0.2 to 0.5 , a crossmatch of totalintensity (Stokes I) and Galactic Gaia Data Release 2 (DR2) sources is dominated by chance co-incidence matches 23 .
To form a reliable sample of Galactic Gaia-LoTSS counterparts, we restricted the crossmatching sample to LoTSS sources that were 4σ in Stokes V and 4σ in Stokes I, where σ is the local rms noise. Restricting our search to only circularly-polarised sources also suited our science goal of identifying coherent lowfrequency-emitting stellar systems. The source density of the Stokes V radio sky is orders of magnitudes lower than the Stokes I radio sky because > 99% of extragalactic sources are not circularly polarised. Only stellar systems, ultracool dwarfs, and pulsars are known to have circular-polarisation fractions > 10% 61,62,2 . We find ≈ 0.1 circularly-polarised source per sq. degree with a circularly-polarised fraction > 2%, compared to ≈ 770 Stokes I sources per sq. deg.
After identifying circularly-polarised LoTSS sources, we cross-matched their Stokes I position to all proper-motion corrected Gaia DR2 sources within 1 that had a parallax over error of / err 3. Since all the LoTSS pointings searched are at Galactic latitudes |b| > 20 • , the average Gaia DR2 source density for sources with / err 3 is ≈ 1300 per sq. degree. Therefore, we would expect the number of chancecoincidence associations between a Gaia Galactic source and the Stokes V LoTSS sample to be ≈ 0.05. As performed by Callingham et al. 23 , the chance-coincidence rate was also empirically confirmed by applying 10 -20 shifts to the Gaia DR2 positions and performing the crossmatch again. All of the sources presented in this paper were also observed to display variability incongruent with noise or proper motion shifts in neighbouring LoTSS fields. Therefore, all detections are considered reliable.
The above Gaia DR2 / err 3 requirement was to ensure we were selecting Galactic sources and to remain as unbiased as possible to the type of stellar systems detected. We did not restrict ourselves to crossmatching to only known M dwarfs so that we could systematically explore the type of sources that occupy the highly-circularly polarised Galactic radio sky. Hence, we have also detected other types of stellar systems such as RS Canum Venaticorum (RS CVn), BY Draconis (BY Dra), FK Comae Berenices (FK Com), and W Ursae Majoris (W Uma) variable stars, millisecond pulsars, and others. For this publication we have chosen to focus on the population of non-degenerate, main-sequence stellar systems that do not have known binary companions with < 1 AU separation. We have not detected any isolated, main-sequence star other than M dwarfs. For the area of sky we have observed, the non-detection of any isolated, main-sequence star is consistent with the number density of M dwarfs relative to G to K dwarfs stars within 50 pc 63 . We also note that we have not detected any previously identified radio-bright ultracool dwarfs, despite there being more than 20 candidates in our survey footprint 27 . This is likely because the majority of the searched ultracool dwarfs only have observed flux densities < 300 µJy at gigahertz frequencies. We are following up highly circularly polarised detections that have no optical counterparts.

Gaia DR2 M dwarf comparison sample
To derive the detection statistics for the M dwarf population presented Figure 1, we have to isolate a general representative population. We take our sensitivity horizon to be 50 pc, our most distant detection.
After selecting all Gaia DR2 sources within 50 pc that are in the area so far surveyed by LoTSS, we followed the quality cuts outlined by Kiman et al. 24 to form a clean sample of M dwarfs contained in Gaia DR2.
Firstly, we ensured that our sample contained stars with accurate astrometric data by requiring / err > 10 and that the number of visibility periods used were > 8. Secondly, to reliably select red stars we made a photometric quality cut in the red part of the Gaia band G RP ([630,1050] nm) by requiring a signal-to-noise ratio of > 10. This resulted in a final sample of 4518 stars with Gaia colour G − G RP > 0.9. All of our radio detected stars are contained in this final general sample.
Note that we did not make a quality cut on the Gaia "unit weight error" (UWE) as we wanted to remain sensitive to detecting unresolved binaries. The majority of the scatter around the main-sequence in the top panel of Figure 1 is because we have excluded this UWE quality cut.

Rossby number calculation
The Rossby number is a common parameter used to compare coronal and chromospheric activity strengths across both mass and rotation period ranges. To calculate the Rossby number for our detections, we require the rotation period of the star P rot and the convective turnover time τ conv . The stellar rotation period can be measured directly from photometric observations 55,56 . To estimate the convective turnover time, we used the empirical values derived using X-ray luminosity and V − K s colour for partially and fully convective active M dwarfs by Wright et al. 57 .

Circular polarisation sign convention
There are several instances in the collection, correlation, and reduction of LOFAR data where a signconvention decision can influence the final handedness of our Stokes V measurements. To add to further to the confusion, there are many competing definitions used in astronomy 88,89 . Pulsar astronomy commonly defines the sign of circularly-polarised light as left-hand circularly-polarised light minus right-hand circularly-polarised light. To ensure that the handedness of our Stokes V measurements are consistent with that definition 89 , we checked whether the Stokes V signs of pulsars observed in LoTSS were consistent with their circular polarisation handedness reported in the literature 61,62,89 . While we are integrating over many pulse profiles and comparing the pulsar's polarisation properties at a lower frequency than originally measured, pulsars are one of the few sources outside of stellar systems that can be strongly circularly polarised at low-frequencies.
In the LoTSS data we have reduced, we have detected significant circular polarisation from nine pulsars that are also in Han et al. 62 , Gould et al. 61 , and van Straten et al. 89 : B0917+63, B0025+33, B2210+29, B0751+32, B2315+21, B0655+64, B0823+26, B1737+13, and B0301+19. We exclude B0751+32 and B2315+21 from our analysis as their circularly-polarised handedness reverses between two of the references 61,62 . For the remaining seven pulsars, we are in agreement with the handedness reported in the literature, implying that our Stokes V sign convention is consistent with the definition applied in pulsar astronomy. The multi-wavelength properties pertinent to this manuscript for our sample are presented in Supplementary Table 1.
Supplementary Information Table 1. Multi-wavelength properties of the 144-MHz detected M dwarfs in increasing distance. M * , R * , L ν,rad , L X , L Hα /L bol , T B,plasma and Epoch correspond to stellar mass, stellar radius, 144 MHz radio luminosity, 0.2-2 keV X-ray luminosity, ratio of Hα to bolometric luminosity, the brightness temperature of fundamental plasma emission when a hot component 30 times the relative coronal temperature is injected, and the epoch of radio detection, respectively. If multiple detections are known, the earliest detection is listed. Stars with known close (< 10 day orbit) exoplanets are identified by an asterisk.

Radio emission mechanisms
There are several emission mechanisms that can produce the observed low-frequency radiation from M dwarfs. The principal radiation characteristics that are used to differentiate between the plausible emission mechanisms are the observed brightness temperature, luminosity, degree of circular polarisation, duration, bandwidth, and coronal temperature of the star. These emission characteristics for our sample are provided in Table 1 and Supplementary Table 1.
In this section we only outline coherent emission mechanisms since all of our sources display a circularly polarised fraction 50% and brightness temperature 10 12 K at 144 MHz. Therefore, incoherent gyrosynchrotron radiation can be confidently excluded 1 . More detailed information about the emission physics can be found in Vedantham 64 .

Plasma Emission
Plasma emission results from the conversion of Langmuir waves to transverse electromagnetic waves via non-linear processes 1 . The commonly encountered mechanism that generates Langmuir waves is the bump-on-tail and loss-cone instabilities. In a stellar corona they can form when impulsively heated plasma (T h ∼ 10 8 K) is injected into an ambient colder plasma (T c ∼ 10 6 K). The impulsively heated plasma is a proxy for the characteristic velocity of the beam that emits Langmuir waves via the bump-in-tail instability. The resulting radio emission can occur at the fundamental plasma frequency and its second harmonic.
We use the expressions derived by Stepanov et al. 65 to calculate the brightness temperature T B,plasma of the radio emission at the fundamental and the harmonic. The Langmuir wave spectrum was limited to wavenumbers ranging from those in resonance with the hot particles up to the Landau damping scale 8 . We also assume that the total energy density of the Langmuir waves is 10 −5 of the kinetic energy of the ambient plasma. Such an assumption corresponds to the canonical level at which the turbulence alters the dispersion in the background plasma, leading to a runaway collapse of the Langmuir wave packets 66,67 . We estimate the ambient coronal temperatures for our detected stars via T c = 0.11F 0.26 X , where F X is the X-ray surface flux calculated using a star's soft X-ray luminosity and the size of its photosphere 68 . Such a relation has been shown to hold for both active and quiescent M dwarfs. The coronal temperatures calculated for the X-ray detected stars range between 2 × 10 6 K (GJ 625) to 1 × 10 7 K (LP 259-39), with a mean of 5×10 6 K. For LP 169-22 and LP 259-39, which are not detected in X-rays, we assumed an X-ray luminosity equal to the three-sigma upper limit derived from our XMM-Newton observation and 2RXS, respectively. Finally, the hydrostatic density structure of the coronae of the stars was assumed to have a scale height 8 L n = 6 × 10 9 (T h /10 6 K)(R * /R ) 2 (M * /M ) −1 .
To calculate the brightness temperature of any potential plasma emission from our detected stars, we used Equations 15 to 22 of Stepanov et al. 65 and varied the temperature of the injected hot plasma to be up to 30 times the coronal temperature. Such temperature ranges are consistent with the observed temperatures of coronal loops 69  To make the comparison fair to the observed brightness temperatures reported in Table 1, we assume the size of the emitter is the stellar disk. Such an assumption is likely accurate to at least within a factor of two, as plasma source regions observed on the Sun generally do not exceed the size of the Solar photometric disk 72 . The estimated brightness temperature of fundamental plasma emission T B,plasma , when the injected heated plasma is equal to 30 times the coronal temperature, is provided in Supplementary Table 1.
We note that the assumptions required in this calculation imply that the estimated brightness temperature is only accurate to an order of magnitude for most of the sources. This is largely because it is unclear what is the total energy density of the Langmuir waves for each star. If we allow extreme values for the energy density, such as ≈ 10 −4 , fundamental plasma emission can produce the observed brightness temperature for about half of the sample. However, even in this extreme case, the plasma emission from the fundamental cannot achieve the measured brightness temperatures for our best quiescent candidates (GJ 625, GJ 450, LP 169-22, G 240-45, and GJ 1151) and more distant systems (LP 212-62, 2MASS J10534129+5253040, and 2MASS J14333139+3417472).
Plasma emission from the second harmonic can reach the observed brightness temperatures for all our detections. However, harmonic emission can not generate the high circularly polarised fraction detected for our sample. For example, Solar plasma harmonic emission has only been observed to have a circular polarised fraction < 20% 66 . For the most extreme configurations, where coronal loops occupy the entire stellar disk, the circularly polarised fraction of harmonic plasma emission is limited to ≈60% 73,74,75 . Only our detections of AD Leo and DO Cep have circularly polarised fractions < 60%. Furthermore, even if the entire stellar disk was occupied by coronal loops, it is likely that the emission from regions with opposite magnetic field orientation would result in substantially lower circularly polarised fraction. Therefore, we argue that plasma emission from the fundamental can plausibly generate the high observed brightness temperature and circularly polarised fraction for some, but not all, of our detections. Only the observed emission from AD Leo and DO Cep are potentially consistent with harmonic plasma emission due to their low circular polarised fractions.
We caution that the theoretical plausibility of plasma emission does not mean that it is the true mechanism. Plasma emission from flare stars is usually attributed to flaring coronal loops close to the surface 76 whereas, in our brightness temperature calculation, we have implicitly assumed that the whole surface of the star is filled with loops that are simultaneously flaring. There is no precedent for such a scenario. In addition, induced emission at the fundamental is significantly more efficient at higher frequencies (ν ∼ 1 GHz). Therefore, if disc-averaged T b ∼ 10 12 − 10 13 K from our sample is attributed to fundamental plasma emission, then we would expect to long-lived emission with T b 10 13 K at ν 1 GHz. Such emission has not been observed even on the most active stars 2 . Finally, if we are indeed observing plasma emission from short flaring loops, then we expect the polarity of the emission to change over time as loops with different orientation with respect to the line of sight come in and out of view. We do not observe this. For example, all 20 of CR Dra's low-frequency radio detections have a constant circularly-polarised handedness. The constant positive handedness of the circularly polarised light over 6.5 days of monitoring, taken in two observing blocks separated by a year 31 . For all our stars in which we have multiple detections, no star has ever flipped polarity (with follow-up observations occurring weeks and up to three years apart).
The formalism presented above also allows us to place some key coherent and low-frequency detections in the context of our detections. For example, Slee et al. 77 detected coherent, narrowband emission from Proxima Centauri that persisted for over two-days at ∼ 1.4 GHz. While the emission was 100% circularly polarised, the detected brightness temperature of 3 × 10 9 K is over three orders of magnitude lower than the brightness temperatures we measure for the majority of our sources. An application of the equations in Stepanov et al. 65 shows that the observed brightness temperature could readily be supplied by plasma emission. Additionally, Lynch et al. 9 recently detected coherent emission from the flare star UV Ceti at 154 MHz with the Murchison Widefield Array (MWA). While the emission achieved brightness temperatures of 10 13 − 10 14 K, the non-detection in total intensity only provided a lower limit to the circularly polarised fraction of 27%. Such emission is therefore also consistent with harmonic plasma emission, similar to our detections of AD Leo and DO Cep. However, we note that the MWA detections of UV Ceti appear consistent with flares lasting ∼30 min, while AD Leo and DO Cep have persistent emission lasting 4 hrs.

Electron-cyclotron maser emission
For low density, highly magnetised plasmas, coherent radio emission can be produced via the electron cyclotron maser instability (ECMI). There are several possible astrophysical configurations that can produce the population inversion necessary for the classical maser to operate. One common configuration is the unstable loss-cone distribution of electron energies set up in a stellar coronal loops 1 . In this case, the brightness temperature of a continuously operating maser is 78,8 where λ is the wavelength of radiation, R * is the stellar radius, L is the length-scale of the magnetic trap, and β = v/c, corresponding to the velocity of the electrons v divided by the speed of light c. ECMI emission will occur at a frequency dependent on the opening angle of the loss-cone and energy of the electrons. While the bandwidth of a single pulse of the maser will be small, the radiation will occur at a range of heights in the flux tube, generating broadband emission as high as 4:1 1 .
To check if the broadband emission of the observed brightness temperature can be generated by a compact coronal loop, we consider a heuristic model of a coronal loop of length L and a characteristic width of W = 0.1L 8 . The instantaneous bandwidth of the loss cone maser is ∆ν/ν ≈ β 2 α 2 where α is the opening angle of the loss cone 78 . Hence, to produce the observed broadband emission, the loop can be conceptually thought of as ν/∆ν maser sites each radiating with a brightness temperature given by Equation 1. For an observed flux density of 1 mJy at the distance of 20 pc (the median distance to our population), we arrive at a rough constraint on the loop length of L ∼ 25R * for characteristics values of β = α = 0.5. This shows that a loss cone maser operating in a compact coronal loop (L R * ) cannot account for the observed emission. Even in cases of the most active M dwarfs, coronal loops have only been observed to extend to ∼ 4 stellar radii 79 . We note that we could also have modelled this situation as many small-scale loops participating in the maser emission. However, that would likely result in a low circular polarisation fraction, or changing handedness of the polarisation. Additionally, if the observed emission was generated from coronal loops, we would have expected correlations with other coronal activity indicators, such as X-rays. Therefore, the observed emission is likely generated within a larger global magnetospheric trap.
As elaborated in the main text, an alternative way to drive ECMI radiation is through the breakdown of co-rotation. For example, the volcanic activity of Io in the Jovian system produces an equatorial plasma sheet that rotates with Jupiter. The plasma begins to lag behind the magnetic field of Jupiter at a radius of ≈6 Jovian radii, generating a magnetosphere-ionosphere coupling current that accelerates electrons into the neutral atmosphere of Jupiter 80 . The deposition of high-energy electrons at the poles of Jupiter set-up the population inversion necessary for ECMI to occur 81 . A similar particle system can also be established for M dwarfs, where continual flaring behaviour could deposit a torus of plasma around the star or, potentially, a volcanically active satellite.
Following the formulation of Nichols et al. 34 , ECMI radiation generated by the breakdown of co-rotation can only achieve the observed radio spectral luminosities (Supplementary Table 1) for stellar rotation periods 2 days. Excluding the systems in which we think plasma emission is a plausible origin of the low-frequency radiation, and stars without known rotation periods, the breakdown of co-rotation model can generate the observed radio luminosities for WX UMa, DG CVn, CR Dra, and (potentially) HAT 182-00605.
The physics of such systems share a similarity to the detection of ECMI radiation from rapidly rotating ultracool dwarfs 26 and V374 Peg 82 . Additionally, the linearly-polarised ≈1 GHz bursts identified by Zic et al. 37 from UV Ceti are likely produced by ECMI from the breakdown of co-rotation. However, the bursts that can be confidently associated with ECMI by Zic et al. 37 appear to last 20 min and spaced at least several hours apart. In comparison, the majority of our sample are observed to emit highly-circularly polarised light for at least eight hours. This may imply that we have a preferential line of sight to the emission site(s).
Finally, if breakdown of co-rotation was the sole mechanism driving the coherent radiation emission in our sample, we would expect rapidly-rotating systems to be preferentially detected. We show in Supplementary Information Section "Rotation period and detection rate" that that is not the case.
An alternative mechanism that can drive ECMI radiation is a Sub-Alfvénic star-planet interaction. Expressions for the radio power generated by a sub-Alfvénic interaction of a star and its exoplanet have been originally provided by Zarka et al. 4 . Here we use the slightly modified expressions of Saur et al. 83 and Turnpenney et al. 5 to compute the starward energy flux from a sub-Alfvénic interaction with an exoplanet: in Gaussian units, where α is the interaction strength 83 , v rel is the velocity of the stellar wind flow in the frame of the planet, B and ρ are the stellar magnetic field strength and wind density at the location of the planet, and R eff is the effective radius of the planetary obstacle to the stellar wind flow. The starward Poynting flux is converted to ECMI radiation with an efficiency rad . The flux density of the observed emission is then where Ω is the beam solid angle of the emitter, d is the distance to the stellar system, and ∆ν is the total bandwidth of emission. Because the emission has a cut-off at a frequency corresponding to the cyclotron frequency at the stellar surface, the bandwidth of emission can be approximated as ∆ν = eB * /(2πm e ) where B * is the surface field strength of the star, e and m e are the electron charge and mass respectively. Substituting for the starward energy flux, we get We note that because B is proportional to B * , to first order, the radio flux density is largely insensitive to the stellar magnetic field, although the frequency band over which a particular star radiates depends on the stellar magnetic field strength. This is consistent with the lack of dependence of the radio luminosity seen in our sample with proxies for the stellar magnetic field strength such as rotation rate and/or Rossby number. We examined the feasibility of the exoplanet-induced model in terms of the energy requirements as follows. We adopt a radiation efficiency of rad = 0.01, a conservative estimate 34 The theoretically expected Poynting flux can be computed from Equation 2. We adopt v rel = 500 km/s which is typical of winds driven by a thermal expansion of ∼ 10 6 K coronae. We computed the density ρ by conservatively adopting a base coronal density of n b = 10 7 cm −3 , a value that is an order of magnitude lower than the solar value, and assuming that radial expansion leads to a n b r −2 fall-off in density with distance r from the star in units of the stellar radius. We assume a Parker spiral configuration for the magnetic field with a surface field strength of B * = 100 G. We emphasise that the Alfvénic point only depends on the radial structure of the field strength and plasma density, rather than the azimuthal structure of the field. The radial profile of the angle θ in Equation 2 depends on the rotation rate of the star, the orbital period of the planet and the wind speed all of which will change substantially from system to system. As a benchmark value, we adopt θ = 0.2 which corresponds to the case of a planet orbiting at a distance of r = 10 (units of stellar radii) around a star with a 500 km/s wind and rotation period of 10 days. We assume a terrestrial Earth-like planet with R eff = 6500 km. Finally, we assume α = 1, which corresponds to the case of a planet with a highly conductive atmosphere (see Equation 17 of Saur et al. 83 ). The starward Poynting flux for our benchmark case is then Applying the aforementioned values to this equation allows us to recover the radio flux densities and luminosities listed in Table 1 and Supplementary Table 1. Therefore, ECMI radiation driven by a Sub-Alfvénic flux from star-planet interaction can produce the appropriate luminosity, duration, bandwidth, and brightness temperature for our sample, especially since we observe no correlation with other stellar activity markers. The model is particularly well suited to explain the origin of bright radio emission from sources that before to this study, and that of Vedantham et al. 8 , would not have been considered potential radio-emitters, such as GJ 625, GJ 450, GJ 1151, LP 169-22, G 240-45, and LP 212-62.
Finally, since we argue it is plausible that breakdown of co-rotation and star-planet interactions are driving the radio emission for some of our detections, beaming effects may substantially vary the flux density of our detections epoch to epoch. However, such intrinsic variations does not necessarily preclude observable correlations because stars with all activity levels are equally affected by beaming. In other words, the intrinsic variations affect quiescent and active stars in the same way so underlying correlations (or lack thereof) are not destroyed.

Güdel-Benz relationship
One test of whether our sample belongs to the class of canonical radio stars characterised by gigahertzfrequency observations 10 is to compare 144 MHz radio L ν,rad and 0.2-2.0 keV X-ray L X luminosity, as shown in Supplementary Information Figure 2. The correlation between quiescent gyrosynchrotron radio and soft X-ray luminosity, referred to as the Güdel-Benz relation where L X ∝ L 0.73 ν,rad , holds over 10 orders of magnitude in radio luminosity despite the fundamentally different emission mechanisms for the radio and soft X-ray radiation 84,85 . The standard model used to explain the Güdel-Benz relation posits that magnetic reconnection events in the stellar corona produce a population of non-thermal electrons that emit in the radio. This non-thermal population deposits some of its energy in the chromosphere of the star, producing an enhancement in soft X-rays as some of the newly ablated material concentrates in coronal loops 85,30 .
Since we know that a coherent emission mechanism is operating in our detected stars, we expect the sample should deviate from the established Güdel-Benz relation. As observed in Supplementary Information Figure 2, our sample significantly deviates from the established Güdel-Benz relation, with some stellar systems (e.g. GJ 1151, GJ 625, LP 212-62) departing from the relation by two to three orders of magnitudes. While it is not surprising that coherent auroral radio emission does not follow this canonical relationship 26 , it is noteworthy that the observed low-frequency radio luminosity does not show a clear dependence on the coronal temperature and density traced by the X-ray luminosity. For example, the seven detected stars at similar radio luminosities around ∼ 10 14 erg s −1 Hz −1 display a spread of over two orders of magnitude in their X-ray luminosity. Additionally, the systems do not establish a clear new relationship, suggesting that chromospheric activity has limited influence on the detected radio luminosity.
Such an inference is also supported by the constant low-frequency detection rate from spectral types M1.5 to M6.0, representing substantially different coronal temperatures and dynamo mechanisms. The detection of spectral types from M1.5 to M6.0 also implies that the reduction of coronal plasma heating efficiency can not resolve the violation of the canonical Güdel-Benz relation in our low-frequency detected population, as sometimes used to explain why ultra-cool dwarfs deviate from the Güdel-Benz relation 86,87 .
Finally, we note that there are underlying uncertainties when comparing our detections to the literature Güdel-Benz relation. The relation was established by observations between 5 and 8 GHz, while our observing frequency is significantly lower. Therefore, any spectral structure will impact the comparison. Supplementary Information Figure 2. Soft (0.2-2.0 keV) X-ray luminosity L X against radio luminosity L ν,rad . The literature data for the chromospherically-active stars used to derive the original Güdel-Benz relation are plotted as coloured triangles or squares 84,85 , with the best fit to the literature data indicated by the red line (L X ∝ L 0.73 ν,rad ). Some ultracool dwarfs with X-ray and radio detections are shown in yellow 30 . There are other ultracool dwarf systems with lower X-ray luminosities that are off the scale of this plot. Our 144 MHz radio sample (all with spectral types earlier than M7) is represented by black symbols, with triangles, diamonds, and circles indicating chromospherically active, quiescent, and unknown chromospheric activity, respectively. The chromospheric activity classifications are derived from Hα measurements (see Supplementary Information Figure 6). Non-detected X-ray sources are shown as 3-σ upper limits. While the radio luminosities plotted for the literature population are derived from 5 to 8 GHz observations, all of our detections have flat spectra.

Variability
The detection strategy of our survey involves searching LOFAR fields that have been synthesised over 8 hours of observation. Therefore, our survey can be biased towards detecting sources that show large, short duration bursts, or relatively constant emission that continues for multiple hours. We are confident that the emission we detect from our sources are of the latter type for two reasons. Firstly, if the sources were emitting bright (> 10 mJy) short (< 1 hr) bursts, they would have been previously detected by similar lower-sensitivity all-sky low-frequency Stokes V searches 7 . Secondly, we can form 144 MHz lightcurves for our sources, provided we integrate over the available bandwidth. The lightcurves for all our detected sources are shown in Supplementary Information Figures 3 and 4, where we have presented the absolute value of Stokes V emission. The sign of the circularly polarised emission can be inferred from Table 1. Figure 3. 144 MHz lightcurves for the detections listed in Table 1. Total intensity and absolute circular polarisation measurements are represented by red squares and black circles, respectively. To infer the sign of the circular polarisation, refer to Table 1. The title of each panel corresponds to the name of the stellar system. The dashed grey line represents zero level. Any measurement consistent with zero should be considered a non-detection. The error bars represent 1-σ uncertainty.

Supplementary Information
For all of our detected sources we can sample the lightcurves in three or more time bins. The flux density for each time sample was measured from images formed from integrating over 48 MHz bandwidth, and using the position of the star as a prior in forced-photometry fitting. The uncertainties represent 1-σ, and any uncertainty consistent with zero should be considered a non-detection. Different signal-to-noise properties influenced the variable sampling of the lightcurves, particularly the position of the source relative to the observing pointing centre or any non-Gaussian noise properties. Note that the data for the last ≈ 20% of LP 169-22's and middle of GJ 3729's observation are not included as the data was heavily impacted by radio frequency interference.
All of the lightcurves show highly circularly polarised emission that lasts at least 4 hours, and in the majority of cases, the entirety of the 8 hours of observation. Furthermore, no lightcurve shows characteristics of flare-like behaviour or a statistically significant flip in the handedness of the circular polarisation. We also made lightcurves with time spans of minutes to tens of minutes, but no stellar system showed temporal structure expected from Type II or III bursts. The lightcurves and the spectra of the detected sources will be studied in detail in a follow-up manuscript.
Finally, the variability of our sources between observing epochs, as listed in Table 1, could be used to test the favoured emission mechanism. For example, ECMI radiation will be beamed along the surface of a cone that will be modulated by the rotation of the star or the orbit of a putative satellite. However, the sampling currently available for most our stellar systems is not adequate to provide meaningful constraints considering the number of degrees of freedom in the radiation model, such as magnetic field and orbital orientation. Additionally, we note that our observing strategy is naturally biased to systems in which the geometry of the emission is preferential to our line of sight.

Rotation period and detection rate
Previous gigahertz-frequency studies of M and ultracool dwarfs have shown a clear increase in the fraction of radio detected stars with rotation period 27 -a star with less than ∼2 days rotation period is at least five times more likely to be detected at gigahertz-frequencies than a star rotating with period greater than ∼2 days. The relationship between rotation rate and radio detection is used as evidence that stellar rotation regulates the magnetic field strength and topology in the corona via the stellar dynamo 10,27 .
If the 144 MHz stellar system emission that we detect is driven by an amalgamation of breakdown of corotation and star-planet interactions, we would not expect a strong dependence of detection fraction with stellar rotation. While stars with rotation 2 day will be able to generate coherent low-frequency radio emission from the breakdown of corotation, stars with longer rotational periods far outnumber such rapid rotators, and these more plausibly generate their emission from a star-planet interaction. In Supplementary Information Figure 5 we show our 144 MHz detection rate with rotation period for all the M dwarfs reported by Newton et al. 55 that were in our survey footprint and at a distance less than 33 pc, the limiting distance of the study by Newton et al. 55 . We do not observe a trend with rotation period. While this figure is impacted by small number statistics, the number of detections implies we are sensitive to a change in detection rate of a factor of five or more.
We note that there are several important caveats to Supplementary Information Figure 5. In particular, while the work of Newton et al. 55 is the most comprehensive study of rotation in M dwarfs to date, it is inhomogeneously formed and almost certainly incomplete for stars with rotation periods 80 days. However, the lack of an observed trend between 0.1 and 80 days is robust to incompleteness. This incompleteness and inhomogeneity also implies that the absolute value of the detection rate should not be extrapolated to other surveys.  Figure 5. 144 MHz detection rate as a function of rotation period for the M dwarfs studied by Newton et al. 55 . While the sample produced by Newton et al. 55 is inhomogeneous and incomplete at periods 80 days, the lack of a trend in the detection rate with rotation is in contrast to what is observed by gigahertzfrequency studies 27 . The lack of an observed trend supports the conclusion that the 144 MHz emission being driven by a combination of breakdown of co-rotation and star-planet interactions. Due to incompleteness in the Newton et al. 55 sample, the absolute detection percentages should not be extrapolated to other surveys. The error bars represent 1-σ uncertainty.

Hα characteristics of the population
Chromospheric Hα emission is a tracer of magnetic activity, and is known to have a have a strong relationship with rotation in M dwarfs 55 . If the radio emission we have detected is related to the underlying dynamo, we could expect a relationship with Hα emission. As shown in Supplementary Information Figure 6, we see no correlation with Hα normalised by bolometric luminosity and 144 MHz radio luminosity − sources with equivalent Hα emission can vary by nearly 2 orders of magnitude in radio luminosity. Also, our sample spans the whole chromospheric activity range from non-detections in Hα to stars with Hα equivalent widths of ∼10Å.

144 MHz source counts of M dwarf stellar systems
Radio source counts N of a flux density-limited S sample can provide information about survey completeness and an estimate of the number of detections that a more sensitive survey could achieve 90 . The Stokes V and Stokes I differential source counts for our 19 M dwarf stellar systems detected at 144 MHz in LoTSS fields that were observed for 8 hours are presented in Supplementary Information Figure 7. Despite the sources being of Galactic origin, we perform a Euclidean Universe normalisation by multiplying the differential counts by S 2.5 to be consistent with previously presented radio source counts 90 . Supplementary Information Figure 6. 144 MHz radio luminosity against Hα luminosity by bolometric luminosity (left panel) and Hα equivalent width (right panel). The ratio of Hα luminosity to bolometric is relative to "inactive" stars defined by Newton et al. 55 . Upper-limits and error bars represent 1-σ uncertainty.
Supplementary Information Figure 7 demonstrates that we have not detected sources with total flux density S > 10 mJy in an 8 hour syntehsised image, and suggests our survey of nearby stellar-systems becomes significantly incomplete at ≈ 400 µJy. Such a value is consistent with our requirement that a Stokes V detection be 4σ and the fact the rms noise σ of a LoTSS field can vary between ≈60 and ≈120 µJy. Due to this rms noise variation per LoTSS field, it is also likely that the bin centred at ≈ 700 µJy is affected by incompleteness.
If we assume that the three bins at > 1 mJy are complete, we can estimate how many stellar systems a next-generation survey conducted with Square Kilometre Array (SKA1)-Low will detect. Assuming that the SKA1-Low survey will consist of 8 hour pointings at 148 MHz with a bandwidth of ≈ 50 MHz, a 4σ Stokes V detection will need to be at least 28 µJy based on current expected array sensitivity 91,92 . We predict that such a survey of the entire Southern sky (declination 0 • and below) with a Stokes V detection threshold of 28 µJy will achieve 7900±2500 detections of M dwarf systems. Note that this estimate assumes there is no transition in emission mechanism or type of emitter beyond 50 pc, which appears to be the sensitivity horizon of this work, and no difference in the volume density of detections when surveying the Galactic Plane. If we also include active binaries, such as RS CVns, the number of detections will at least double. Finally, the source counts also imply that we will have a total of 100±30 detections of M dwarf stellar systems when LoTSS finishes observing and processing the remaining area the Northern sky over the next ≈4 years. Supplementary Information Figure 7. 144 MHz Stokes V (left panel) and Stokes I (right panel) differential source counts for the 17 detected M dwarf stellar systems. Since we conducted for the initial source detection in circular polarisation, the Stokes V source counts are shifted to lower flux densities than the Stokes I source counts. The error bars represent 1-σ uncertainty.