# The population of M dwarfs observed at low radio frequencies

J.R. Callingham, H.K. Vedantham, T.W. Shimwell, B.J.S. Pope, I.E. Davis, P.N. Best, M.J. Hardcastle, H.J.A. Rottgering, J. Sabater, C. Tasse, et al.

► To cite this version:

J.R. Callingham, H.K. Vedantham, T.W. Shimwell, B.J.S. Pope, I.E. Davis, et al. The population of M dwarfs observed at low radio frequencies. Nature Astronomy, 2021, 5 (12), pp.1233-1239.
10.1038/s41550-021-01483-0. hal-03403123

HAL Id: hal-03403123

https://hal.science/hal-03403123v1

Submitted on 15 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
THE POPULATION OF M DWARFS OBSERVED AT LOW RADIO FREQUENCIES

J. R. CALLINGHAM,<fnref n="1" />,<fnref n="2" /> H. K. VEDANTHAM,<fnref n="2" />,<fnref n="3" /> T. W. SHIMWELL,<fnref n="2" />,<fnref n="1" /> B. J. S. POPE,<fnref n="4" />,<fnref n="5" /> I. E. DAVIS,<fnref n="6" />,<fnref n="2" />

P. N. BEST,<fnref n="7" /> M. J. HARDCASTLE,<fnref n="8" /> H. J. A. RÖTTGERING,<fnref n="1" /> J. SABATER,<fnref n="7" /> C. TASSE,<fnref n="9" />,<fnref n="10" /> R. J. VAN WEEREN,<fnref n="1" />
W. L. WILLIAMS,<fnref n="1" /> P. ZARKA,<fnref n="11" />,<fnref n="12" /> F. DE GASPERIN,<fnref n="13" /> AND A. DRABENT<fnref n="14" />

(Accepted to *Nature Astronomy*)

INTRODUCTORY PARAGRAPH

Coherent low-frequency ($\lesssim$ 200 MHz) radio emission from stars encodes the conditions of the outer corona, mass-ejection events, and space weather<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,<fnref n="4" />,<fnref n="5" />. 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<fnref n="2" />,<fnref n="6" />,<fnref n="7" />,<fnref n="8" />,<fnref n="9" />. 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<fnref n="1" />,<fnref n="10" />, 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<fnref n="4" />,<fnref n="11" />,<fnref n="12" /> 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.

<footnote n="1">¹Leiden Observatory, Leiden University, PO Box 9513, 2300 RA, Leiden, The Netherlands</footnote>

<footnote n="2">²ASTRON, Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, Dwingeloo, 7991 PD, The Netherlands</footnote>

<footnote n="3">³Kapteyn Astronomical Institute, University of Groningen, PO Box 72, 97200 AB, Groningen, The Netherlands</footnote>

<footnote n="4">⁴School of Mathematics and Physics, The University of Queensland, St Lucia, QLD 4072, Australia</footnote>

<footnote n="5">⁵Centre for Astrophysics, University of Southern Queensland, West Street, Toowoomba, QLD 4350, Australia</footnote>

<footnote n="6">⁶Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="7">⁷SUPA, Institute for Astronomy, Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ, UK</footnote>

<footnote n="8">⁸Centre for Astrophysics Research, University of Hertfordshire, College Lane, Hatfield AL10 9AB, UK</footnote>

<footnote n="9">⁹GEPI, Observatoire de Paris, Université PSL, CNRS, 5 Place Jules Janssen, 92190 Meudon, France</footnote>

<footnote n="10">¹⁰Department of Physics & Electronics, Rhodes University, PO Box 94, Grahamstown 6140, South Africa</footnote>

<footnote n="11">¹¹Station de Radioastronomie de Nançay, Observatoire de Paris, PSL Research University, CNRS, Université Orléans, OSUC, 18330 Nançay, France</footnote>

<footnote n="12">¹²LESIA, Observatoire de Paris, CNRS, PSL, Meudon, France</footnote>

<footnote n="13">¹³Hamburger Sternwarte, Universität Hamburg, Gojenbergsweg 112, 21029, Hamburg, Germany</footnote>

<footnote n="14">¹⁴Thüringer Landessternwarte, Sternwarte 5, D-07778 Tautenburg, Germany</footnote>
MAIN BODY

Low-frequency ($\lesssim 200$ MHz) radio observations are one of the few tools available that can probe the conditions of the outer coronae of stars$^{1,2}$. For example, magnetic eruptive events, such as coronal mass ejections (CMEs), are often accompanied by low-frequency emission as the ejected plasma traverses the corona into interplanetary space$^{13}$. In such situations it is expected that circularly polarised low-frequency emission would be generated by plasma mechanisms in stellar coronae$^{1}$, providing a measure of coronal density.

It has also been predicted that we should expect low-frequency radio emission to be generated by a satel-
lite interacting with its host star via the electron-cyclotron maser instability (ECMI)$^{4,18}$. The expected sign-
post of a star-planet interaction associated with M dwarfs – which are significantly less massive but more
magnetically-active than our Sun, and are the most common stars in our Milky Way$^{14,15,16,17}$ – is emission
similar to that produced by magnetospheric processes in gas-giants of the Solar System$^{4}$. Such emission
is expected to be low-frequency and highly circularly polarised ($\gtrsim 60\%$), with a brightness temperature
exceeding $10^{12}$ K$^{1,8,19}$.

While low-frequency radio emission provides unique information about the magnetic field and environ-
mental conditions around M dwarfs, particularly about high coronal layers, previous low-frequency blind
searches have lacked the required sub-mJy sensitivity necessary to achieve a detection$^7$. In addition, prior
low-frequency observations of M dwarfs have focused on the most chromospherically active stars, largely
because it has been hypothesised that the radio emission from CMEs might correlate with flare activity$^{2,9,20}$.
Such stars are unique and it is unclear if they are representative of the general low-frequency radio-emitting
population. With the LOw Frequency Array (LOFAR), we are currently conducting the LOFAR Two-meter
Sky Survey (LoTSS)$^{21}$ – the deepest, wide-field low-frequency radio survey performed to date. At 144 MHz,
we can achieve an optimal root-mean-square noise level of $\approx 60$ µJy per 8 hr exposure in circular polarisa-
tion (Stokes V)$^{21}$. With $\approx 20\%$ of the Northern sky processed, we have examined the LoTSS Stokes V maps
for $\ge 4\sigma$ sources and crossmatched those detections to sources in the Gaia Data Release 2 (DR2)$^{22}$ to search
for radio counterparts to stellar systems$^{23}$ (see methods).

In Table 1, we present 19 new low-frequency detections of M dwarfs, with only the LoTSS detection
of GJ 1151 having been presented in a previous publication⁸. The size of the sample represents over an
order of magnitude increase in the number of M dwarfs detected blindly at low frequencies, allowing us to
perform the first population study not biased by targeting strategies. As shown in Figure 1, the sample spans
the main-sequence spectral range M1.5 to M6.0, implying that we have detected systems with M dwarfs
that are both partially and fully convective. We detect 0.5 ± 0.1% of the stars in Gaia DR2 that are within
50 pc and which have spectral types M1.5 to M6.0. Each spectral type has an equal probability of detection,
within uncertainty, although we note that the Gaia catalogue begins to become significantly incomplete
for spectral types later than ~M7²⁴. Three of our detections (DG CVn, CR Dra, and GJ 3729) are known
binaries with less than 3 AU separation, which is consistent with the known binarity rate of M dwarfs²⁵.

All of our sources, except AD Leo and DO Cep, show a high degree of circular polarisation ($> 60\%$), high
brightness temperature ($T_B \gtrsim 10^{12}$ K), broadband emission ($\Delta\nu/\nu \sim 1$, where $\nu$ is frequency), and are
observed to be emitting for $\gtrsim 4$ hours in multiple fields. We demonstrate (see Supplementary Information
section “Radio emission mechanisms”) that depending on the coronal and chromospheric activity of the
star, plasma or ECMI mechanisms are the only processes that can generate the long-lived low-frequency
emission with such high brightness temperature and circularly polarised fraction$^8$. The detected emission
properties of our sample show a clear departure from the characteristics of gigahertz-frequency emission
<table><caption>**Table 1. Radio and optical properties of the 144 MHz detected M dwarfs.** Sp. Type, *d*, *S*<sub>V</sub>, *S*<sub>V</sub>/*S*<sub>I</sub>, *M*, *N*, *P*<sub>rot</sub>, and *T*<sub>B</sub> 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*<sub>B</sub> 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<sup>11,19</sup>. Stellar rotation periods are not reported if unknown. Uncertainties on rotational periods are not listed if they are smaller than the last significant figure.</caption><thead><tr><th>Name</th><th>Sp. Type</th><th>d<br>(pc)</th><th>S<sub>V</sub><br>(mJy)</th><th>|S<sub>V</sub>/S<sub>I</sub>|<br>%</th><th>M, N</th><th>P<sub>rot</sub><br>(days)</th><th>T<sub>B</sub><br>(×10<sup>12</sup> K)</th></tr></thead><tbody><tr><td>DO Cep</td><td>M4.0</td><td>4.00</td><td>-1.83±0.13</td><td>38±5</td><td>1, 1</td><td>0.41</td><td>0.7±0.1</td></tr><tr><td>WX UMa</td><td>M6.0</td><td>4.90</td><td>-1.52±0.15</td><td>96±4</td><td>3, 3</td><td>0.78</td><td>2.6±0.4</td></tr><tr><td>AD Leo</td><td>M3.0</td><td>4.97</td><td>-1.34±0.11</td><td>41±10</td><td>2, 2</td><td>2.23</td><td>0.4±0.1</td></tr><tr><td>GJ 625</td><td>M1.5</td><td>6.47</td><td>-1.20±0.19</td><td>75±15</td><td>6, 18</td><td>79.8±0.1</td><td>1.2±0.1</td></tr><tr><td>GJ 1151</td><td>M4.5</td><td>8.04</td><td>0.52±0.06</td><td>64±6</td><td>1, 4</td><td>125±23</td><td>1.4±0.2</td></tr><tr><td>GJ 450</td><td>M1.5</td><td>8.76</td><td>0.37±0.05</td><td>62±8</td><td>1, 2</td><td>23±1</td><td>2.0±0.3</td></tr><tr><td>LP 169-22</td><td>M5.5</td><td>10.49</td><td>-0.48±0.07</td><td>62±11</td><td>1, 3</td><td></td><td>4.4±0.9</td></tr><tr><td>CW UMa</td><td>M3.5</td><td>13.36</td><td>-1.21±0.11</td><td>61±4</td><td>1, 4</td><td>7.77</td><td>3.4±0.3</td></tr><tr><td>HAT 182-00605</td><td>M4.0</td><td>17.87</td><td>-0.59±0.08</td><td>76±7</td><td>1, 3</td><td>2.21</td><td>1.4±0.2</td></tr><tr><td>LP 212-62</td><td>M5.0</td><td>18.20</td><td>-6.43±0.17</td><td>90±6</td><td>4, 4</td><td>60.75</td><td>69.7±0.1</td></tr><tr><td>DG CVn*</td><td>M4.0</td><td>18.29</td><td>-0.57±0.06</td><td>91±9</td><td>2, 2</td><td>0.11</td><td>1.0±0.6</td></tr><tr><td>GJ 3861</td><td>M2.5</td><td>18.47</td><td>0.62±0.04</td><td>70±6</td><td>2, 2</td><td></td><td>1.6±0.2</td></tr><tr><td>CR Dra*</td><td>M1.5</td><td>20.26</td><td>8.12±0.14</td><td>92±8</td><td>20, 21</td><td>1.98</td><td>5.2±0.6</td></tr><tr><td>GJ 3729*</td><td>M3.5</td><td>23.57</td><td>-1.04±0.09</td><td>73±3</td><td>1, 3</td><td>13.59</td><td>3.4±0.5</td></tr><tr><td>G 240-45</td><td>M4.0</td><td>27.59</td><td>1.08±0.08</td><td>80±4</td><td>1, 3</td><td></td><td>40±5</td></tr><tr><td>2MASS J09481615+5114518</td><td>M4.5</td><td>36.17</td><td>1.76±0.09</td><td>96±4</td><td>1, 3</td><td></td><td>92±9</td></tr><tr><td>LP 259-39</td><td>M5.0</td><td>36.93</td><td>-0.41 ± 0.09</td><td>67±6</td><td>1, 1</td><td></td><td>20±4</td></tr><tr><td>2MASS J10534129+5253040</td><td>M4.0</td><td>45.19</td><td>-0.70±0.10</td><td>60±7</td><td>1, 3</td><td></td><td>15±3</td></tr><tr><td>2MASS J14333139+3417472</td><td>M5.0</td><td>47.84</td><td>0.68±0.11</td><td>61±5</td><td>1, 1</td><td></td><td>140±30</td></tr></tbody></table>

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<sup>2,26,37</sup>.

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<sup>11,19,26</sup>. 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.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Our low-frequency radio detections highlighted on a *Gaia* colour ($G - G_{RP}$)-magnitude ($M_G$) diagram (top panel), and our low-frequency radio stellar system detection statistics for *Gaia* DR2 (bottom panel). The *Gaia* source population plotted consists of a clean sample of M and L dwarfs that are within 50 pc and above a Galactic latitude of 20°24 (see methods). In the top panel, the LoTSS detected stellar systems are indicated by red circles, with red squares highlighting known close (< 3 AU separation) binaries. The colour of the density map represents the number of sources in each coloured pixel, as set by the colour bar in the top-right of the plot. The contour levels correspond to 25, 75, and 150 sources, respectively. In the bottom panel, our LoTSS-*Gaia* DR2 detection statistics for different *Gaia* colour bins are shown as black circles, with the mean of all our detections plotted as a red triangle. Upper-limits and error bars represent 1-$\sigma$ uncertainty.</figcaption></figure>

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¹,¹⁰.
<figure><img src="image_3.png" /></figure>

**Figure 2.** 144 MHz radio luminosity $L_{\nu,rad}$ as a function of rotation period $P_{\text{rot}}$ (left panel), and radio surface flux density ($L_{\nu,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_{\nu,rad}$ and $P_{\text{rot}}$, or with $L_{\nu,rad}/R_*^2$ and $R_o$, with 95% confidence. The uncertainties in $P_{\text{rot}}$ and $R_o$ are smaller than the symbol size unless otherwise indicated. The error bars represent $1-\sigma$ uncertainty.

In stellar coronae, the electrons are accelerated by the chromospheric and coronal dissipation of magnetic fields¹,¹⁰. 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¹⁰. 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¹²,²⁶.

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\alpha$ luminosity, rotation period, and Rossby number as shown in Figure 2 and the Supplementary Information. The Kendall-$\tau$ 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 $\approx 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²⁷,²⁸,²⁹. 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<sup>27</sup>. 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<sup>27,29</sup>. 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<sup>27</sup>. 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<sup>30</sup> (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<sup>31</sup>. Plasma emission would be expected to flip in circularly polarised handedness as it emerges from areas of complex magnetic geometry<sup>2</sup>. 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 gigahertz-frequencies, and has been used to support ECMI as the emission mechanism<sup>2</sup>.

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<sup>32</sup>. 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, WXUMa, provides direct evidence for a distant source site. WXUMa has a large-scale dipolar magnetic field with a strength of ~ 3.5 kG at the stellar surface<sup>32</sup> (cyclotron frequency of ≈ 10 GHz) which means that the 120 MHz emission from WXUMa must originate at a distance of at least 5.5 stellar radii from the star<sup>33</sup>. 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<sup>26,34</sup>. 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)<sup>26</sup>. 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 $\lesssim 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 X-ray 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.

## REFERENCES

1. Dulk, G. A. Radio emission from the sun and stars. ARA&A **23**, 169–224 (1985).

2. Villadsen, J. & Hallinan, G. Ultra-wideband Detection of 22 Coherent Radio Bursts on M dwarfs. ApJ **871**, 214 (2019).
3. Lazio, W., T. Joseph *et al*. The Radiometric Bode's Law and Extrasolar Planets. *ApJ* **612**, 511–518 (2004).

4. Zarka, P. Plasma interactions of exoplanets with their parent star and associated radio emissions. *Planet. Space Sci.* **55**, 598–617 (2007).

5. Turnpenney, S., Nichols, J. D., Wynn, G. A. & Burleigh, M. R. Exoplanet-induced Radio Emission from M dwarfs. *ApJ* **854**, 72 (2018).

6. Crosley, M. K. *et al*. The Search for Signatures of Transient Mass Loss in Active Stars. *ApJ* **830**, 24 (2016).

7. Lenc, E., Murphy, T., Lynch, C. R., Kaplan, D. L. & Zhang, S. N. An all-sky survey of circular polarization at 200 MHz. *MNRAS* **478**, 2835–2849 (2018).

8. Vedantham, H. V., Callingham, J. R., Shimwell, T. & Pope, B. J. S. Coherent radio emission from a quiescent red dwarf indicative of star-planet interaction. *Nature Astronomy* **4**, 577–583 (2020).

9. Lynch, C. R., Lenc, E., Kaplan, D. L., Murphy, T., & Anderson, G. E. 154 MHz Detection of Faint, Polarized Flares from UV Ceti. *ApJL* **836**, L30 (2017).

10. Güdel, M. Stellar Radio Astronomy: Probing Stellar Atmospheres from Protostars to Giants. *ARA&A* **40**, 217–261 (2002).

11. Zarka, P., Cecconi, B. & Kurth, W. S. Jupiter's low-frequency radio spectrum from Cassini/Radio and Plasma Wave Science (RPWS) absolute flux density measurements. *Journal of Geophysical Research (Space Physics)* **109**, A09S15 (2004).

12. Zarka, P. *et al*. Jupiter radio emission induced by Ganymede and consequences for the radio detection of exoplanets. *A&A* **618**, A84 (2018).

13. Wild, J. P. & McCready, L. L. Observations of the Spectrum of High-Intensity Solar Radiation at Metre Wavelengths. I. The Apparatus and Spectral Types of Solar Burst Observed. *Australian Journal of Scientific Research A Physical Sciences* **3**, 387 (1950).

14. Cohen, O. *et al*. The Interaction of Venus-like, M-dwarf Planets with the Stellar Wind of Their Host Star. *ApJ* **806**, 41 (2015).

15. Blackman, E. G. & Tarduno, J. A. Mass, energy, and momentum capture from stellar winds by magnetized and unmagnetized planets: implications for atmospheric erosion and habitability. *MNRAS* **481**, 5146–5155 (2018).

16. Khodachenko, M. L. *et al*. Coronal Mass Ejection (CME) Activity of Low Mass M Stars as an Important Factor for The Habitability of Terrestrial Exoplanets. I. CME Impact on Expected Magnetospheres of Earth-Like Exoplanets in Close-In Habitable Zones. *Astrobiology* **7**, 167–184 (2007).

17. Cauley, P. W., Shkolnik, E. L., Llama, J. & Lanza, A. F. Magnetic field strengths of hot Jupiters from signals of star-planet interactions. *Nature Astronomy* **3**,1128-1134 (2019).

18. Lazio, T. J. W., Shkolnik, E., Hallinan, G. & Planetary Habitability Study Team. Planetary Magnetic Fields: Planetary Interiors and Habitability. Tech. Rep. (2016).

19. Zarka, P. Auroral radio emissions at the outer planets: Observations and theories. *J. Geophys. Res.* **103**, 20159–20194 (1998).

20. Crosley, M. K. & Osten, R. A. Low-frequency Radio Transients on the Active M-dwarf EQ Peg and the Search for Coronal Mass Ejections. *ApJ* **862**, 113 (2018).

21. Shimwell, T. W. *et al*. The LOFAR Two-metre Sky Survey. II. First data release. *A&A* **622**, A1 (2019).

22. Gaia Collaboration *et al*. Gaia Data Release 2. Summary of the contents and survey properties. *A&A* **616**, A1 (2018).

23. Callingham, J. R., Vedantham, H. K., Pope, B. J. S., Shimwell, T. W. & the LoTSS Team. LoTSS-HETDEX and Gaia: Blind Search for Radio Emission from Stellar Systems Dominated by False Positives. *Research Notes of the American Astronomical Society* **3**, 37 (2019).

24. Kiman, R. *et al*. Exploring the Age-dependent Properties of M and L Dwarfs Using Gaia and SDSS. *AJ* **157**, 231 (2019).

25. Shan, Y., Johnson, J. A. & Morton, T. D. Measuring the Number of M dwarfs per M dwarf Using Kepler Eclipsing Binaries. *ApJ* **813**, 75 (2015).

26. Hallinan, G. *et al*. Confirmation of the Electron Cyclotron Maser Instability as the Dominant Source of Radio Emission from Very Low Mass Stars and Brown Dwarfs. *ApJ* **684**, 644–653 (2008).

27. McLean, M., Berger, E. & Reiners, A. The Radio Activity-Rotation Relation of Ultracool Dwarfs. *ApJ* **746**, 23 (2012).

28. Slee, O. B. & Stewart, R. T. Stellar radio luminosity and stellar rotation. *MNRAS* **236**, 129–138 (1989).
29.Stewart, R. T., Innis, J. L., Slee, O. B., Nelson, G. J. & Wright, A. E. A Relation Between Radio Luminosity and Rotation for Late-Type Stars. *AJ* **96**, 371 (1988).

30.Williams, P. K. G., Cook, B. A. & Berger, E. Trends in Ultracool Dwarf Magnetism. I. X-Ray Suppression and Radio Enhancement. *ApJ* **785**, 9 (2014).

31.Callingham, J. R. *et al*. Low-frequency monitoring of flare star binary CR Draconis: Long-term electron-cyclotron maser emission. *A&A* **648**, A13 (2021).

32.Morin, J. *et al*. Large-scale magnetic topologies of late M dwarfs*. *MNRAS* **407**, 2269–2286 (2010).

33.Davis, I. E. *et al*. Large closed-field corona of WX UMa evidenced from radio observations. *A&A* **650**, L20 (2021).

34.Nichols, J. D. *et al*. Origin of Electron Cyclotron Maser Induced Radio Emissions at Ultracool Dwarfs: Magnetosphere-Ionosphere Coupling Currents. *ApJ* **760**, 59 (2012).

35.Linsky, J. L. *et al*. Radio Emission from Chemically Peculiar Stars. *ApJ* **393**, 341 (1992).

36.Vidotto, A. A., Feeney, N. & Groh, J. H. Can we detect aurora in exoplanets orbiting M dwarfs? *MNRAS* **488**, 633–644 (2019).

37.Zic, A. *et al*. ASKAP detection of periodic and elliptically polarized radio pulses from UV Ceti. *MNRAS* **488**, 559–571 (2019).

38.Strugarek, A. Assessing Magnetic Torques and Energy Fluxes in Close-in Star-Planet Systems. *ApJ* **833**, 140 (2016).
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³⁹. The LoTSS direction dependent calibration and imaging pipeline (http://github.com/mhardcastle/ddf-pipeline/) was run on compute clusters at Leiden Observatory and the University of Hertfordshire, which are supported by a European Research Council Advanced Grant [NEWCLUSTERS-321271] and the UK Science and Technology Funding Council [ST/P000096/1]. The Jülich LOFAR Long Term Archive and the German LOFAR network are both coordinated and operated by the Jülich Supercomputing Centre (JSC), and computing resources on the supercomputer JUWELS at JSC were provided by the Gauss Centre for Supercomputing e.V. (grant CHTB00) through the John von Neumann Institute for Computing (NIC). This work was performed in part under contract with the Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute. JRC thanks the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) for support via the Talent Programme Veni grant. PNB and JS are grateful for support from the UK STFC via grant ST/R000972/1. RJvW acknowledges support from the Vidi research programme with project number 639.042.729, which is financed by the Netherlands Organisation for Scientific Research (NWO). AD acknowledges support by the BMBF Verbundforschung under the grant 05A17STA. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France, and NASA's Astrophysics Data System. This work has also made use of TOPCAT⁴⁰, the IPYTHON package⁴¹; SciPy⁴²; MATPLOTLIB, a PYTHON library for publication quality graphics⁴³; ASTROPY, a community-developed core PYTHON package for astronomy⁴⁴; and NUMPY⁴⁵.

**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.

**Code availability:** The important codes used to analyse the data are available at the following sites: WS-Clean (https://gitlab.com/aroffringa/wsclean) and DDF pipeline (https://github.com/mhardcastle/ddf-pipeline).

**Correspondence:** Correspondence and requests for materials should be addressed to JRC (email: jcal@strw.leidenuniv.nl).
# 1. METHODS

## 1.1. 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.²¹ and by Tasse et al.⁴⁶) 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⁻¹.

Each individual 20" Stokes V image was searched for $\ge 4\sigma$ sources using the Aegean source finder⁴⁷,⁴⁸, where $\sigma$ 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 \pm 0.07)\%$. Based on this, we consider any $\ge 4\sigma$ 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⁴⁹ and the images were made using WSClean⁵⁰. 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 $\ge 3\sigma$ 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.
<figure><img src="image_4.png" /></figure>

**Supplementary Information Figure 1.** Total intensity (Stokes I) and circularly polarised (Stokes V) (insets) 144 MHz images of the detected M dwarf systems. The red circle in each image represents the position of the Gaia DR2 counterpart at the time of the LOFAR observation. The radius of the circle is 1". The intensity scale of each image is different, with dark-blue representing positive flux density. All detections are ≥ 4σ in both Stokes I and V. The Stokes I images are deconvolved and have a 5" point-spread function. The Stokes V images are not deconvolved and have a 20" diameter point-spread function. The horizontal and vertical axes of each image are right-ascension and declination, respectively.
## 1.2. 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⁵¹. We converted the ROSAT count rate to a 0.1 to 2.4 keV flux using the conversion factor⁵²,⁵³ $CF = (5.30\text{HR1} + 8.31) \times 10^{-12} \text{ erg cm}^{-2} \text{ 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⁵⁴.

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⁵⁸,⁵⁹. 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⁻¹, but for ≈ 30 min the source flared to a count rate of ≈ 0.25 count s⁻¹. 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 photoelectronically absorbed optically-thin plasma model (XSPEC⁶⁰ 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⁻¹⁴ erg s⁻¹ cm⁻², 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⁻¹, and we did not detect any X-ray flares. The spectrum was reasonably fit by an apec model with a temperature of 0.25<sup>+0.09</sup><sub>-0.06</sub> keV. The estimated 0.2 to 2.4 keV flux from such a model is (2.0±1.1) × 10<sup>-14</sup> erg s<sup>-1</sup> cm<sup>-2</sup>, consistent with the non-detection of GJ 1151 in an earlier 3 ks Chandra pointing⁵⁷.

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<sup>+0.3</sup><sub>-0.1</sub> keV. The estimated 0.2 to 2.4 keV flux from such a model is (1.1±0.9) × 10<sup>-15</sup> erg s<sup>-1</sup> cm<sup>-2</sup>. 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<sup>+0.04</sup><sub>-0.05</sub> keV. The estimated 0.2 to 2.4 keV flux from such a model is (2.1±0.8) × 10<sup>-14</sup> erg s<sup>-1</sup> cm<sup>-2</sup>.

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⁻¹⁴ erg s⁻¹ cm⁻². This corresponds to an X-ray luminosity <2.3 × 10²⁶ erg s⁻¹, making LP 169-22 a quiescent candidate of our sample.

### 1.3. 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¹⁰. However, even with the LoTSS astrometric precision of ≈0.2" to 0.5", a crossmatch of total-intensity (Stokes I) and Galactic Gaia Data Release 2 (DR2) sources is dominated by chance co-incidence matches²³.

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 low-frequency-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%⁶¹,⁶²,². 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 $\pi/\pi_{err} \ge 3$. Since all the LoTSS pointings searched are at Galactic latitudes $|b| > 20°$, the average Gaia DR2 source density for sources with $\pi/\pi_{err} \ge 3$ is ≈ 1300 per sq. degree. Therefore, we would expect the number of chance-coincidence associations between a Gaia Galactic source and the Stokes V LoTSS sample to be ≈ 0.05. As performed by Callingham et al.²³, 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 $\pi/\pi_{err} \ge 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 cross-matching 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⁶³. 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²⁷. 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.

### 1.4. 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.²⁴ 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 $\bar{\omega}/\omega_{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.

## 1.5. 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 $\tau_{conv}$. The stellar rotation period can be measured directly from photometric observations⁵⁵,⁵⁶. 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.⁵⁷.

## 1.6. Circular polarisation sign convention

There are several instances in the collection, correlation, and reduction of LOFAR data where a sign-convention 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⁸⁸,⁸⁹. 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⁸⁹, we checked whether the Stokes V signs of pulsars observed in LoTSS were consistent with their circular polarisation handedness reported in the literature⁶¹,⁶²,⁸⁹. 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.⁶², Gould et al.⁶¹, and van Straten et al.⁸⁹: 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⁶¹,⁶². 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.

# 2. SUPPLEMENTARY INFORMATION

## 2.1. Literature information on detected stellar systems:

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_{\nu,rad}$, $L_X$, $L_{H\alpha}/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\alpha$ 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.

<table><thead><tr><th>Name</th><th>M<sub>*</sub><br>(M<sub>☉</sub>)</th><th>R<sub>*</sub><br>(R<sub>☉</sub>)</th><th>L<sub>ν,rad</sub><br>(×10<sup>14</sup> erg s<sup>-1</sup> Hz<sup>-1</sup>)</th><th>L<sub>X</sub><br>(×10<sup>28</sup> erg s<sup>-1</sup>)</th><th>L<sub>Hα</sub>/L<sub>bol</sub><br>(×10<sup>-4</sup>)</th><th>T<sub>B,plasma</sub><br>(×10<sup>12</sup> K)</th><th>Epoch<br>YYYY-MM-DD</th></tr></thead><tbody><tr><td>DO Cep</td><td>0.316</td><td>0.332</td><td>0.92±0.10</td><td>0.23</td><td>0.64</td><td>0.02</td><td>2019-02-07</td></tr><tr><td>WX UMa</td><td>0.095</td><td>0.121</td><td>0.45±0.09</td><td>0.36</td><td>1.21</td><td>0.04</td><td>2016-03-22</td></tr><tr><td>AD Leo*</td><td>0.420</td><td>0.431</td><td>0.96±0.10</td><td>3.20</td><td>1.72</td><td>0.15</td><td>2019-07-13</td></tr><tr><td>GJ 625*</td><td>0.317</td><td>0.332</td><td>0.80±0.09</td><td>0.04</td><td>0.06</td><td>0.004</td><td>2014-05-26</td></tr><tr><td>GJ 1151</td><td>0.167</td><td>0.190</td><td>0.63±0.15</td><td>0.02</td><td>0.07</td><td>0.005</td><td>2014-06-15</td></tr><tr><td>GJ 450</td><td>0.460</td><td>0.474</td><td>0.54±0.20</td><td>0.66</td><td>&lt;0.09</td><td>0.04</td><td>2017-09-21</td></tr><tr><td>LP 169-22</td><td>0.111</td><td>0.138</td><td>1.03±0.48</td><td>&lt;0.03</td><td></td><td>&lt;0.05</td><td>2014-06-08</td></tr><tr><td>CW UMa</td><td>0.306</td><td>0.322</td><td>4.23±0.44</td><td>5.37</td><td>1.82</td><td>0.22</td><td>2019-09-23</td></tr><tr><td>HAT 182-00605</td><td>0.442</td><td>0.422</td><td>2.94±0.57</td><td>3.40</td><td></td><td>0.13</td><td>2015-07-08</td></tr><tr><td>LP 212-62</td><td>0.161</td><td>0.183</td><td>28.34±1.55</td><td>0.38</td><td>0.84</td><td>0.03</td><td>2015-07-29</td></tr><tr><td>DG CVn</td><td>0.560</td><td>0.567</td><td>2.50±0.81</td><td>10.72</td><td>3.03</td><td>0.41</td><td>2018-04-17</td></tr><tr><td>GJ 3861</td><td>0.419</td><td>0.430</td><td>3.64±0.57</td><td>3.36</td><td>0.67</td><td>0.15</td><td>2017-04-19</td></tr><tr><td>CR Dra</td><td>0.823</td><td>0.827</td><td>43.38±2.46</td><td>36.65</td><td>1.38</td><td>0.91</td><td>2016-11-30</td></tr><tr><td>GJ 3729</td><td>0.472</td><td>0.481</td><td>9.49±0.80</td><td>7.54</td><td>1.64</td><td>0.32</td><td>2015-01-21</td></tr><tr><td>G 240-45</td><td>0.125</td><td>0.162</td><td>12.30±1.57</td><td>0.02</td><td></td><td>0.01</td><td>2018-05-02</td></tr><tr><td>2MASS J09481615+5114518</td><td>0.122</td><td>0.161</td><td>28.71±2.27</td><td>0.28</td><td></td><td>0.03</td><td>2018-08-18</td></tr><tr><td>LP 259-39</td><td>0.173</td><td>0.202</td><td>10.11±2.75</td><td>&lt;18.70</td><td></td><td>0.52</td><td>2018-06-11</td></tr><tr><td>2MASS J10534129+5253040</td><td>0.408</td><td>0.379</td><td>26.39±2.55</td><td>28.01</td><td>3.74</td><td>0.30</td><td>2015-03-14</td></tr><tr><td>2MASS J14333139+3417472</td><td>0.101</td><td>0.136</td><td>30.82±4.88</td><td>0.83</td><td>5.95</td><td>0.07</td><td>2015-10-10</td></tr></tbody></table>

## 2.2. 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 $\ge 50\%$ and brightness temperature $\ge 10^{12}$ K at 144 MHz. Therefore, incoherent gyrosynchrotron radiation can be confidently excluded¹. More detailed information about the emission physics can be found in Vedantham⁶⁴.

### 2.2.1. Plasma Emission

Plasma emission results from the conversion of Langmuir waves to transverse electromagnetic waves via non-linear processes¹. 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 \sim 10^8 K)$ is injected into an ambient colder plasma ($T_c \sim 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.⁶⁵ to calculate the brightness temperature $T_{B,\text{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⁸. 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⁶⁶,⁶⁷. We estimate the ambient coronal temperatures for our detected stars via $T_c = 0.11 F_X^{0.26}$, where $F_X$ is the X-ray surface flux calculated using a star's soft X-ray luminosity and the size of its photosphere⁶⁸. 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 \times 10^6$ K (GJ 625) to $1 \times 10^7$ K (LP 259-39), with a mean of $5 \times 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⁸ $L_n = 6 \times 10^9 (T_h/10^6 \text{ K})(R_*/R_\odot)^2(M_*/M_\odot)^{-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.⁶⁵ 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⁶⁹ and flares⁷⁰ on M dwarfs. In general, the plasma temperatures of M dwarf flares are less than 30 MK⁷⁰,⁶⁵,⁷¹.

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⁷². The estimated brightness temperature of fundamental plasma emission $T_{B,\text{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 $\approx 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, GJ 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%⁶⁶. For the most extreme configurations, where coronal loops occupy the entire stellar disk, the circularly polarised fraction of harmonic plasma emission is limited to $\approx 60\%$⁷³,⁷⁴,⁷⁵. 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⁷⁶ 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 ($\nu \sim 1$ GHz). Therefore, if disc-averaged $T_b \sim 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 \gg 10^{13}$ K at $\nu \gtrsim 1$ GHz. Such emission has not been observed even on the most active stars². 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³¹. 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.⁷⁷ detected coherent, narrowband emission from Proxima Centauri that persisted for over two-days at $\sim 1.4$ GHz. While the emission was 100% circularly polarised, the detected brightness temperature of $3 \times 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.⁶⁵ shows that the observed brightness temperature could readily be supplied by plasma emission. Additionally, Lynch et al.⁹ 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 $\gtrsim 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 $\sim 30$ min, while AD Leo and DO Cep have persistent emission lasting $\gtrsim 4$ hrs.

### 2.2.2. 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¹. In this case, the brightness temperature of a continuously operating maser is⁷⁸,⁸

$$T_{B,\text{flare}} \approx 2 \times 10^{14} (\lambda/200 \text{ cm})^2 (L/R_*)^{-1} (\beta/0.2)^2 \text{ K}, \quad (1)$$

where $\lambda$ is the wavelength of radiation, $R_*$ is the stellar radius, $L$ is the length-scale of the magnetic trap, and $\beta = 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¹.

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⁸. The instantaneous bandwidth of the loss cone maser is Δν/ν ≈ β²α² where α is the opening angle of the loss cone⁷⁸. 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⁷⁹. 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⁸⁰. The deposition of high-energy electrons at the poles of Jupiter set-up the population inversion necessary for ECMI to occur⁸¹. 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.³⁴, 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²⁶ and V374 Peg⁸². Additionally, the linearly-polarised ≈1 GHz bursts identified by Zic et al.³⁷ 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.³⁷ 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.⁴. Here we use the slightly modified expressions of Saur et al.⁸³ and Turnpenney et al.⁵ to compute the starward energy flux from a sub-Alfvénic interaction with an exoplanet:

$$S^{\text{th}} = \frac{\alpha^2}{2} \sin^2 \theta v_{\text{rel}}^2 B R_{\text{eff}}^2 \sqrt{4\pi\rho}, \quad (2)$$
in Gaussian units, where $\alpha$ is the interaction strength$^{83}$, $v_{rel}$ is the velocity of the stellar wind flow in the frame of the planet, $B$ and $\rho$ 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 $\epsilon_{rad}$. The flux density of the observed emission is then

$$F_{rad} = S\epsilon_{rad}/(\Omega d^2\Delta\nu), \qquad (3)$$

where $\Omega$ is the beam solid angle of the emitter, $d$ is the distance to the stellar system, and $\Delta\nu$ 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 $\Delta\nu = eB_*/(2\pi 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

$$F_{rad} = \frac{\pi m_e \alpha^2 \epsilon_{rad}}{e d^2 \Omega} R_{eff}^2 \sin^2 \theta v_{rel}^2 \sqrt{4\pi\rho} \frac{B}{B_*} \qquad (4)$$

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 $\epsilon_{rad} = 0.01$, a conservative estimate$^{34}$ as numerical simulations allow for $\epsilon_{rad} \approx 0.1$. We assume a beam-solid angle for the emission of $\Omega = 0.15$ as expected for an electron speed of $0.2c$. The starward Poynting flux inferred from observations can then be computed using equation 3 as

$$S^{\text{obs}} \approx 4 \times 10^{22} \left( \frac{F_{\text{rad}}}{\text{mJy}} \right) \left( \frac{d}{10 \text{ pc}} \right)^2 \left( \frac{\Omega}{0.15 \text{ sr}} \right) \left( \frac{B_*}{100 \text{ G}} \right) \text{ erg/s.} \qquad (5)$$

The theoretically expected Poynting flux can be computed from Equation 2. We adopt $v_{rel} = 500 \text{ km/s}$ which is typical of winds driven by a thermal expansion of $\sim 10^6 \text{ K coronae}$. We computed the density $\rho$ by conservatively adopting a base coronal density of $n_b = 10^7 \text{ 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 \text{ 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 $\theta$ 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 $\theta = 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 \text{ km}$. Finally, we assume $\alpha = 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

$$S^{\text{th}} \approx 6 \times 10^{22} \left( \frac{\sin \theta}{0.2} \right)^2 \left( \frac{v_{\text{rel}}}{500 \text{ km/s}} \right)^2 \left( \frac{B_*}{100 \text{ G}} \right) \left( \frac{R_{\text{eff}}}{6500 \text{ km}} \right)^2 \left( \frac{n_b}{10^7 \text{ cm}^{-3}} \right)^{1/2} \left( \frac{r}{10} \right)^{-3} \text{ erg/s.} \qquad (6)$$
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 are before to this study, and that of Vedantham et al.⁸, 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.

## 2.3. Güdel-Benz relationship

One test of whether our sample belongs to the class of canonical radio stars characterised by gigahertz-frequency observations¹⁰ is to compare 144 MHz radio $L_{\nu,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 \propto L_{\nu,rad}^{0.73}$, holds over 10 orders of magnitude in radio luminosity despite the fundamentally different emission mechanisms for the radio and soft X-ray radiation⁸⁴,⁸⁵. 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⁸⁵,³⁰.

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²⁶, 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⁻¹ Hz⁻¹ 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⁸⁶,⁸⁷.

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.
<figure><img src="image_5.png" /></figure>

**Supplementary Information Figure 2.** Soft (0.2-2.0 keV) X-ray luminosity $L_X$ against radio luminosity $L_{\nu,\text{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 \propto L_{\nu,\text{rad}}^{0.73}$). 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\alpha$ measurements (see Supplementary Information Figure 6). Non-detected X-ray sources are shown as 3-$\sigma$ 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.

## 2.4. *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><img src="image_6.png" /><figcaption>**Supplementary Information 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-\sigma$ uncertainty.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Supplementary Information Figure 4.** 144 MHz lightcurves for the detections listed in Table 1 (continued).</figcaption></figure>

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-\sigma$ uncertainty.
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-\sigma$, 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.

## 2.5. *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><img src="image_9.png" /></figure>

**Supplementary Information Figure 5.** 144 MHz detection rate as a function of rotation period for the M dwarfs studied by Newton et al.⁵⁵. While the sample produced by Newton et al.⁵⁵ 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 gigahertz-frequency studies²⁷. 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.⁵⁵ sample, the absolute detection percentages should not be extrapolated to other surveys. The error bars represent 1-σ uncertainty.

## 2.6. Hα characteristics of the population

Chromospheric Hα emission is a tracer of magnetic activity, and is known to have a strong relationship with rotation in M dwarfs⁵⁵. 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 Å.

## 2.7. 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⁹⁰. 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*².⁵ to be consistent with previously presented radio source counts⁹⁰.
<figure><img src="image_10.png" /><figcaption>**Supplementary Information Figure 6.** 144 MHz radio luminosity against $H\alpha$ luminosity by bolometric luminosity (left panel) and $H\alpha$ equivalent width (right panel). The ratio of $H\alpha$ luminosity to bolometric is relative to "inactive" stars defined by Newton et al.⁵⁵. Upper-limits and error bars represent $1-\sigma$ uncertainty.</figcaption></figure>

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⁹¹,⁹². 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.
<figure><img src="image_11.png" /><figcaption>**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.</figcaption></figure>

REFERENCES

39. Mechev, A. *et al*. An Automated Scalable Framework for Distributing Radio Astronomy Processing Across Clusters and Clouds. In *Proceedings of the International Symposium on Grids and Clouds (ISGC) 2017*, 2 (2017).

40. Taylor, M. B. TOPCAT and STIL: Starlink Table/VOTable Processing Software. In Shopbell, P., Britton, M. & Ebert, R. (eds.) *Astronomical Data Analysis Software and Systems XIV*, vol. 347 of Astronomical Society of the Pacific Conference Series, 29 (2005).

41. Pérez, F. & Granger, B. E. IPython: a system for interactive scientific computing. *Computing in Science and Engineering* 9, 21–29 (2007). URL http://ipython.org.

42. Jones, E., Oliphant, T., Peterson, P. & Others. SciPy: Open source scientific tools for python (2001). URL http://www.scipy.org/.

43. Hunter, J. D. Matplotlib: A 2d graphics environment. *Computing In Science & Engineering* 9, 90–95 (2007).

44. Astropy Collaboration *et al*. Astropy: A community Python package for astronomy. A&A **558**, A33 (2013).

45. Van Der Walt, S., Colbert, S. C. & Varoquaux, G. The numpy array: a structure for efficient numerical computation. *Computing in Science & Engineering* **13**, 22–30 (2011).

46. Tasse, C. *et al*. The LOFAR Two-meter Sky Survey: Deep Fields Data Release 1. I. Direction-dependent calibration and imaging. A&A **648**, A1 (2021).

47. Hancock, P. J., Murphy, T., Gaensler, B. M., Hopkins, A. & Curran, J. R. Compact continuum source finding for next generation radio surveys. MNRAS **422**, 1812–1824 (2012).

48. Hancock, P. J., Trott, C. M. & Hurley-Walker, N. Source Finding in the Era of the SKA (Precursors): Aegean 2.0. PASA **35**, e011 (2018).

49. van Diepen, G., Dijkema, T. J. & Offringa, A. DPPP: Default Pre-Processing Pipeline (2018). 1804.003.

50. Offringa, A. R. *et al*. WSCLEAN: an implementation of a fast, generic wide-field imager for radio astronomy. MNRAS **444**, 606–619 (2014).
51. Boller, T. et al. Second ROSAT all-sky survey (2RXS) source catalogue. A&A **588**, A103 (2016).

52. Fleming, T. A., Schmitt, J. H. M. M. & Giampapa, M. S. Correlations of Coronal X-Ray Emission with Activity, Mass, and Age of the Nearby K and M dwarfs. ApJ **450**, 401 (1995).

53. Schmitt, J. H. M. M., Fleming, T. A. & Giampapa, M. S. The X-Ray View of the Low-Mass Stars in the Solar Neighborhood. ApJ **450**, 392 (1995).

54. Evans, I. N. et al. The Chandra Source Catalog. ApJS **189**, 37–82 (2010).

55. Newton, E. R. et al. The Hα Emission of Nearby M Dwarfs and its Relation to Stellar Rotation. ApJ **834**, 85 (2017).

56. Díez Alonso, E. et al. CARMENES input catalogue of M dwarfs. IV. New rotation periods from photometric time series. A&A **621**, A126 (2019).

57. Wright, N. J., Newton, E. R., Williams, P. K. G., Drake, J. J. & Yadav, R. K. The stellar rotation-activity relationship in fully convective M dwarfs. MNRAS **479**, 2351–2360 (2018).

58. Callingham, J. R., Farrell, S. A., Gaensler, B. M., Lewis, G. F. & Middleton, M. J. The X-Ray Transient 2XMMi J003833.3+402133: A Candidate Magnetar at High Galactic Latitude. ApJ **757**, 169 (2012).

59. Callingham, J. R. et al. Anisotropic winds in a Wolf-Rayet binary identify a potential gamma-ray burst progenitor. Nature Astronomy **3**, 82–87 (2019).

60. Arnaud, K. A. XSPEC: The First Ten Years. In Jacoby, G. H. & Barnes, J. (eds.) *Astronomical Data Analysis Software and Systems V*, vol. 101 of Astronomical Society of the Pacific Conference Series, 17 (1996).

61. Gould, D. M. & Lyne, A. G. Multifrequency polarimetry of 300 radio pulsars. MNRAS **301**, 235–260 (1998).

62. Han, J. L., Manchester, R. N., Xu, R. X. & Qiao, G. J. Circular polarization in pulsar integrated profiles. MNRAS **300**, 373–387 (1998).

63. Henry, T. J. et al. The Solar Neighborhood. XVII. Parallax Results from the CTIOPI 0.9 m Program: 20 New Members of the RECONS 10 Parsec Sample. AJ **132**, 2360–2371 (2006).

64. Vedantham, H. K. On the mechanism of polarised metrewave stellar emission. MNRAS **500**, 3898–3907 (2021).

65. Stepanov, A. V., Kleim, B., Fürst, E., Jessner, A., Krüger, A., Hildebrand, J., & Schmittr, J. H. M. M. Microwave plasma emission of a flare on AD Leo. A&A **374**, 1072-1084 (2001).

66. Benz, A. O. Plasma astrophysics: Kinetic processes in solar and stellar coronae, vol. 184 (1993).

67. Reid, H. A. S. & Kontar, E. P. Langmuir wave electric fields induced by electron beams in the heliosphere. A&A **598**, A44 (2017).

68. Johnstone, C. P. & Güdel, M. The coronal temperatures of low-mass main-sequence stars. A&A **578**, A129 (2015).

69. Giampapa, M. S. et al. The Coronae of Low-Mass Dwarf Stars. ApJ **463**, 707 (1996).

70. Robrade, J. et al. Quiescent and flaring X-ray emission from the nearby M/T dwarf binary SCR 1845-6357. A&A **513**, A12 (2010).

71. Osten, R. A. et al. The Radio Spectrum of TVLM 513-46546: Constraints on the Coronal Properties of a Late M dwarf. ApJ **637**, 518 (2006).

72. Saint-Hilaire, P. et al. A Decade of Solar Type III Radio Bursts Observed by the Nançay Radioheliograph 1998-2008. ApJ **762**, 60 (2013).

73. Melrose, D. B., Dulk, G. A. & Smerd, S. F. The polarization of second harmonic plasma emission. A&A **66**, 315–324 (1978).

74. Melrose, D. B., Dulk, G. A. & Gary, D. E. Corrected formula for the polarization of second harmonic plasma emission. *Proceedings of the Astronomical Society of Australia* **4**, 50–53 (1980).

75. Ledenev, V. G. The directivity and polarization of radio emission at the second harmonic of plasma frequency from coronal magnetic loops. A&A **285**, 1019–1022 (1994).

76. Zaitsev, V. V. & Stepanov, A. V. The Plasma Radiation of Flare Kernelss. Solar Physics **88**, 297–313 (1983).

77. Slee, O. B., Willes, A. J. & Robinson, R. D. Long-duration Coherent Radio Emission from the dMe Star Proxima Centauri. PASA **20**, 257–262 (2003).

78. Melrose, D. B. & Dulk, G. A. Electron-cyclotron masers as the source of certain solar and stellar radio bursts. ApJ **259**, 844–858 (1982).

79. Benz, A. O., Conway, J. & Gudel, M. First VLBI images of a main-sequence star. A&A **331**, 596–600 (1998).
80. Hill, T. W. & Michel, F. C. Heavy ions from the Galilean satellites and the centrifugal distortion of the Jovian magnetosphere. *J. Geophys. Res.* **81**, 4561 (1976).

81. Cowley, S. & Bunce, E. Origin of the main auroral oval in Jupiter's coupled magnetosphere-ionosphere system. *Planetary and Space Science* **49**, 1067 – 1088 (2001). URL http://www.sciencedirect.com/science/article/pii/S0032063300001677. Magnetosphere of the Outer Planets Part II.

82. Llama, J., Jardine, M. M., Wood, K., Hallinan, G. & Morin, J. Simulating Radio Emission from Low-mass Stars. *ApJ* **854**, 7 (2018).

83. Saur, J., Grambusch, T., Duling, S., Neubauer, F. M. & Simon, S. Magnetic energy fluxes in sub-Alfvénic planet star and moon planet interactions. *A&A* **552**, A119 (2013).

84. Guedel, M. & Benz, A. O. X-Ray/Microwave Relation of Different Types of Active Stars. *ApJL* **405**, L63 (1993).

85. Benz, A. O. & Guedel, M. X-ray/microwave ratio of flares and coronae. *A&A* **285**, 621–630 (1994).

86. Berger, E. Flaring up All Over-Radio Activity in Rapidly Rotating Late M and L Dwarfs. *ApJ* **572**, 503-513 (2002).

87. Berger, E. et al. Simultaneous Multi-Wavelength Observations of Magnetic Activity in Ultracool Dwarfs. III. X-ray, Radio, and Hα Activity Trends in M and L dwarfss. *ApJ* **709**, 332-341 (2010).

88. Robishaw, T. & Heiles, C. The Measurement of Polarization in Radio Astronomy. *arXiv e-prints* arXiv:1806.07391 (2018).

89. van Straten, W., Manchester, R. N., Johnston, S. & Reynolds, J. E. PSRCHIVE and PSRFITS: Definition of the Stokes Parameters and Instrumental Basis Conventions. *PASA* **27**, 104–119 (2010).

90. Kellermann, K. I. & Wall, J. V. Radio Source Counts and Their Interpretation. In Hewitt, A., Burbidge, G. & Fang, L. Z. (eds.) *Observational Cosmology*, vol. 124 of *IAU Symposium*, 545 (1987).

91. Braun, R., Bonaldi, A., Bourke, T., Keane, E. F. & Wagg, J. Anticipated SKA1 Science Performance. SKA-TEL-SKO-0000818 (2017).

92. Pope, B. J. S., Withers, P., Callingham, J. R. & Vogt, M. F. Exoplanet transits with next-generation radio telescopes. *MNRAS* **484**, 648–658 (2019).