Radio emission as a stellar activity indicator

Radio observations of stars trace the plasma conditions and magnetic field properties of stellar magnetospheres and coronae. Depending on the plasma conditions at the emitter site, radio emission in the metre-and decimetre-wave bands is generated via different mechanisms, such as gyrosynchrotron, electron cyclotron maser instability, and plasma radiation processes. The ongoing LOFAR Two-metre Sky Survey (LoTSS) and VLA Sky Survey (VLASS) are currently the most sensitive wide-field radio sky surveys ever conducted. Because these surveys are untargeted, they provide an opportunity to study the statistical properties of the radio-emitting stellar population in an unbiased manner. Here we perform an untargeted search for stellar radio sources down to sub-mJy level using these radio surveys. We find that the population of radio-emitting stellar systems is mainly composed of two distinct categories: chromospherically active stellar (CAS) systems and M dwarfs. We also seek to identify signatures of a gradual transition within the M-dwarf population, from chromospheric or coronal acceleration close to the stellar surface similar to that observed on the Sun to magnetospheric acceleration occurring far from the stellar surface similar to that observed on Jupiter. We determine that radio detectability evolves with spectral type, and we identify a transition in radio detectability around spectral type M4, where stars become fully convective. Furthermore, we compare the radio detectability versus spectra type with X-ray and optical flare (observed by TESS) incidence statistics. We find that the radio efficiency of X-ray and optical flares, which is the fraction of flare energy channelled into radio-emitting charges, increases with spectral type. These results motivate us to conjecture that the emergence of large-scale magnetic fields in CAS systems and later M dwarfs leads to an increase in radio efficiency.


Introduction
Radio observations are excellent tracers of supra-thermal plasma in the coronae of stars, magnetic activity of ultracool dwarfs (UCDs; including brown dwarfs), and the mangetospheres of planets.This is because plasma oscillation and electric charges moving in a magnetic field emit at radio frequencies, implying a radio detection of a stellar system is sensitive to the magnetic fields and plasma dynamics of the source.Therefore, radio observations probe phenomena on these objects that are otherwise not easily accessible in other wavelengths (Bastian et al. 1998;Güdel 2002).For example, radio emission from exoplanets provides the only known technique to directly measure the planet's magnetic field strength and topology, which is an important aspect of exo-habitability (e.g.Driscoll & Barnes 2015;Airapetian et al. 2020;López-Morales et al. 2011;Khodachenko et al. 2007;See et al. 2014;Grießmeier et al. 2015).Additionally, plasma emission from stars can be used to constrain fundamental coronal plasma parameters such as plasma density and the radial structure and dynamics of coronal plasma (Toet et al. 2021), giving insights on the impact of stellar plasma on exoplanet atmospheres.
Depending on the plasma parameters at the emitter, radio emission in the metre-and decametre-wave bands can be generated via a variety of mechanisms: from incoherent ones such as free-free, gyromagnetic, and gyrosynchrotron processes (Nin-⋆ Email: yiu@astron.nldos 2020) to coherent ones such as plasma radiation (Dulk 1985) and electron cyclotron maser instability (ECMI; Wu & Lee 1979;Melrose & Dulk 1982, 1984;Treumann 2006).For the Sun, there exist two principal radio components: the unpolarised or weakly polarised broadband solar flares at gigahertz (GHz) frequencies, and the narrowband coherent solar bursts at megahertz (MHz) frequencies.The unpolarised radio emission from solar flares have generally been attributed to gyrosynchrotron radiation (Alissandrakis 1986), produced by mildly relativistic electrons (γ ≲ 2−3) either arising from the Maxwellian tail of a thermal distribution (Dulk 1985), or from a non-thermal energy distribution produced via flare acceleration and/or magnetic reconnection (Parker 1988;Sweet 1958a,b;Matthews 2019).Many stellar radio flares on isolated main-sequence stars (e.g.Bastian 1990;Güdel & Benz 1993) as well as very close binaries (e.g.Drake et al. 1989;Mutel et al. 1985) are often observed to be broadband and unpolarised.On the contrary, coherent solar bursts with narrow instantaneous bandwidths are much more luminous and more circularly polarised compared to gyrosynchrotron flares.These bursts are typically produced by plasma radiation processes (Dulk 1985;Melrose 1980Melrose , 2009Melrose , 2017;;Güdel 2002), thus emitting at plasma frequencies or its harmonics that are typically below the GHz regime.Regardless of the mechanism, all solar radio-emitting charges are expected to be created via chromospheric/coronal acceleration close to the stellar surface due to magnetic reconnection, driven by immense shear in the solar radiative-convective interface (known as tachocline; Spiegel & Zahn 1992).We refer to this broad class of acceleration mechanisms collectively as the "Sun-like radio engine".
In comparison, radio observations of Jupiter show that Jupiter can generate luminous auroral radio bursts with strong circular polarisation (reaching 100%; Zarka 1998Zarka , 2004Zarka , 2007) ) at the ambient cyclotron frequency and its harmonics through the ECMI mechanism.Therefore, Jupiter's plasma dynamics are fundamentally different to the Sun due to its fast rotation rate, strong magnetic field, and low-density magnetosphere.Jovian radio emission is thus magneto-rotation-driven instead, with the magnetospheric acceleration occurring far from the planetary surface as a result of a breakdown of co-rotation between the magnetic field and plasma.An alternate acceleration mechanism known as centrifugal breakout, which also operates far from the surface, has been suggested for radio emission from hot magnetic stars with dipole-dominated magnetopspheres (Shultz et al. 2020;Owocki et al. 2020Owocki et al. , 2022)).We shall use the umbrella term "Jupiter-like engine" for acceleration in a magnetosphere far from the surface of the central body.
In contrast to the Sun, the coherent radio emission of UCDs are generally interpreted as electron cyclotron maser (ECM) emission powered by Jupiter-like engines (e.g.Hallinan et al. 2007Hallinan et al. , 2008;;Kao et al. 2016;Turnpenney et al. 2017), since their highly circular-polarised radio periodic bursts share many properties with that of Jupiter and other magnetised solar-system planets (Zarka 2007).The change in radio emission mechanism between solar-like stars and UCDs suggests that there exists a transition in radio-emitting engine of stellar systems -from the Sun-like paradigm to the Jupiter-like paradigm -somewhere in the M-dwarf sequence.Finding this transition is vital in understanding the radio detectability of stars and the predominant emission mechanism of radio-bright stellar systems.It may also shed light on the effect of the transition from partially convective to fully convective interiors (i.e. the tachocline no longer exists) on stellar activity (around spectral type M4; Chabrier & Baraffe 2000;Dorman et al. 1989;Stassun et al. 2011;Williams et al. 2014).
There have been previous studies that investigate and characterise the radio properties of stellar systems, from the hottest stars to the coldest brown dwarfs, in order to determine trends (see e.g.Bieging et al. 1989;Güdel & Benz 1993;Williams et al. 2014;Kao et al. 2016;Villadsen & Hallinan 2019).However, all of the above studies employed some target selection strategy, and thus could be implicitly biased towards special cases of stars and UCDs that may not reflect the general stellar population as a whole.One benefit of searching for radio-bright stellar systems in wide-field radio surveys is therefore to bypass these biases.However, this method is not without some major challenges.Helfand et al. (1999) performed the first large unbiased study of radio stars using the Faint Images of the Radio Sky at Twenty-Centimeters (FIRST; Becker et al. 1995), and concluded that FIRST astrometry alone was insufficient to avoid crippling chance-coincidence rates, especially when good proper-motion information was unavailable.Kimball et al. (2009) conducted a search for radio stars by combining FIRST with optical data from the Sloan Digital Sky Survey (SDSS), but they estimated that 108±13 out of 112 candidate radio stars were contamination, i.e. optical-faint radio quasars in chance alignment with a foreground star.Thus, they concluded a radio survey with a much higher sensitivity and resolution compared to FIRST was necessary to confidently identify radio stars.Antonova et al. (2013) conducted a 4.9-GHz volume-limited radio survey of 32 nearby UCDs with spectral types M7-T8, but failed to detect any ra-dio emission from them.They thus highlighted that low rotation rates and long-term variability of these targets might be the cause of the non-detections.Lenc et al. (2018) conducted the first lowfrequency circular-polarised all-sky survey using the Murchison Widefield Array (MWA; Tingay et al. 2013).However, with its relatively low sensitivity of ≈ 3mJy beam −1 , the survey only detected pulsars which were previously known to be radio-bright.
Studies with a newer crop of radio surveys have been more successful.Pritchard et al. (2021) presented results from a circular polarisation survey for radio stars in the Rapid Australian Square Kilometre Array Pathfinder (ASKAP) Continuum Survey (RACS;McConnell et al. 2020), and they identified M dwarfs, close binaries, young stellar objects, and chemically peculiar A-and B-type stars in the radio sample.Callingham et al. (2021b) detected coherent radio emission of M dwarfs from an untargeted flux-limited low-frequency survey using LOw-Frequency ARray (LOFAR; van Haarlem et al. 2013).They highlighted the benefit of utilising circular polarisation (Stokes V) information to suppress false associations (Callingham et al. 2019).More recently, Driessen et al. (2023) identified radio stellar sources from multiple radio surveys using their proper motions provided by Gaia Data Release 3 (DR3; Gaia Collaboration et al. 2022), and showed that the transient nature of radio stellar sources makes surveys with multiple epochs more preferable for searching them.
Having a wide-field radio survey with high sensitivity and angular resolution is therefore paramount in order to identify radio-bright stellar systems from the radio data, lest nondetections and extragalactic contamination dominate the sample.The ongoing LOFAR Two-metre Sky Survey (LoTSS; Shimwell et al. 2017) and the Karl G. Jansky Very Large Array (VLA; Perley et al. 2011) Sky Survey (VLASS; Lacy et al. 2020) are some of the most sensitive wide-field radio sky surveys ever conducted.Because the surveys are untargeted, they provide an opportunity to study the statistical properties of the radio-emitting stellar population in an unbiased manner.Moreover, Callingham et al. (2023) recently published V-LoTSS: the circular-polarised component of LoTSS.This provides yet another advantage as the Stokes V radio sky is orders of magnitudes sparser than its Stokes I counterpart (Callingham et al. 2023).In this paper, we present our latest efforts to identify and study the radio emission from stellar systems in these surveys as a way to understand if and how radio activity evolves with spectral type, and see whether a transition in acceleration mechanism (the so-called radio-emitting engine) exists, specifically from the aforementioned Sun-like to Jupiter-like engine as one goes from earlier to later spectral types.Additionally, to understand how stellar activity impacts radio luminosity and thus detectability, we also compare the radio detection rate in our sample to canonical activity indicators such as optical and X-ray flares.
The paper is structured as follows: in Sect.2, we present details of the radio survey catalogues, Gaia Catalogue of Nearby Stars (GCNS) and the TESS/X-ray flare statistics used in our analysis.Sect. 3 contains the description of our methodology, including the radio × GCNS crossmatching and TESS-flare-rate debiasing procedure.We present our results in Sect. 4 and conclude in Sect. 5.This paper contains many acronyms, and thus we include a table of acronyms (Table D.1) in the Appendix for clarity.Notes.The symbols ν survey , AA, t exp , and Ω survey represent radio survey frequency coverage, absolute astrometric accuracy, exposure time, and the sky coverage respectively.Note that the final uncertainty in the radio position must also factor the astrometric precision that depends on the signal-to-noise ratio (SNR) of the detection.For example, LoTSS source with <1 mJy has an average astrometric precision of 0.5 ′′ (Callingham et al. 2019).The Ω survey of V-LoTSS is slightly larger since Callingham et al. (2023) searched for sources further down the primary beam of single LoTSS pointing, where the noise would have been higher than if the data was mosaicked.

Radio sky surveys
In order to draw conclusions regarding the radio-bright stellar population as a whole, we begin by crossmatching different radio catalogues to optical positions of known stars in the solar neighbourhood.The following section details the radio sky surveys used in this work.A summary of the important parameters of these surveys can be found in Table 1.

LOFAR Two-metre Sky Survey
The LOFAR Two-metre Sky Survey (LoTSS; Shimwell et al. 2017) is an ongoing low-frequency (120-168 MHz) radio survey of the Northern Sky.Each LoTSS pointing is observed for 8 hours, reaching a median 1σ rms noise of around 83 µJy/beam.
Here we use the LoTSS radio catalogue from the second data release (DR2; Shimwell et al. 2022) covering 27% of the Northern Sky (5634 square degrees) containing around 4 million Stokes I sources.The unprecedented depth of LoTSS allows us to reveal a low-frequency radio stellar and substellar population neverbefore-seen as previous searches below ∼300 MHz generally lacked the required sub-mJy sensitivity necessary to detect the general population (Callingham et al. 2021b;Vedantham et al. 2020Vedantham et al. , 2023)).However, using the LoTSS Stokes I catalogue alone to identify radio stellar system candidates is not straightforward due to the potentially crippling rate of chance-coincidence associations (i.e. the so called false-positive matches; more details in Sect.3.1.1).To circumvent this, we also include the available circular polarisation (Stokes V) information of LOFAR-detected stellar sources, described in the following section.

Circular-polarised sky of LoTSS
The source density of the Stokes V radio sky is > 5 orders of magnitude lower than the Stokes I radio sky because most radio sources are extragalactic objects powered by the synchrotron mechanism (Begelman et al. 1984), and therefore do not have a significant degree of circular polarisation (Rayner et al. 2000;Beckert & Falcke 2002).As anticipated, only 1 extragalactic source (an active galactic nucleus) has a detectable circular polarisation (≈ 1%) in LoTSS-DR2 (Callingham et al. 2023).In rare cases where extragalactic radio sources are circularly polarised, the known radio-emission mechanisms of these objects would only allow them to have at most ≈ 1% circular polarised fraction (Saikia & Salter 1988;Valtaoja 1984;Wardle & Homan 2003), which matches with observations (e.g.Macquart et al. 2003;Agudo & Thum 2022).
On the other hand, the circular polarisation of some stellar radio sources are known to be able to reach ≈ 100% (e.g.Hallinan et al. 2006Hallinan et al. , 2007Hallinan et al. , 2008;;Williams & Berger 2015) owing to plasma radiation processes and ECMI (Dulk 1985;Vedantham 2021).Callingham et al. (2023) presented V-LoTSS, which consists of Stokes V maps of LoTSS-DR2 with median 140 µJybeam −1 and a resolution of 20 ′′ .As we expect the radio emission of stellar systems to be time variable, an overall mosaic of the fields may wash out genuine detections as these sources may not be emitting in Stokes V in an adjoining field that was observed on a different date.Therefore, using V-LoTSS is again advantageous since unlike LoTSS-DR2, V-LoTSS performs the Stokes V search on individual LoTSS pointings rather than on a mosaicked image.

VLA Sky Survey
As the aim of this paper is to determine how much radio emission is indicative of stellar activity, we are also interested in the radio population that has a higher frequency than the LOFAR band, as different radio frequencies probe different radio emission mechanisms (Güdel 2002).For example, a flaring star can emit intense gyrosynchrotron radiation which typically falls under the decimetre-wave regime (Nindos 2020), a frequency range observable by the VLA.However, such a star might not have detectable emission in the LOFAR band due to self-absorption and brightness temperature limitations: the brightness temperature of gyrosynchrotron emissions cannot exceed inverse Compton limit of T B ∼ 10 12 K (Kellermann & Pauliny-Toth 1969).Since T B ∝ ν −2 , it can be trivially shown that gyrosynchrotron cannot produce emission in the LOFAR band detectable with the LoTSS sensitivity.Conversely, stellar systems detected in the metrewave regime are much more likely to be due to plasma radiation processes and/or ECMI (Vedantham 2021).These stars may have a large-scale magnetic field of few hundred gauss which can produce strong coherent ECM radiation that peaks at LO-FAR band but decays sharply at VLA band.Therefore, we also utilise the VLA Sky Survey (VLASS; Lacy et al. 2020) to select radio-bright stellar systems, and see how different the VLA-detected stellar population is when compared to the LOFAR-detected population.VLASS is an ongoing multi-epoch S-band (2)(3)(4) continuum radio survey covering the whole sky visible to the VLA (i.e.δ > −40 • ; 33885 square degrees) and aims to produce images with high angular resolution (≈ 2.5 ′′ ) and 1σ rms noise of ≈ 130µJy.Three epochs of observation (with two sub-epochs each) are planned and currently the first two epochs (VLASS1 & VLASS2) have been fully observed and processed.We use the latest release of VLASS which contains Epochs 1 and 2 'Quick Look' (QL) Catalogues1 , each of them containing around 3 million sources.We note that currently the VLASS data releases do not include Stokes V information.

Gaia Catalogue of Nearby Stars
To identify radio-bright stellar systems, we crossmatch the aforementioned radio catalogues with Gaia, the largest available star catalogue with detailed astrometric information.Previously, Callingham et al. (2021b) examined the LoTSS Stokes V maps for ≥ 4σ sources and crossmatched those detections to sources in the Gaia Data Release 2 (DR2) to search for radio counterparts to stellar systems.However, the Gaia DR2 catalogue becomes significantly incomplete for spectral types later than ∼M7 (Kiman et al. 2019).Therefore, in this work, we instead use the Gaia Catalogue of Nearby Stars (GCNS; Gaia Collaboration et al. 2021b).The GCNS is a catalogue of well-characterised objects within 100 pc of the Sun from the Gaia Early Data Release 3 (EDR3; Gaia Collaboration et al. 2021a).The GCNS contains 331312 objects, 40234 of which are within 50 pc.The fact that the catalogue is volume-complete for all objects earlier than M8 at the nominal G = 20.7 magnitude limit2 of Gaia (Gaia Collaboration et al. 2021b) is essential as our aim is to see whether radio detection rates of stellar systems evolves with spectral type across the stellar main sequence.

Flare statistics
As we aim to understand the correlation between the radio emission of stars and stellar activity, statistics of stellar flares occurrence are also needed in order to analyse if a higher flaring activity implies a higher chance of radio detection.The motivation for searching such a correlation is that flares may be the essential in accelerating the radio-emitting electrons.To do this, we use the study by Günther et al. (2020) for optical stellar flares from the stars observed by NASA's Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015), and by Johnstone et al. (2021) for Xray stellar flares statistics using the NEXXUS database (Schmitt & Liefke 2004).

TESS flares
TESS is an all-sky survey mission designed to discover exoplanets orbiting bright nearby main-sequence stars using transit photometry.As TESS has observed numerous stars (continuously for days, some even months) and measured their light curve in order to find exoplanet transits, TESS data also contains valuable information on the incidence of stellar flares.Günther et al. (2020) performed a study of optical stellar flares for the 24809 stars observed with 2-minute cadence during the first two months of the TESS mission.Most importantly, they studied the flare rate and energy as a function of stellar type and rotation period.
As seen in Figs. 4 and 5 reported by Günther et al. (2020), flares are most commonly detected on M dwarfs, especially on mid to late spectral type (beyond M4), where more than 40% of the stars showing observable flares lie.Stars of spectral type earlier than M-type have significantly lower observable flare occurrences of less than 10%.This general picture is also consistent with past findings from e.g.Kepler (Davenport 2016;Van Doorsselaere et al. 2017) and MEarth (Mondrik et al. 2019) cat-alogues of stellar flares: stars of later spectral types (especially fast-rotating, young M dwarfs) are the most likely to flare, and that their flare amplitude is independent of the rotation period (cf.Maehara et al. 2012).
Although Günther et al. (2020) only present findings derived from the first two months (i.e.Sectors 1 and 2) of TESS data, here we use the data from the first 2 years (Sectors 1-26) of TESS mission instead (Günther et al., in prep.).Such a database increases the number of flaring stars detected by a factor of ≈ 13, greatly improving the flare statistics.The data for the TESS short-cadence targets from Sectors 1-26 is publicly available3 .
In addition, for our analysis in Sect.4.4.1, we now also account for the average TESS flare rate ν TESS (i.e.number of stellar flares per unit time) in each spectral-type bin, instead of merely considering the fraction of flaring stars provided by Günther et al. (2020).The motivation is to compare optical flare rate to the radio detectability as a function of spectral type.If stellar activity has direct correlation to radio detectability, then it would be expected that the more frequently a star flares, the more likely it is detectable in radio within any given observation window.We describe the definition of a flaring star and the method of how TESS flares are counted in Sect.3.2.
The additional data obtained from year 1 and 2 further reinforce the original conclusion by Günther et al. (2020), where they found M dwarfs of type M4-M6 dominate the TESS sample of flaring stars than any other star (Günther et al., in prep.).Note that, however, Günther et al. (2020) pointed out the several biases in their TESS flare study as a consequence of TESS only being able to observe in relative flux units.And so, a common flare energy threshold for all stars is thus necessary for an unbiased comparison of flare rates between different spectral types of stars.The several TESS flare biases and subsequently our method for debiasing the TESS flare statistics are described in Sect.3.3.

X-ray flares
Besides optical stellar flares like the ones observed by TESS, Xray flares are also indicative of stellar activity.The X-ray luminosity of a star is empirically determined by its mass, age, and rotation rate (Pizzolato et al. 2003).This empirical correlation implies stellar activity evolution is linked to its rotational evolution, since many activity parameters, such as magnetic dipole field strength and mass loss rate, are tightly correlated with the Rossby number R o = P rot /τ c , where P rot is the rotation period and τ c is the convective turnover time (Wright et al. 2011).
Therefore, we also utilise the model by Johnstone et al. (2021) for X-ray (0.1-2.4 keV) flare rates of F, G, K, and M dwarfs, with masses between 0.1 and 1.2 M ⊙ .As seen in Fig. 19 in their work, the rates of X-ray flares above a fixed energy threshold monotonically decline with declining stellar effective surface temperature.Such a conclusion is the exact opposite of that from TESS flare study by Günther et al. (2020), in which lower mass stars seemingly flare more than higher mass ones.One reason behind such discrepancy may stem from the different definition of flares in these two studies.Johnstone et al. (2021) only includes flares of energy above an energy threshold of 10 32 erg in their analysis, unlike the TESS flare statistics where all impulsive changes in relative flux are classified as flares based on an inference framework called allesfitter (Günther & Daylan 2021) and complementary criteria such as high signalto-noise ratio (SNR).Therefore, although debiasing is required to compare TESS stellar flare rates of different spectral types (more details in Sect.3.3), there is no such need for any debiasing in the X-ray flare statistics owing to the fact that Johnstone et al. (2021) only counts flares above an energy threshold.

Crossmatching method
Now that we have introduced the different all-sky surveys used in this study, we outline the method we applied for crossmatching these catalogues.

False alarm rate
As mentioned in Sects. 1 and 2.1.1,one key challenge of using the Stokes I catalogue with Gaia alone to identify radio stellar systems is that the radio catalogue is composed mostly of galaxies.With a high density of radio sources (≈ 780 sources per square degree on average), the LoTSS sky is far denser than any previous wide-area radio survey such as the NRAO VLA Sky Survey (NVSS; Condon et al. 1998), FIRST (Becker et al. 1995), TIFR GMRT Sky Survey first alternative data release (TGSS ADR1; Intema et al. 2017), and GaLactic and Extragalactic Allsky MWA (GLEAM; Hurley-Walker et al. 2017) survey.Therefore, despite the LoTSS's high astrometric precision of 0.2 ′′ to 0.5 ′′ (Shimwell et al. 2022), true associations between dense optical surveys and radio sources cannot be confidently done by simple blind crossmatching; Callingham et al. (2019) showed that a blind search for radio-bright stellar systems in Gaia and LoTSS is dominated by false positives, and thus either additional observational or physically motivated information is needed in order to form a reliable sample of Galactic Gaia-LoTSS counterparts.For example, assuming the radio sources in the sky surveys are homogeneous, the following equation gives us an approximate estimate of false matches with the Gaia catalogue: where N false is the number of chance-coincidence associations, n radio and Ω radio are the source density and sky coverage of the radio survey respectively, and θ is the crossmatching radius.For LoTSS-Gaia DR3 crossmatching, this gives us N false ∼ 10 5 , hence highlighting the difficulty of confidently associating stellar systems with radio emission from LoTSS alone.Similarly, the VLASS catalogue also faces the same issue of false positive as it has a source density only a factor of ≈ 8 lower than that of LoTSS.

LoTSS × GCNS match
Nevertheless, identifying radio-bright stellar systems in the LoTSS Stokes I catalogue while keeping chance-coincidence associations small is possible.We achieve this by crossmatching with GCNS within 50 pc and by setting a small crossmatching radius of 1.4 ′′ with proper motion correction based on Gaia information.This crossmatching radius corresponds to ≈ 1 chancecoincidence association on average.

V-LoTSS × GCNS match
On the other hand, as mentioned in Sect.2.1.2,the circularpolarised sky of LoTSS allows us to be significantly more confident in the radio sources' association with stellar systems.Since the number of detected sources in V-LoTSS is less than a hundred, we need not restrict the crossmatching with GCNS to be within 50 pc.Moreover, the low source density also allows us to set the crossmatching radius generously to be 6 ′′ (which corresponds to a false association rate of 10

Average TESS ensemble flare rates
We use the catalogue for all individual flares found in TESS year 1 and 2 (Günther et al. 2020 and in prep.).Each entry contains information on a candidate flare from a TESS short-cadencetargeted star, such as the flare's peak time, full width at half maximum (FWHM), and flare amplitude (see Günther et al. 2020).
The new catalogue also contains extra information on qualitycontrol filters and a probability of the candidate being a true flare, which is computed by the convolutional neural network stella (Feinstein et al. 2020a,b).This extra information is also detailed in Feinstein et al. (2022).In this work, we follow the definition of TESS flare rate per star ν TESS by Feinstein et al. (2022): where N is the number of flares for a given star, p i is the stella probability for each flare candidate from that star, and t obs is the total observed time of the star by TESS.By binning the flare rates into effective temperature T eff bins, we obtain a relationship between stellar spectral type and the average ν TESS .This curve, when multiplied by the fraction of flaring stars in each T eff bin, gives us the average likelihood of TESS detecting a flare on a star of a particular spectral type (see Fig. 1a).We shall refer to this quantity, i.e.N flaring /N total × ν TESS , as the "average TESS ensemble flare rate".Intuitively, this quantity tells us that if one were to pick a star at random and observe it with TESS, how many flares on average would be seen on that star given its particular spectral type.

Debiasing TESS flare statistics
As mentioned in Sect.2.3.1, the TESS flare statistics used here are inherently biased since they defined "flaring" as a sufficient increase in flux amplitude compared to its quiescent level (i.e.sufficient increase in relative flux units), rather than a lower energy threshold above which all flares are counted.These TESS flare statistics biases are described in depth in the following sections.

TESS sampling bias
Firstly, TESS' short-cadence observations are biased towards brighter stars to ensure follow-up spectroscopic observations of transiting exoplanets (Collins et al. 2018).Indeed, TESS is essentially a targeted survey with a complicated target selection (for light curves) determined by the myriad of short-cadence proposals submitted by the scientific community.We can minimise this TESS sampling bias by excluding stars beyond a cer-tain distance such that TESS has observed every star within such volume; around 15 pc is where the TESS sample is volumecomplete (Günther et al. 2020 and in prep.).However, the number of flaring stars with early spectral type within 15 pc is too small to perform proper statistics.In particular, there is almost no star earlier than K0 within 15 pc in the TESS data.Therefore, we compromise on considering TESS stars within 50 pc instead, as this achieves a good balance between sample completeness and a sufficient sample size for proper statistical analysis.As shown in Fig. 1a, the orange curve -representing the average-TESSensemble-flare-rate that excludes TESS stars beyond 50 pc -has a very different shape when compared to the green curve which does not exclude any star.In particular, hotter stars with spectral type earlier than ∼K5 have significant amount of flaring events in the green curve compared to the orange curve.In addition, along with the usual mid-M peak, there is an additional peak around K0 which is absent in the orange curve.However, M dwarfs (especially around M4) still show the highest average TESS ensemble flare rate (N flaring /N total × ν TESS ) in both curves.The deviations between the curves are the manifestation of the TESS sampling bias from which the full TESS flare dataset (i.e.green curve) suffers.For example, Fig. 1a suggests that flaring solar-type stars are oversampled by TESS.One of the possible reason for this oversampling is that the scientific community is interested in studying the superflares of solar-type stars in order to obtain insights on space weather around solar analogs (e.g.Namekata et al. 2017;Maehara et al. 2017;Tu et al. 2020).In any case, we shall only consider TESS stars within 50 pc in subsequent analysis, as such a 50-pc sample achieves a good balance between minimising the TESS sampling bias and having enough flaring stars to avoid small-number statistics.

TESS flare detection biases
Secondly, there exists an interplay of two opposing flaredetection biases: (i) since TESS detects flares in relative flux units, a flare of a given energy is more readily detected on cooler stars because of a larger contrast between the flare and the stellar quiescent flux density; and (ii) a higher photometric noise in cooler stars decreases the SNR of their relative flare flux densities.As shown in Fig. 2f of Günther et al. (2020), the recovery rates of flares in injection tests is significantly lower for cooler stars in the M-dwarf range compared to F/G/K dwarfs.And so, in order to have a fair flare-rate comparison between stars of different spectral types, these biases must be taken into account.In the following paragraphs, we describe the procedure for the debiasing of the TESS flare statistics.Here, we start by converting the known flare completeness in relative flux units (Günther et al. 2020; in Fig. 2f), to a flare completeness in absolute flux units: Let A S be the relative flux of a stellar flare observed by TESS with respect to the corresponding quiescent stellar flux.We define A min = L ′ f,min /L ′ * as the minimum relative flux of which a stellar flare of luminosity L f,min can still be detectable from a star with an effective temperature T eff and a stellar quiescent luminosity L * , i.e.A min = A min (T eff ), with the prime symbol ′ denoting the luminosity in the TESS observing bandpass.This function encapsulates the interplay of the two aforementioned detection biases of TESS flares.The value of A min (T eff ) is already known from Fig. 2f by Günther et al. (2020).The stellar quiescent luminosity L ′ * can be expressed as where R * is the stellar radius, B * is the Planck function at the T eff of the star, and S TESS is the TESS spectral response function which is defined as the product of the long-pass filter transmission curve and the detector quantum efficiency curve (Ricker et al. 2015).To relate the stellar radius R * to its T eff , we use the Basti-IAC isochrone model (Hidalgo et al. 2018) with a solartype metallicity for stars with T eff ≳ 3000 K and the BHAC15 isochrone model (Baraffe et al. 2015) for main sequence lowmass stars down to T eff = 2600 K.We assume a stellar age of 1 Gyr (around the same order of magnitude of most stars in our galaxy, e.g.Haywood et al. 2013;Snaith et al. 2015) for these models.Fig. 1b shows R * as a function of T eff and the corresponding smooth spline interpolation.
We can determine the minimum flare luminosity observable by TESS: Now, to relate this quantity to its unprimed counterpart, we assume that flares have a blackbody spectrum with the same effective temperature no matter the spectral type of the star (e.g.Shibayama et al. 2013;Davenport 2016;Yang et al. 2018;Günther et al. 2020;Gao et al. 2022).The precise value of effective temperature does not matter in subsequent analysis 4 as we only are only interested in the ratio between the flare luminosities from stars of different spectral types, not their absolute values.This assumption implies that for two stars with stellar effective temperatures T 1 and T 2 respectively, we have . Therefore, we now have L f,min as a function of minimum relative flux A min and T eff .Moreover, we assume the average duration of a stellar flare is the same for stars of all spectral type, thus , where E f (T ) is the flare energy from a star with T eff = T .Since the curve of A min vs T eff is already given by Günther et al. (2020), we can determine the minimum flare energy E f,min detectable by TESS on a particular star with T eff = T .As seen in Fig. 1c, the weakest TESS-detected flare on a early-type star is much stronger than its late-type star counterpart.For example, the E f,min of a solar-like G star is around ∼2 orders of magnitude stronger than the E f,min of an M dwarf.
Lastly, we need to convert this flare energy threshold into a flare rate ν.To do this, we utilise the study by Gao et al. (2022) regarding cumulative flare frequency distributions (FFDs; e.g.Lacy et al. 1976;Hawley et al. 2014;Günther et al. 2020;Jackman et al. 2021) detected by TESS and Kepler.The cumulative FFD represents the number of flares in unit time with an energy greater than a particular value E, i.e. how often a flare of energy E f ≥ E occurs on a star.Mathematically, where α cum < 0 and β are empirically determined parameters.The debiased flare rate for a star with spectral type S is then log 10 (ν S /ν TESS,S ) = α cum log 10 (E 0 /E f,min,S ), (6) where E 0 is the common energy threshold above which all flares are counted for all stars, ν TESS,S and E f,min,S are the TESSobserved flare rate and the minimum TESS-detected flare energy on star of type S respectively.Gao et al. (2022, in Fig. 12) empirically determined that α cum from solar-type stars to mid M dwarfs is approximately consistent with a value of −1, and so we shall assume α cum = −1 in subsequent calculations.We can now convert the ratio of observed flare rates for stars of types S 1 and S 2 into a ratio of their true flare rates above some energy threshold as log 10 (ν S 1 /ν TESS,S 1 ) − log 10 (ν S 2 /ν TESS,S 2 ) = − log 10 (E 0 /E f,min,S 1 ) + log 10 (E 0 /E f,min,S 2 ).( 7) Hence, Finally, multiplying the two blue curves in Figs.1a and 1c gives us Fig. 1d, which shows the the debiased average TESS ensemble flare rate, i.e. the fraction of flaring stars N flaring /N total weighted by the debiased flare rate ν.Here, given a particular energy threshold above which flares are counted, we can see a general trend of the ensemble flare rate decreasing as we move from early-type stars towards K dwarfs.Then, at T eff ≈ 4800 K (around spectral type K3), the ensemble flare rate increases, peaking at spectral type M0 before decreasing again in the range of mid-M dwarfs.Specifically, early-type (≲G0) stars are actually around 2 − 4 times more likely to flare than early-K and mid-M dwarfs, while the ensemble flare rate around spectral type M0 is quite comparable.This debiased TESS flare trend is not exactly in line with the X-ray flare rate statistics by Johnstone et al. (2021).Regardless, we have now obtained a debiased average TESS ensemble flare rate curve, which shall be used in the TESS flare statistics analysis in Sect.4.4.1.

Crossmatching Results
In Tables A.1, A.2, and A.3 (all shown in Appendix A), we present the radio × GCNS crossmatching samples from the three aforementioned radio surveys.There are 22 V-LoTSS-detected sources that have optical counterparts in the Gaia Catalogue for Nearby Star (GCNS), while there are 25 LoTSS-detected sources and 65 VLASS-detected sources that have an optical counterpart in GCNS within 50 pc.Figures 2 and 3 show each of three radiodetected stellar populations in the Hertzsprung-Russell (HR) diagram.
For the V-LoTSS × GCNS sample, note that out of the 22 matches, all except two have an angular separation between the V-LoTSS source (based on Stokes I astrometry) and Gaia counterpart of ≲ 1 ′′ .The remaining two sources, HAT 182-00605 and II Pegasi, have an angular separation of 1.22 ′′ and 1.83 ′′ respectively.However, since both of them are already known to be radio-bright (Callingham et al. 2021b;Toet et al. 2021), we are confident that none of the matches in our V-LoTSS × GCNS sample are false positives.Moreover, all of these matches are consistent with the stellar sample in V-LoTSS, as Callingham et al. (2023) also crossmatched V-LoTSS sources to the Gaia Data Release 2 & 3 (DR2 & DR3; Gaia Collaboration et al. 2018Collaboration et al. , 2022) ) catalogues.However, there are two stellar systems in V-LoTSS that are not in our radio sample despite being within 100 pc: DG CVn and i Boötis.They are not present in the GCNS likely due to unreliable astrometry and thus their Gaia DR2 source identifier are no longer available in Gaia DR3.
Therefore, we choose to ignore these stellar systems in our analysis.
Not every source in the V-LoTSS × GCNS sample makes an appearance in the LoTSS × GCNS sample, and vice versa.The latter is expected since some radio-bright stellar systems may have fractional polarisation values that makes them drop below the detection threshold of V-LoTSS.In addition, due to time variability of stellar radio sources, a genuine V-LoTSS source that appears in a particular LoTSS field may get "washed out" during mosaicing of the fields in LoTSS-DR2 (Callingham et al. 2023).The mosaicing also explains why even when a stellar system shows up in both the Stokes I and Stokes V catalogues, the quoted Stokes I flux is slightly different in each catalogue.
As for the VLASS × GCNS sample, despite the fact that VLASS and LoTSS have very similar sensitivity and that the VLASS sky fully overlaps with that of LoTSS, very few LoTSS/V-LoTSS detections can be actually found in the VLASS × GCNS sample, and vice versa.The reasons for this discrepancy include the transient nature of stellar radio emission and the vastly different frequencies of the two radio surveys (see Sect. 2.1.3on the importance of survey frequency coverage).Therefore, there is no guarantee that a stellar system that emit in LOFAR band must also emit in VLA band, and vice versa.
There is one source that we remove from the VLASS sample, despite a crossmatching association with GCNS.The source is a M2 dwarf named HD 9770C (Gaia DR3 5022972468944971648), which is part of the visual triple system HD 9770 (also known as BB Scl).Only the C component of this system is in the GCNS, as both the primary A star and secondary B star are not present in the GCNS due to incomplete astrometry in Gaia EDR3 and DR3 (e.g.missing proper motion value).Moreover, the Gaia colour and Gaia absolute magnitude of HD 9770C makes it deviate from the HR main sequence by > 4 mag.All of these peculiar Gaia properties are probably caused by the angular proximity of the three stars: the two stars A and B are in a well-defined 4.559-year orbit with a semi-major axis of 0.171 ′′ , and the AB × C system is in a 111.8-year orbit with a semi-major axis of 1.419 ′′ (Hirshfield & Sinnott 1985).Moreover, both A and B are themselves binaries, with B being an eclipsing binary of the BY Dra type (Watson et al. 2001).Using Gaia DR2 astrometry, we found that the radio emission is actually closer to the AB system than the C component (around 0.3 ′′ vs 1.4 ′′ ).Thus, we conclude that the radio emission most likely stems from the BY Dra variable from the secondary B component of HD 9770, and we remove HD 9770C from our analysis.
As shown in Figs. 2 and 3, both the LoTSS (V-LoTSS included) and VLASS populations can be generally classified into two categories: chromospherically active stellar (CAS) systems and M dwarfs.The CAS systems can be further subdivided into RS Canum Venaticorum (RS CVn) variables and BY Draconis (BY Dra) variables.These are close stellar binaries5 that consist of late spectral types (F to M). Indicated by the presence of strong Ca II H and K emission lines, the stars in these systems have active chromospheres due to strong magnetic field generated by rapid stellar rotation (P rot ∼ days) twisting magnetic flux loops, as the binaries are generally in very close proximity (≲ 0.01 AU) and thus tidally locked.This causes extreme degrees of solar-type activity in these systems, manifested in the form of large stellar spots and magnetic interactions.In the case of the latter, they coherently accelerate electrons (Toet et al.On the other hand, M dwarfs compose the majority of the radio detections, with their spectral type ranging from M1.5 6 to M6 for the LOFAR sample, and M1.5 to M8.5 for the VLA sample.All of these M dwarfs in our V-LoTSS sample were already presented in previous publication (Callingham et al. 2021b), where it was shown that coherent radio emission is ubiquitous across the M-dwarf main sequence, as each spectral type has an equal probability of detection.This implies V-LoTSS has detected both partially convective M dwarfs and fully convective ones.The same applies to both LoTSS and VLASS as well.In particular, VLASS dives into the realm of ultracool dwarfs (UCDs) as the survey detects a few late-M dwarfs as well.This includes the M8.5V ultracool dwarf LSR J1835+3259, which is recently discovered to have a Jupiter-like radiation belt (Climent et al. 2023;Kao et al. 2023).
The underlying reason why M dwarfs are more radio-bright than other isolated stars could be that M dwarfs generally exhibit not just Sun-like activity, but also Jupiter-like properties such as fast rotation, strong large-scale magnetic fields, and auroral radio emission (e.g.Hallinan et al. 2015;Callingham et al. 2021a).Therefore, unlike earlier spectral types, their radio emission could be predicated on the existence of their strong largescale stellar magnetic fields.The M dwarfs should then have a higher efficiency of converting the available energy to radio emission compared to early-type stars, thus making them easier to detect in an untargeted survey.
The only LoTSS/V-LoTSS × GCNS objects that do not fall into these two categories are the prototypical α 2 Canum Venaticorum (α 2 CVn) and HD 220242, the latter of which is an isolated F5 star situated 69 pc away.This makes HD 220242 by far the furthest stellar radio detection from our sample.It is the only main-sequence star other than an M dwarf without known companion stars and is detected with high circular polarisation (78 ± 16).Generally, one does not expect a main-sequence Ftype star to harbour a large scale field strong enough to generate cyclotron maser emission at 144 MHz.As for solar-type bursts powered by plasma emission, the brightest solar burst ever recorded would only be detectable up to ≈10 pc (Vedantham 2020).Therefore, the strong circularly polarised radio emission from such a faraway F-type star is unexpected and more details about this peculiar star are discussed in Appendix B.1.
However, unlike the LoTSS/V-LoTSS × GCNS population, where we mostly detect either CAS systems or M dwarfs, we see a larger variety of stellar systems in the VLASS detections (see Fig. 3).Granted, M dwarfs are still the majority of the detections, along with the chromospherically active stars such as RS CVn and BY Dra systems.However, the VLASS population also consists some unique stellar systems, including a few UCDs, one T Tauri star (TTS), the prototypical α 2 CVn, one chemically peculiar Am star, and a few K-type main-sequence stars which are probably chromospherically active.We describe more details regarding some of the more interesting stellar systems in Appendix B.
One peculiarity regarding the V-LoTSS × GCNS sample is the small amount of V-LoTSS stellar detections outside of the 50-parsec sphere; only 2 out of the 22 sources are beyond 50 pc: EV Draconis which is a RS CVn system and the isolated F star HD 220242.See Appendix C for more details on the lack of detections and on the analysis of the V-LoTSS incomplete- ness.In short, there may be be some incompleteness in the V-LoTSS M-dwarfs population, but the statistical evidence for this is marginal.

Pre-main-sequence stars in the VLASS × GCNS sample
Another peculiarity regarding the VLASS × GCNS sample is the large number of stellar systems that deviate from the main sequence.Shown in Fig. 3, some stellar systems such as RS CVn variables, BY Dra variables, and M dwarfs (including a T Tauri star) are away from the general GCNS population by a few magnitudes of brightness.For the RS CVn systems, this is not surprising since these are all binary systems and some RS CVn systems even host subgiants or giants as one of the components.These make them much more luminous than their main-sequence non-binary counterparts, thus branching off the main sequence.Yet, the fact that so many M dwarfs also exhibit the same phenomenon is unexpected, considering that nothing of such can be seen in the LOFAR sample.A naive explanation could be that these M dwarf systems are unresolved binaries with another star of very similar spectral type, causing the observed luminosity to be up to twice of a typical isolated M dwarf.This alone, however, cannot explain more than half of the M dwarfs with a deviation from the main sequence by more than 0.753 mag, which corresponds to an unresolved binary of two identical stars (see Fig. 4).Since these are M dwarfs, they also cannot be in an evolved state as a typical M-dwarf lifespan far exceeds the current age of the Universe (Laughlin et al. 1997).
An alternative possibility is that these are all young active stars similar to a T Tauri star, and belong to known kinematic groups such as stellar associations, stellar nurseries, and young moving groups (YMGs).Ultimately, we found that around 80% of the stars above the red line shown in Fig. 4 are indeed known members of YMGs and stellar associations, including Beta Pictoris (β Pic) Moving Group (e.g.Shkolnik et al. 2017), Octans-Near Association (e.g.Zuckerman et al. 2013), AB Doradus (AB Dor) moving group (e.g.Malo et al. 2013), and TW Hydrae association (TWA; e.g.Neuhäuser et al. 2010).These stars are all very young (<100 Myr) and are considered to be in the stage when they have not yet reached the main sequence, i.e. pre-main-sequence stars.Therefore, these young stellar objects (YSOs) most likely have extremely fast rotation as they have yet to experience significant rotational spin-down through stellar winds, making them even more magnetically and chromospherically active than the typical M dwarfs due to their convective envelope enabling dynamo processes (Bouvier et al. 2014).This leads to highly variable yet consistent radio emission from YSOs via mechanisms such as gyrosychrotron and ECMI.
The remaining question is why such a YSO population seen in the VLA sample is nowhere to be found in the LOFAR samples.The most likely explanation is that most of the nearby young moving groups (AB Dor, TWA, β Pic, etc.) are located in the Southern Sky and/or the galactic plane, both of which lie outside of the current LoTSS footprint.Another possible explanation is that the predominant radio emission mechanism for these YSOs have most of its power emitted in the decimetre-wave (VLA band) regime rather than metre-wave (LOFAR band).This could happen in the case where most of their radio emissions are gyrosynchrotron emission stemming from flares (Nindos 2020).

Radio evolution with spectral type
As mentioned in Sects.3.1.3and 3.1.4,the vast majority of isolated stars detected in the radio surveys are M dwarfs, so one might be tempted to claim that there is a radio evolution with spectral type and the radio-engine (Sun-like vs Jupiter-like) transition occurs in the realm of M dwarfs.However, M dwarfs are also the most numerous stellar type, and so a proper statistical test is necessary to ascertain this apparent evolution of radio detectability with spectral type.As a statistical test, we compute the cumulative distribution functions (CDFs) with respect to Gaia colours for the 4 populations: the GCNS "background" population, the LoTSS × GCNS population, the V-LoTSS × GCNS population, and the VLASS × GCNS population.We only consider sources within 50 pc of our Solar System in the following analyses since spectral types later than M9 start to become incomplete in GCNS beyond 50 pc (Gaia Collaboration et al. 2021b).The number of sources in GCNS thus decreases to ≈ 40k, and the two sources beyond 50 pc from V-LoTSS × GCNS are removed.The LoTSS × GCNS and VLASS × GCNS samples remain unchanged as we only crossmatch VLASS catalogue with GCNS within 50 pc.
There are many statistical tests widely used for quantifying the discrepancy between two CDFs.Here, we choose the Cramér-von Mises test (CvM test; Anderson 1962) to test the null hypothesis: the two samples being compared are drawn from the same distribution.The CvM test computes a p-value in a two-sample statistical test, which represents how likely the null hypothesis is true.For our analysis, a small p-value signifies a deviation between the radio × GCNS population and the background population, suggesting a transition of radio detection rate with spectral type.We choose the classic threshold of p-value < 0.05 as the condition for a rejection of the null hypothesis.
In the following section, we consider the cases for both the inclusion and exclusion of CAS systems in our CDF analyses.

GCNS background vs radio-detected samples
In Fig. 5, we show the CDFs of the three radio-bright stellar populations (V-LoTSS, LoTSS, and VLASS) compared to the GCNS background CDF, both with and without inclusion of the chromospherically active stars (i.e.RS CVn and BY Dra).The justification for removing the CAS systems is that they likely have a different mechanism powering the radio emission as opposed to isolated stars.We note that with this removal, the GNCS background is practically unchanged due to the rarity of such systems.
We see from the top panel of Fig. 5 that the VLASS × GCNS sample is inconsistent (> 3σ) with the GCNS background distribution, specifically between the spectral types G0 and M3.This shows that RS CVn and BY Dra systems are much more luminous in radio wavelengths compared to any other stellar systems of similar spectral type, likely due to magnetic interactions, tidal locking, and high chromospheric activity (Dulk 1985;  Slee et al. 2008;Toet et al. 2021).Indeed, despite constituting less than 1% of the stellar population (Eker et al. 2008), these CAS systems are so radio-luminous that they dominate the radio source counts.On the other hand, there is agreement between the GCNS background distribution and the LOFAR-detected population within uncertainties.However, since LoTSS and V-LoTSS both detect a significant amount of CAS systems despite their small number (≲200 in our solar neighbourhood), this agreement should only be treated as a mere coincidence due to large uncertainties from the small LOFAR-detected sample size.We conclude that in terms of stellar systems, radio surveys are much more likely to detect these CAS systems than any other stellar systems.
The bottom panel of Fig. 5 shows the radio populations after removal of the known CAS systems.Now, the conspicuous "G0−M3 bump" in the VLASS × GCNS sample is absent, leaving instead a sudden rise in radio detections around the spectral type M3 to M5.The same rise can be seen in the LOFAR sample, albeit with larger uncertainty.Visually, this region is where the VLASS population deviates from the background distribution the most.The CvM test agrees with the discrepancy as well; the p-value for VLASS vs GCNS is 0.0039, whereas the p-values 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 Fig. 6: Zoomed-in version of the bottom panel of Fig. 5, with V-LoTSS × GCNS sample omitted for the sake of clarity.
for V-LoTSS vs GCNS and for LoTSS vs GCNS are 0.039 and 0.198 respectively.

M4 transition in stellar radio engine
The low p-values demonstrate an evolution in radio properties with spectral type.Specifically, Fig. 6 shows where the largest discrepancy between the VLASS (grey line) and background population (blue line) is located.As one can see, while the error bar for the LoTSS population (orange region) is too large to make any definite statements, the grey line deviates from the blue line the most around spectral type M3-M4.This is approximately the spectral type at which most theoretical stellar evolution models predict the transition from partially convective regime (i.e. with a radiative-convective interface akin to the solar tachocline) into the fully convective regime (e.g.Dorman et al. 1989;Clemens et al. 1998;Ribas 2006;Morales et al. 2009).This so-called "M4 transition" (see Stassun et al. 2011) is also supported by observational evidence from the properties of M dwarfs -stellar parameters, activity lifetime, rotation, magnetic field strength, and topology (e.g.West et al. 2008;Donati et al. 2008;Morin et al. 2008).As such, the radio evolution with spectral type may also be a direct cause of such a regime transition.This motivates us to conjecture that the M4 transition also leads to a gradual shift from a Sun-like radio engine to a Jupiter-like radio engine.Originally, we also expect the LOFAR-detected and VLAdetected populations would be inconsistent with each other, with reasons previously mentioned in Sect.2.1.3.However, in both cases where CAS systems are included/excluded, we cannot claim any significant deviations between the radio samples as they all overlap each other within 95% confidence interval.

Comparison to multiwavelength flare rates
One natural question that arises from the detected radio evolution with spectral type is whether the evolution follows other known stellar activity indicators.This leads us to the investigation on the impact of stellar flare occurrence on radio detectability.

GCNS × TESS flare CDF
If flares are ultimately the prime source of radio energy, and if the ratio between the flare energy that goes into the optical band and that into the radio band does not change with spectral type, then we expect the radio detections and the TESS flare detection to have the same CDF shapes.By comparing the two CDF curves, we can therefore determine if stars of some spectral type are more efficient than others at generating radio emission from flaring events.
To compare the CDFs of the optical flare rate and radio detectability, we start with the debiased average TESS flare rates versus T eff from Sect.3.3.2.The CDF of TESS optical flares as a function of spectral type is proportional to the multiplication of two CDF curves: the debiased average-TESS-ensemble-flarerate CDF curve (recall Sects.3.2 and 3.3), and the CDF of the GCNS background distribution.
To cast every CDF curve as a function of T eff , we determine the values of T eff for each star in our sample from the TESS Input Catalog version 8 (TICv8; Stassun et al. 2018) catalogue if available, or from Gaia DR3 catalogue otherwise.If neither catalogue provides a valid T eff , then we estimate T eff from the Gaia colour G BP − G RP using a polynomial fit (from numpy.polyfit;deg = 10) between colour and temperature to the data from Pecaut & Mamajek (2013) 7 .
The result is shown in Fig. 7.Note that here we did not plot the CDF without CAS systems (RS CVn and BY Dra variables) since the TESS targeting strategy is not expected to discriminate against close binaries or active flaring stars, thus there is no compelling reason to exclude them in our analysis.As mentioned 2800 3000 3200 3400 3600 3800 Fig. 8: Same as Fig. 7, but with only M dwarfs considered.Stars before spectral type M0 are excluded from this specific analysis.
during the discussion of Fig. 1d in Sect.3.3.2,early-type stars (earlier than ∼G0) and late-K/early-M dwarfs are much more likely to flare compared to early-to mid-K dwarfs and mid-M dwarfs.Conversely, late-type stars (later than ∼M3) are much more numerous in the Gaia stellar population.Hence, the interplay of these two factors yields a CDF (solid blue curve) not too dissimilar to the CDF of the GCNS population (dashed blue curve), as illustrated in Fig. 7.
Here, the GCNS × TESS flare CDF curve (solid blue line) is visually consistent with all three radio-detected populations at the > 95% confidence level for T eff ≳ 3500 K, i.e. most of the solid blue line is within the shaded regions of the radio × GCNS samples for T eff ≳ 3500 K.More interestingly, it is only within the realm of M dwarfs (later than spectral type ∼M3 in particular) where the GCNS × TESS flare curve CDF starts to significant deviate from all three radio-detected populations.This inconsistency is confirmed by the CvM test as well; the p-values for GCNS × TESS flare vs V-LoTSS, LoTSS, and VLASS are 0.00770, 0.00398, and 0.00013 respectively.
To confirm that the deviation is indeed the most significant in the M-dwarf population seen in Fig. 6, we also compare the CDFs in the T eff range of 3850K to 2800K only, as illustrated in Fig. 8.Such a T eff range corresponds to M0-M6 dwarfs.Now, we can clearly see that the discrepancy between GCNS × TESS flare curve and the radio-detection population is the largest around spectral type M4-M5.Here, the CvM-test p-values are 0.00300, 0.00252, and 3.66 × 10 −9 , for V-LoTSS, LoTSS, and VLASS respectively.These p-values are all much smaller than their counterparts that consider the entire T eff range.
Hence, we can confidently conclude that all three radiodetected population do not fully follow the TESS optical flare statistics.Specifically, the radio samples are consistent with TESS flare CDF for stars with T ≳ 3500 K, but significantly inconsistent for T ≲ 3500 K.This is somewhat unexpected, because one would expect that flares occur due to a universal underlying process which leads to a consistent fraction of energy channelled into radio-emitting accelerated charges.Our results instead suggest that the radio efficiency of optical flares -the fraction of flare energy channelled into radio-emitting chargesevolves with spectral type.In particular, Figs.7 and 8 suggest that the radio efficiency of optical flares remains roughly constant throughout the stellar population until the M dwarfs, where the M4 transition (mentioned in Sect.4.3.2) seems to manifest itself once again, this time in the form of an increase in radio efficiency of optical flares.We hypothesise that this evolution is related to the evolution of magnetic field properties with spectral type owing to the transition from a Sun-like engine to a Jupiterlike engine.
We are aware that, despite our best debiasing efforts described in Sect.3.3, the TESS target selection is still inherently biased since the TESS 50-pc sample used in our analysis is not volume-complete, and thus the degree of biases depends on the topics of interest in the scientific community back during TESS year 1 and 2. To the best of our knowledge, however, a full analysis on the TESS target selection biases does not exist.Therefore, we assume that TESS target selection bias no longer plays a significant role when we only consider the TESS 50-pc sample.

GCNS × X-ray flare CDF
Similarly, we follow the X-ray flare rate trend presented by Johnstone et al. (2021) to construct a correlated CDF from the GCNS background distribution and X-ray flare statistics.Specifically, Johnstone et al. (2021, in bottom-most panel of Fig. 19) showed the evolution of flare rate with total emitted X-ray energies above 10 32 erg vs stellar mass for stars older than 5 billion years old.This gives us the linear relation: N(> 10 32 erg) = 1.11M − 0.12, where N is the number of stellar flares per day and M is the stellar mass in units of solar mass.From this relation, we simply construct the CDF of GCNS background distribution × X-ray flare statistics with very similar procedure from Sect.2.3.1 for the GCNS × TESS flare, but in stellar mass rather than effective temperature.In fact, the procedure is even simpler as there is no need to debias the X-ray flare statistics as previously stated in Sect.2.3.2.
We present the result in Fig. 9.As expected, the GCNS × Xray flare CDF has quite a different shape compared to the GCNS × TESS flare CDF, i.e. the X-ray flare CDF follows a more gently rising trend not seen in the TESS flare CDF.One reason for such a discrepancy could be of astrophysical origin.Audard et al. (2000) found that the rate of X-ray flares is almost linearly proportional to X-ray luminosity among stars of all spectral types.Therefore, solar-mass stars appear to flare in X-ray more frequently than low-mass stars as the former are significantly more X-ray luminous on average.On the other hand, stellar activity is also strongly correlated to rotation period.M dwarfs generally rotate much faster than early-type stars since the latter experience more angular momentum loss owing to stronger stellar winds, causing a spin down of the star.This may imply that there is no universal flare energy partition between optical and X-ray bands, thus leading to disparate observed flare rates between the two wavelengths.
Back to the comparison of X-ray flare CDF and radiopopulation CDF: Fig. 9 shows that early-type stars get significantly boosted in the CDF since early-type stars have higher rate of X-ray flares than the low-mass stars.We can see that visually the GCNS × X-ray flare sample is completely inconsistent with all three radio-detected populations.Using the CvM test, we determine the p-values for GCNS × X-ray flare vs V-LoTSS, LoTSS, and VLASS are 2.72 × 10 −4 , 1.07 × 10 −4 and 8.03 × 10 −8 respectively.Therefore, we can conclude that the radio efficiency of X-ray flares -the fraction of X-ray flare energy channelled into radio-emitting charges -evolves with spectral type as well, similar to the conclusion drawn from the TESS flare CDF analysis in Sect.4.4.1.

Güdel-Benz relationship
We also investigate if our radio-detected samples adhere to the empirical quasi-linear relationship which correlates quiescent gyrosynchrotron 5-GHz radio spectral luminosity L ν,rad and soft X-ray luminosity L X .This is known as the canonical Güdel-Benz relationship (GBR; Güdel & Benz 1993;Benz & Güdel 1994), and the L X ∝ L 0.73 ν,rad law holds for over 10 orders of magnitude in L ν,rad , from the energetic RS CVn systems all the way down to nanoflares from the Sun.Despite their vastly different emission mechanisms, where thermal Bremsstrahlung in the coronal plasma is responsible for the X-ray emission, and incoherent gyrosynchrotron emission from relativistic electrons gyrating in the coronal magnetic field for the radio emission, the presence of this relationship implies there must exist a mechanism that deposits a consistent fraction of flare energy into accelerating charges that emit in the radio and X-ray bands.
However, there are also stellar systems that violate GBR.For example, Hallinan et al. (2008) found that the 5-GHz mildly circularly polarised quasi-quiescent radio emission of UCDs, including brown dwarfs, do not obey GBR.They argued that the radio emission for these systems stems from ECMI instead of gyrosynchrotron, making them radio-overluminous (with respect to GBR) as ECMI is more efficient and capable of creating luminous radio bursts (e.g.Yu et al. 2011).Alternatively, it could be that the energetics of these systems or their magnetospheric structure do not support a stable thermal corona to form, which leads to their X-ray underluminosity (Berger et al. 2010;Williams et al. 2014).Callingham et al. (2021b) found that the highly polarised radio emission at 144 MHz from M dwarfs also violates GBR; these objects are very unlikely to be operating in an incoherent mechanism that is gyrosynchrotron due to their high degree of circular polarisation, and thus a coherent emission mechanism must be responsible for 144-MHz radio emission of these objects.
And so, if one were to accept the explanation that the 5-GHz radio emission from the stellar systems that established the original GBR stems from gyrosynchrotron in the coronal magnetic field, one must expect our V-LoTSS sample to depart from the GBR.Conversely, the VLASS sample should obey GBR, except for stars of very late spectral type -into the realm of UCDswhere their corona starts disappearing.

GBR vs radio-detected samples
Here, we plot the soft X-ray luminosity L X against radio spectral luminosity L ν,rad of the LOFAR-detected and VLA-detected populations, as shown in Figs. 10 and 11 respectively, to see which stellar system obeys GBR and which deviates from it.As previously stated, none of the V-LoTSS sample should obey GBR as these stellar systems operate in coherent mechanisms rather than gyrosynchrotron processes.As expected, such is the case for most M dwarfs and the F-type star; most of these stellar systems deviate from the relationship by a few orders of magnitude.However, the chromospherically active stellar systems in our V-LoTSS sample surprisingly adhere to GBR.This peculiar phenomenon was already discussed by Vedantham et al. (2022), where they report that highly polarised 144-MHz radio emission from RS CVn/BY Dra systems and other high-activity stars unexpectedly follows the GBR.In the other two populations (LoTSS and VLASS), CAS systems also consistently obey GBR, though this is less surprising as previously mentioned.
As for the Güdel-Benz (GB) deviators, we can see from Figs. 10 and 11 that virtually every single M dwarf in all the radio-survey population are consistently radio-bright/X-ray-dim, and most of them by a few orders of magnitude as well.In particular, the M4 transition also seems to manifests itself in the GBR plot, as M dwarfs beyond ∼M4 significantly deviate from the rest of the population on average.Intriguingly, there exist two stellar sources (LSR J1835+3259 and HD 43162C) that deviate from GBR by more than 5 orders of magnitude, and both of them are later than M4 as well.In Fig. 10, there is a clear gap between the early-M-dwarf and the late-M-dwarf populations, especially in V-LoTSS.The deviation from the GBR with respect to spectral type is better illustrated in Fig 12 .Here, we can see that the radio detection becomes significantly more radio-overluminous/Xray-underluminous than the values predicted by GBR in the region of M4-M6.This gap might also signify a bimodality in the dynamo of the M-dwarf population (e.g.McLean et al. 2012;Cook et al. 2014;Williams et al. 2014), supported by geodynamo simulations and Zeeman-Doppler imaging (ZDI) analysis (Morin et al. 2010(Morin et al. , 2011)).
To confirm the monotonic inverse relationship between the L X residual from GBR and spectral type for each of our radio Mdwarf sample in a statistical manner, we calculate the Kendall's τ (Kendall 1938), which is a measure of the relationship between two variables -L X / GB fit and G BP −G RP in our case.This allows us to perform a τ hypothesis test on each of our radio-detected samples, where the alternative hypothesis is that the rank correlation is less than zero, i.e. a negative correlation between L X / GB fit and G BP − G RP .To be succinct, we obtain a p-value much smaller than 0.05 and a negative τ coefficient (around −0.3 to −0.4) for all three radio M-dwarf samples.We thus conclude there exists a negative correlation between goodness of GB fit to the radio data and spectral type in M dwarfs from all three radio surveys.(Williams et al. 2014), where L X is the soft (0.1-2.4 keV) X-ray luminosity and L ν,rad is the 144-MHz total-flux radio luminosity for the LOFAR population.The radio sources without a detectable X-ray luminosity are indicated with a downward-pointing arrow and shown as 3σ upper limits.Colours and symbols are as in Fig. 2, with the additional classification of M dwarfs beyond M4 (represented by yellow circles) motivated by the M4 transition mentioned in Sect.4.3.2.

Rotation-related activities
Even though the M dwarfs are consistently radio-bright or Xray-dim with respect to Güdel-Benz relationship, there seems to be no obvious new scaling law that can be established as these GBR violators generally scatter around the L X against L ν,rad plot in Figs. 10 and 11.Such was the original conclusion from Callingham et al. (2021b); the highly polarised radio-emitting M dwarfs do not establish a clear new relationship even when accounting for chromospheric activity indicators such as Hα luminosity, rotation period, and Rossby number.Therefore, they suggest that the chromospheric activity has limited influence on the detected radio luminosity.In contrast, previous GHz-frequency studies discovered that the fastest rotating M dwarfs have a higher radio surface flux and radio detection fraction (McLean et al. 2012).Here, we attempt to apply similar investigation in our VLA-detected sample.
As mentioned in Sect.4.2, the majority of the non-CAS stellar objects are members of young moving groups (YMGs).This already indicates that they are most likely very fast-rotating stars since they have yet to lose significant angular momentum due to stellar wind.The implication, therefore, is that perhaps the fastest rotators do have a higher radio detection fraction in the 3-GHz band, similar to previous findings by McLean et al. (2012).Unfortunately, finding literature values of rotation period P rot for the VLASS-detected M dwarfs proves to be difficult, as the studies of most of these M dwarfs are lacking in general since they are otherwise run-of-the-mill objects of little interest to the scientific community.Additionally, most of them are too distant and/or too active to be traced by detailed M-dwarf surveys such as CARMENES 8 and MEarth (Irwin et al. 2009).However, for the stellar systems which have measured rotation period (e.g.GJ 3237, GJ 4185, V374 Peg, and V1274 Her), they indeed rotate very fast (typically P rot < 1 day).Therefore, although it is currently not possible to make any definite statements, we believe that rotation period plays an important role on the radio detectability of these stellar objects.

Conclusion
We conducted a search for radio-bright stellar systems using existing radio sky survey catalogues -LoTSS, V-LoTSS, and VLASS -and crossmatched them with the Gaia data along with proper-motion correction.The high sensitivities and resolutions of these radio surveys provided datasets to isolate an untargeted sample of radio stellar systems such that the statistical properties of the population could be determined accurately.We ascertained our stellar radio samples are not dominated by chancecoincidence associations by using the Gaia Catalogue of Nearby Stars.To the best of our knowledge, this is the first untargeted search for ALL stellar radio sources down to sub-mJy level, and thus the first statistical analysis of radio activity trend with spectral type using an unbiased radio-emitting stellar population.
We found that both the LoTSS (V-LoTSS included) and VLASS populations can be generally classified into two categories: CAS systems and M dwarfs.In particular, the VLAdetected sample also contains many more peculiar stellar systems such as YSOs, seemingly isolated K dwarfs, and chemically peculiar stars.The reason why these systems are not present in the LOFAR-detected samples is most likely due to the 3 different survey parameters between LoTSS and VLASS.Firstly, the survey sky coverage: VLASS covers much more of the sky than LoTSS-DR2 (33885 deg 2 vs 5634 deg 2 ), thus it is currently easier to find rare stellar objects in VLASS compared to LoTSS.In addition, VLASS also covers the entire galactic plane visible to the VLA (δ > −40 • ), while LoTSS-DR2 focused on imaging regions at high Galactic latitude and thus deliberately avoids the galactic plane, in which most of the nearby young moving groups (AB Dor, TWA, β Pic, etc.) are located.Secondly, the survey frequency coverage: these two surveys are probing very different radio regimes, as LoTSS is more likely to  11: Same as Fig. 10, but with the VLASS × GCNS population instead.Colours and symbols are as in Fig. 3, with the additional classification of M dwarfs beyond M4 (represented by yellow circles) motivated by the M4 transition mentioned in Sect.4.3.2.The one CAS system that significantly deviates from GBR is a BY Dra system called V1274 Her, which is a mid-M-dwarf binary.Note that two M dwarfs from the VLASS sample are not shown in this figure for the sake of clarity, since they have an X-ray luminosity of ≪ 10 26 erg s −1 .These are LSR J1835+3259 and HD 43162C, with spectral type M8.5 and M5 respectively.capture coherent emission compared to VLASS, as mentioned in Sects.2.1.3and 4.2.The predominant radio emission mechanism for these peculiar systems may have most of its power emitted in the decimetre-wave (VLA band) regime rather than metrewave (LOFAR band).This could happen in the case where most of their radio emissions are gyrosynchrotron emission stemming from flares (Nindos 2020).Thirdly, the survey exposure time: each LOFAR pointing in LoTSS is observed for a total of 8 consecutive hours, whereas VLASS creates 3-minute images.Hence, VLASS is much more susceptible to stellar radio sources with high variability.
Using CDF visualisation and the CvM statistical test, we demonstrated an evolution in radio properties with spectral type.The discrepancy between the radio-detected and the GCNS CDF curves is the largest for the VLASS × GCNS population, simply due to the fact that the VLASS sample size is larger -more than twice the sample size of LoTSS/V-LoTSS.In particular, we showed that radio-detected populations deviate the most from the GCNS background population around spectral type M3-M4.This is remarkable due to its location at the so-called "M4 transition", where most theoretical stellar evolution models predict the transition from partially convective into the fully convective regime.This transition drastically changes the dynamo of stars that is responsible for the generation of magnetic field.As such, we suggest that the radio evolution with spectral type may be related to the changing nature of the stellar dynamo.
We also studied the evolution of radio efficiency with spectral type via comparison to other known stellar activity indicators.In particular, we tested the impact of stellar flare occurrence on radio detectability using the TESS and X-ray flare statistics.The flare statistics in these two wavelengths do not show the same trend: given a particular energy threshold above which flares were counted, we saw the general trend of flare rate decreasing as we moved towards stars of later spectral type in the X-ray flare statistics, and yet the trend seen in TESS flare statistics is much more complicated.Despite that, the radio detectability does not seem to follow neither of these trends.To test if the ratio between the flare energy that goes into the optical/X-ray band and that into the radio band also changes with spectral type, we compared the CDF curves of the radio detection and the TESS/X-ray flare detection.We discovered that the GCNS × X-ray flare curve is clearly inconsistent with all Fig. 12: Residuals from the Güdel-Benz relationship versus Gaia G BP − G RP colour.For the sake of clarity, the V-LoTSS × GCNS sample is omitted, as well as the two M dwarfs from the VLASS sample (grey circle) with residuals ≪ 10 −3 .They are LSR J1835+3259 and HD 43162C, with spectral type M8.5 and M5 respectively.The radio sources without a detectable Xray luminosity are indicated with a downward-pointing arrow and shown as 3σ upper limits.The top axis of the CDF plots indicates the nominal stellar spectral types (Pecaut & Mamajek 2013).
three radio-detected populations at the ≫ 95% confidence level, both visually and statistically.The inconsistency in the GCNS × X-ray flare curve is much more subtle: the GCNS × TESS flare curve agrees with the radio populations throughout most of the stellar population until we reached the M-dwarf population.There, the TESS flare CDF significantly deviates from the radio-population CDFs, especially around spectral type M4-M5.Therefore, the radio efficiency of flares in both optical and Xray -the fraction of flare energy channelled into radio-emitting charges -evolves with spectral type, with later M-dwarfs being more radio-efficient than the rest of the stellar population.We hypothesise that such an evolution may be owing to two factors (or a combination thereof): (i) the efficiency with which charges are accelerated and/or the lifetime of accelerated charges may itself increase for stars of later spectral types; and (ii) the emergence of large-scale axisymmetric magnetic field in fully convective stars may assist more efficient radio emission mechanisms such as ECMI and/or provide new magnetospheric acceleration mechanisms in addition to flare-related processes.We also found that nearly every single M dwarf in all the radio-survey population are consistently and significantly more radio-overluminous or X-ray-underluminous with respect to GBR.In particular, the M4 transition is alluded once again, as M dwarfs beyond ∼M4 significantly deviate from the rest of the population on average.This gap might be signifying a bimodality in the dynamo of the M-dwarf population (e.g.McLean et al. 2012;Cook et al. 2014;Williams et al. 2014), supported by geodynamo simulations and ZDI analysis (Morin et al. 2010(Morin et al. , 2011)).We also confirmed a negative correlation between goodness of GB fit to the radio data and spectral type in M dwarfs from all three radio surveys, i.e. late-M dwarfs deviate from GBR more than early M dwarfs.This was done by performing the Kendall's τ test.We propose follow-up ZDI observations on our radio stel-lar population to see if their radio emission is dominated by the large-scale axisymmetric magnetic field.Additional light curve analysis of the M dwarfs mentioned in Sect.4.2 to get their rotation periods is also vital in helping us prove our conjecture of rotation period playing an important role on the radio detectability of the stellar objects in VLA band.
The fact that radio-detected stellar population is dominated by the CAS systems and M dwarfs regardless of the radio surveys implies there must exist some similarities between these two classes of stars such that they are consistently radio bright.Therefore, we postulate that all of these radio-bright stellar systems exhibit strong large-scale stellar magnetic field.CAS systems and M dwarfs (especially young ones) can generate magnetic field strength up to kilogauss (e.g.Giampapa et al. 1983;Saar et al. 1986;Petit et al. 2004;Kochukhov et al. 2013for CAS systems, Johns-Krull & Valenti 1996;Reiners et al. 2009;Morin et al. 2010;Afram & Berdyugina 2019;Kochukhov 2021 for M dwarfs, andGregory et al. 2012 for YSOs).If the hypothesis of radio predominantly tracing magnetic field dynamics is true, then the predominant radio emission mechanism is likely to be ECMI.Since electron cyclotron maser emits mostly at the fundamental and the second harmonic of the cyclotron frequency: ν B = 2.8B MHz, where B is the local magnetic field strength of the radio emission region (Melrose & Dulk 1982), the main difference between the stellar population in LoTSS and VLASS would be that they have different magnetic field strength; B ∼ 50 G would peak at the LOFAR band, while B ∼ 1 kG at VLA band.
Currently, the LOFAR-detected sample size is severely limited by the incomplete sky coverage of LoTSS/V-LoTSS, as mentioned in Sect 4.2.This makes interpreting the results of statistical test and drawing conclusions from the LOFAR sample more difficult due to small number statistics.By the completion of the LoTSS survey of the Northern Sky, this will effectively increase the current LoTSS/V-LoTSS sample size threefold.Having the new radio stellar detections will most certainly give us definite conclusion on e.g.whether there exists a discrepancy between the LOFAR-detected and VLA-detected stellar population.Moreover, as mentioned in Sect.2.1.2,looking at the Stokes V sky for radio stellar sources is statistically motivated, since the chance-coincidence association is effectively zero in a circularpolarised radio sky survey.Looking at the Stokes V sky is also physically motivated: circular polarisation information may provide a measure of coronal density for stars, and magnetic field strength for brown dwarfs and exoplanets.Therefore, an equivalent of V-LoTSS for VLASS would be very useful for understanding the underlying radio emission mechanism of the VLA radio-bright population.
(α For A) is a F-type solar-age slow-rotating subgiant Fuhrmann & Chini (2015) with a long orbital period of 280 ± 11 years (Izmailov 2019).The angular separation between α For A and B is 4.05 ± 0.08 ′′ (Izmailov 2019).We clearly identify α For B as the radio-loud component in this system because the radio source is < 1 ′′ from B and > 4 ′′ from A in both epochs of VLASS (see Figure B.1). Fuhrmann & Chini (2015) discovered that α For B shows significant chromospheric activity with strong Hα and Ca II H&K lines, whereas nothing alike can be seen on α For A. They argued that since α For A is about solar age and its slow-rotating inactive nature is confirmed by spectroscopy, this means α For B cannot be young and thus has reignited its activity, possibly due to some mass transfer or merger events (Fuhrmann & Chini 2015).Finally, Fuhrmann et al. ( 2016) followed up this system and found a close companion (dubbed α For Bb) around the K dwarf (dubbed α For Ba) using radial velocities.They determined α For Bb to be a white dwarf since there is significant barium enhancement in α For Ba compared to the primary A component.Barium is a red-giant nucleosynthesis product (Jorissen et al. 2019), and thus the material must have been deposited on the K dwarf during a mass-transfer episode that occurred back when the now white dwarf was a giant.This mass transfer likely reignited the K dwarf's activity.At a distance of d = 14.1 pc and with a period of 3.75 days (Fuhrmann et al. 2016), this makes the subsystem α For Ba-Bb one of the nearest white dwarf in a close binary.
An important question, therefore, is the role of the white dwarf (WD) in generating the radio emission of α For B. Currently, there are only two known radio-bright WD systems, excluding a class of variable stars known as cataclysmic variables (Mason & Gray 2007).The first system discovered is AR Scorpii (AR Sco; Marsh et al. 2016), which is a WD/M-dwarf close binary with a period of 3.56 hours.The emission from AR Sco pulses over a broad range of wavelengths (from radio to X-ray) on a period of 1.97 min, corresponding to the beat frequency between the 1.95-min spin period of the white dwarf and the 3.56-hour orbital period, thus earning the moniker "white dwarf pulsar".Very recently, Pelisoli et al. (2023) discovered another such system: J191213.72-441045.1 (J1912-4410) consists of a WD and an M dwarf in a 4.03-hour orbit, exhibiting similar pulsed emission with a period of 5.30 mins.Despite a large number of proposed models (e.g.Geng et al. 2016;Katz 2017;Garnavich et al. 2019) to understand its peculiar observed behaviour and follow-up observations of this system (e.g.Marcote et al. 2017;Stanway et al. 2018), however, the exact radio emission mechanism of AR Sco remains unknown.If α For B turns out to be a WD-pulsar system, this makes it an extremely unique system since localising the radio emission region within the α For Ba-Bb subsystem and detecting its orbital motion is possible with VLBI telescopes, unlike AR Sco and J1912-4410.This in turn would allow greater understanding of the underlying radio emission mechanism, validate the theoretical models, and show whether this is an unique phenomenon or rather pulsing radio emission is ubiquitous for close binaries that contain a WD.
Lastly, α For B being a radio source is not new - Güdel et al. (1995) observed this star with the VLA at 8.4 GHz in two different epochs (1993 and 1995; respective legacy ID: AG394 and AG458).There is a clear detection in each epoch despite strong variability: around 0.6 mJy in 1993 and 0.3 mJy in 1995, which is ∼20σ and ∼10σ respectively.More importantly, the radio emission shows strong circular polarisation (∼ 30%; Güdel et al. 1995) and no significant linear polarisation (< 1%), both of which are very similar to that seen in AR Sco (Stanway et al. 2018).Previously, Güdel et al. (1995) failed to ascertain if the radio emission stemmed from α For A or B. With Gaia astrometry and proper-motion values now available, we reimaged the data and can confidently conclude that the α For B is the radio emitter rather than A.Moreover, the K dwarf also exhibits strong X-ray luminosity (Hünsch et al. 1999;Boller et al. 2016) compared to a typical isolated K dwarf (e.g.Fleming et al. 1994), yet another similarity to AR Sco (Takata et al. 2018).
Here, α 2 CVn is present in both LoTSS (with no significant circular polarisation) and VLASS in both epochs, but it was previously detected in other radio frequencies as well (e.g.Helfand et al. 1999;Drake et al. 2006;Leto et al. 2021).Hajduk et al. (2022) already investigated the radio emission of α 2 CVn detected in LoTSS.They modelled its radio spectrum to constrain the distance at which co-rotation of plasma with the largescale stellar magnetosphere breaks down, and confirmed that the spectral behaviour is consistent with gyrosynchrotron mechanism in α 2 CVn.This is also consistent with previous studies by Leto et al. (2021); Owocki et al. (2022); Shultz et al. (2022), where centrifugal breakout is thought to drive magnetic reconnection, which is responsible the acceleration of radio-emitting gyrosynchrotron.However, the centrifugal breakout model has been challenged due to the lack of Hα emission or shell absorption in α 2 CVn (Pfeffer & McSwain 2022).
One intriguing point regarding α 2 CVn is that it deviates from GBR significantly in both 144-MHz and 3-GHz band; it has a radio spectral luminosity and X-ray luminosity very similar to a non-early M dwarf.This might be suggesting that its predominant radio emission mechanism may also be ECMI similar to that of UCDs.Hajduk et al. (2022) determined the circular polarisation at 144 MHz to be significantly below 60% from the LoTSS detection, so α 2 CVn may have mildly circularly polarised emission instead.UCDs have been observed to emit mildly circularly polarised quasi-quiescent radiation in the past (Berger 2002;McLean et al. 2012;Williams et al. 2014), which is thought to be ECM emission that appears only mildly polarised because of depolarisation effect during propagation in the magnetosphere (Hallinan et al. 2008).α 2 CVn deviates from GBR in both LoTSS and VLASS, but the deviation is slightly more significant in the VLA band.This may imply that its cyclotron emission could be stronger in the GHz regime.This is not surprisingly given the strong organised magnetic field of α 2 CVn (Kochukhov & Wade 2010).Stokes V information of this VLASS source would be very useful in discerning the underlying mechanism of α 2 CVn's radio emission.The contours are of levels with a increment of 4σ, i.e. 4σ, 8σ, 12σ, etc.The peak flux of the source is around 2.5 mJy for Epoch 1 and 2.7 mJy for Epoch 2.

Appendix B.4: TWA 5 A
The furthest radio source in our VLASS × GCNS sample (at 49.7 pc) is discovered to be a stellar system called TWA 5 A (= CD-33deg7795), which is one of the five original kinematic members of the TW Hydrae association (TWA; Kastner et al. 1997): a group of very young (∼ 5 − 15 Myr) low-mass stars and substellar objects (e.g.Neuhäuser et al. 2010).The system consists of the primary T Tauri star (TTS) and a brown dwarf companion (TWA 5 B) located approximately 2 ′′ away (Webb et al. 1999;Lowrance et al. 1999).TWA 5 B is of spectral type M8.5 and has an estimated mass between ∼ 15 − 40M J (Neuhäuser et al. 2000).TWA 5 B is not present in the Gaia catalogues.Their very far (> 100 AU) separation precludes any magnetic interactions between them.The VLASS images in both epochs show that the radio source is located much closer to A than to B (see Fig. B.2). Therefore, we conclude the TTS should indeed be the predominant radio emitter.
The primary TTS itself is a spectroscopic binary (Webb et al. 1999;Muzerolle et al. 2000;Torres et al. 2000).Macintosh et al. (2001) resolved the primary into a 60-mas binary system (TWA 5 Aab) using adaptive optics.Konopacky et al. (2007) and later Köhler et al. (2013) refined the semi-major axis to be 63.7 ± 0.2 mas.Finally, Köhler et al. (2016) improved the orbital parameters and estimated the mass ratio of Aa and Ab to be M Ab /M Aa = 0.93 ± 0.13, making them very similar in mass.Therefore, these two components -Aa and Ab -are also equally likely to be radio emitters.
Intriguingly, there are previous attempts at detecting radio emission from this peculiar system.Osten & Jayawardhana (2006) searched for radio emission from young brown dwarfs in X band, with TWA 5 B being one of them, and did not detect any source in the map at the positions of TWA 5 B in 2005.They also looked at an archival observation made in 2002 and found no detection either.These two non-detections imply a radio luminosity upper limit of this system L ν,rad < 2 − 3 × 10 14 erg s −1 Hz −1 at 8.4 GHz (Osten & Jayawardhana 2006).And yet, TWA 5 A has very significant radio emission (> 10σ) in both epochs of VLASS, giving a radio luminosity L ν,rad ∼ 10 16 erg s −1 Hz −1 at 3 GHz.This radio information thus gives us two models for the radio emission of TWA 5 A: (i) gyrosynchrotron emission that peaks around 3 GHz and quickly drops off beyond the turnover frequency; and (ii) ECM emission which cuts off before the X band.
be recovered given the survey's flux density detection threshold.The ≈100% completeness of V-LoTSS was calculated to be 1.0 mJy by injecting and recovering point sources (Callingham et al. 2023), and thus we use this value as the detection threshold.The simulation gives us the expected ratio between detected sources within and beyond 50 pc.
Finally, to determine how likely it is to only detect 1 CAS system beyond 50 pc given that there are 6 detections within 50 pc, we simply compute the p-value based on a Poisson distribution: where k and λ and the observed and expected number of radio detections beyond 50 pc respectively.Here, k = 1 and we compute λ = 13.3.This gives us a probability of 0.0022%, which is the likelihood that incompleteness can explain the lack of sources beyond 50 pc.Following the same procedures for the case of M dwarfs, however, gives a result consistent with observation: we obtain λ = 0, i.e. we expect no radio detections from M dwarfs outside of the 50-pc sphere.
In addition, we also perform the V/V max statistical test to see if it agrees with the hypothesis that survey incompleteness alone cannot explain the lack of CAS detections.This test (also known as luminosity-volume test) was first introduced by Schmidt (1968) as a measure of the uniformity of the space distribution of radio sources.Originally, the test was to study the cosmic evolution prove that quasars are not at all uniformly distributed in space and the density of radio quasar evolves with distance.However, it is now also widely used also to test for volume completeness of a flux-limited population detected in a survey (e.g.Locatelli et al. 2019).
The test is as follows: let V be the volume enclosed by a sphere with a radius equal to the distance d to a radio source, and V max be the volume enclosed within the maximum distance at which the source would be detectable.Then, for each source, one can calculate the ratio V/V max .If the sources are indeed uniformly distributed and their properties do not evolve with distance, the sample should have V/V max uniformly distributed between 0 and 1 with an average value ⟨V/V max ⟩ = 0.5 in Euclidean space.In a flux-limited survey, this ratio can be expressed as where N is the number of objects in the sample, S i is the flux density of the i th source, and S min is the detection threshold of the survey.This powerful test makes no assumption on the luminosity function, and thus should work on any flux-limited population.Now, we calculate the ⟨V/V max ⟩ for the V-LoTSS sample.For the non-CAS systems, we obtain ⟨V/V max ⟩ = 0.355 ± 0.075, an inconsistency at the ≈2σ level.For the CAS systems, we obtain ⟨V/V max ⟩ = 0.450 ± 0.116 which is consistent with the expected value of 0.5.The uncertainties are obtained via bootstrapping.In conclusion, there may be be some incompleteness in the V-LoTSS M-dwarfs population but the result is marginal.As such, we proceed with the assumption that we are not missing a significant population in our crossmatching.Eruptive* M3.0Ve (10)  18.46 2.44 0.80 ± 0.15 28.53 CR Dra * M1.5Ve (11)  20.15 2.17 0.83 ± 0.23 29.56 GJ 9603 * M2V (12)  21.49RS CVn F2IV+K2IV (24)  46.89 0.65 1.88 ± 0.33 30.75 2MASS J14333139+3417472 * M5V (25)  48.20 3.75 0.95 ± 0.13 27.92 (26)   Notes.). (c) Stokes I flux density S I is given by LoTSS-DR2 (Shimwell et al. 2022). (d) X-ray luminosity X L of each stellar system is given by the Second ROSAT all-sky survey source catalogue (2RXS; Boller et al. 2016), unless specified otherwise.The derivation of the X L upper-limit values for some stellar systems in our sample is described in detail by Callingham et al. (2021b).† Both M-dwarf components of the binary system GJ 856 have very similar angular separation (≲ 1 ′′ ) with the LoTSS radio source.Since the GJ 856 A and GJ 856 B themselves are separated by around 1.25 ′′ , we cannot ascertain which one is (or if both are) emitting the radio signal.Their wide separation (> 10 AU) precludes direct magnetic or chromospheric interactions between them (cf.CAS systems).Eruptive* M3.5V (5)  13.35 2.69 −1.53 ± 0.14 77 ± 10 28.76 YY Gem RS CVn M1Ve+M1Ve (6)  15.08 1.93 −0.87 ± 0.13 47 ± 10 29.88 HAT 182-00605 * M4V (7)  17.90 2.92 −0.73 ± 0.11 82 ± 25 28.53 LP 212-62 * M5V (8)  18.19 3.40 −6.63 ± 0.13 75 ± 3 27.58 (23) GJ 3861

Fig. 1 :
Fig.1: TESS flare statistics and debiasing.The top horizontal axis of each subfigure indicates the nominal stellar spectral types(Pecaut & Mamajek 2013).(a): The average TESS ensemble flare rate -defined as the fraction of flaring stars N flaring /N total in the TESS short-cadence year 1 & 2 observations, weighted by the average flare rate ν TESS -as a function of stellar effective temperature T eff .ν TESS follows the definition from Equation2.The green curve shows the average TESS ensemble flare rate using the entire TESS flare statistics which suffers greatly from TESS sampling bias (see Sect. 3.3.1),whereas the orange curve only considers TESS stars within 50 pc, with the blue curve as its polynomial fit (from numpy.polyfit;deg = 5) used for subsequent analysis.M dwarfs later than ∼M6 were not observed in a high enough sample size due to TESS target selection (which avoids very faint stars) and are thus excluded (represented by the red region).(b): Stellar radius R * as a function of T eff .The Basti-IAC isochrone model(Hidalgo et al. 2018) and the BHAC15 isochrone model(Baraffe et al. 2015) are combined to create our interpolated model.(c): Minimum flare energy E f,min detectable by TESS on a particular star with effective temperature T eff .The unit for the flare energy is arbitrary since we are only interested in the shape of the E f,min versus spectral type curve.(d): Same as Fig.1a, but with the debiased flare rate ν instead of TESS-observed flare rate ν TESS .This debiased average-TESS-ensemble-flare-rate curve is obtained by multiplying the two blue curves in Figs.1a and 1c.

Fig. 2 :
Fig. 2: Hertzsprung-Russell (HR) diagrams for the V-LoTSS × GCNS sample (left panel) and LoTSS × GCNS sample (right panel), according to Gaia G BP − G RP colour and Gaia G absolute magnitude.The radio-bright stellar systems are represented by different colours and symbols according to object classification, as shown in the legends.There exist two sources in the V-LoTSS × GCNS population that are beyond 50 pc, whereas only radio detections within 50 pc are included in the LoTSS × GCNS population to suppress chance-coincidence associations.The top axis of the HR diagrams indicates the nominal stellar spectral types (Pecaut & Mamajek 2013).

Fig. 3 :
Fig. 3: Hertzsprung-Russell (HR) diagram for the VLASS × GCNS sample according to Gaia G BP −G RP colour and Gaia G absolute magnitude.The radio-bright stellar systems are represented by different colours and symbols according to object classification, as shown in the legends.Green triangle represents stellar systems that are likely to be chromospherically active from literature.The radio sources here are all within 50 pc, as crossmatching the entire GCNS would lead to significant coincidence associations with VLASS sources.The top axis of the HR diagrams indicates the nominal stellar spectral types (Pecaut & Mamajek 2013).

Fig. 4 :
Fig. 4: Top panel: Zoomed-in version of Fig. 3 with RS CVn variables removed.The cyan line represents the median fiducial and the red line shows the same fiducial shifted by -0.753 mag, which corresponds to an unresolved binary system of two identical stars.Therefore, stellar systems that significantly above the green line cannot be explained simply by binarity.The top axis of the HR diagrams indicates the nominal stellar spectral types (Pecaut & Mamajek 2013).Bottom panel: Residuals of M G relative to the median fiducial (cyan).

Fig. 5 :
Fig.5: Cumulative distribution functions (CDFs) for the GCNS background distribution (dashed blue), the V-LoTSS × GCNS sample (green), the LoTSS × GCNS sample (orange), and VLASS × GCNS sample (grey), shown as a function of Gaia colour G BP − G RP .The top panel shows the two radio-detected populations with the inclusion of chromospherically active stellar (CAS) systems.The bottom panel shows the population without them.The shaded regions correspond to the 95% confidence level based on a binomial distribution.As the number of sources in GCNS is much larger than the radio-detected populations, the confidence interval of GCNS is negligible.The top axis of the CDF plots indicates the nominal stellar spectral types(Pecaut & Mamajek 2013).

Fig. 7 :
Fig.7: Cumulative distribution functions (CDFs) for the GCNS background distribution (dashed blue), the GCNS × TESS flare distribution (solid blue), the V-LoTSS × GCNS sample (green), the LoTSS × GCNS sample (orange), and VLASS × GCNS sample (grey), shown as a function of effective temperature T eff .The solid blue line is obtained from the product of the dashed blue line and TESS flare rate curve according to Fig.1d.The shaded regions correspond to the 95% confidence level based on a binomial distribution.Note that the stars with T eff < 2800 K (around spectral type M6; shaded red box) are excluded in our analysis due to their incompleteness in TESS.The top axis of the CDF plots indicates the nominal stellar spectral types(Pecaut & Mamajek 2013).

Fig. 9 :
Fig.9: Cumulative distribution functions (CDFs) for the GCNS background distribution (dashed blue), the GCNS × X-ray flare distribution (solid blue) according to the study byJohnstone et al. (2021) , the V-LoTSS × GCNS sample (green), the LoTSS × GCNS sample (orange), and the VLASS × GCNS sample (grey), shown as a function of stellar mass in units of solar mass M ⊙ .The shaded regions correspond to the 95% confidence level based on a binomial distribution.Note thatJohnstone et al. (2021) only consider stars of stellar mass between 0.1 and 1.2 M ⊙ , and so we exclude stars with less than 0.1 M ⊙ (around spectral type M6) in this specific analysis.The top axis of the CDF plots indicates the nominal stellar spectral types(Pecaut & Mamajek 2013).

Fig. 10 :
Fig. 10: Radio-detected populations (V-LoTSS on the left panel, LoTSS on the right panel) plotted against the canonical Güdel-Benz relationship (GBR) represented by the dashed blue line: L X = 9.48 × 10 18 L 0.73ν,rad(Williams et al. 2014), where L X is the soft (0.1-2.4 keV) X-ray luminosity and L ν,rad is the 144-MHz total-flux radio luminosity for the LOFAR population.The radio sources without a detectable X-ray luminosity are indicated with a downward-pointing arrow and shown as 3σ upper limits.Colours and symbols are as in Fig.2, with the additional classification of M dwarfs beyond M4 (represented by yellow circles) motivated by the M4 transition mentioned in Sect.4.3.2.

Fig
Fig. B.1: VLASS Epoch 1 (left panel) and Epoch 2 (right panel) images using CIRADA cutout service 10 , centred at the corresponding radio source.In both images, the lime dot represents the Gaia position of α For B with proper-motion correction, and the blue dot represents that of α For A. The FWHM beam size is shown in the bottom left corner of each image as a grey crossed circle.The rms noise σ is around 140 µJy in both epochs.The contours are of levels with a increment of 4σ, i.e. 4σ, 8σ, 12σ, etc.The peak flux of the source is around 3.3 mJy for Epoch 1 and 1.5 mJy for Epoch 2.

Fig
Fig. B.2: VLASS Epoch 1 (left panel) and Epoch 2 (right panel) images using CIRADA cutout service 11 , centred at the corresponding radio source.In both images, the lime dot represents the Gaia position of TWA 5 A with proper-motion correction, and the blue dot represents the position of TWA 5 B using astrometry by Neuhäuser et al. (2000).The FWHM beam size is shown in the bottom left corner of each image as a grey crossed circle.The rms noise σ is around 150 µJy in Epoch 1 and 130 µJy in Epoch 2.The contours are of levels with a increment of 4σ, i.e. 4σ, 8σ, 12σ, etc.The peak flux of the source is around 2.5 mJy for Epoch 1 and 2.7 mJy for Epoch 2.
(a) Asterisk * denotes M dwarfs unless specified otherwise."Eruptive*" denotes eruptive variable stars classified by SIMBAD.(b)Distance d (derived from the parallaxes) and Gaia colour G BP − G RP of each stellar systems are given by the Gaia Catalogue of Nearby stars(Gaia  Collaboration et al. 2021b) which quotes measurements from Gaia Early Data Release 3(Gaia Collaboration et al. 2021a Notes.(a)Asterisk * denotes M dwarfs unless specified otherwise."Eruptive*" denotes eruptive variable stars classified by SIMBAD.(b)Distance d (derived from the parallaxes) and Gaia colour G BP − G RP of each stellar systems are given by the Gaia Catalogue of Nearby stars(Gaia  Collaboration et al. 2021b) which quotes measurements from Gaia Early Data Release 3(Gaia Collaboration et al. 2021a

Table 1 :
Parameters of the radio surveys utilised for identifying radio-bright stellar systems.

Table D .
1 contains the definitions of all the acronyms frequently used in this work.

Table D .
1: Table of acronyms frequently used in this work.Table A.1: Multi-wavelength properties of the 144-MHz detected LoTSS × GCNS population sorted by increasing distance d.The symbols G BP − G RP , S I , and log L X correspond to Gaia BP-RP colour, Stokes I flux density, and logarithmic soft (0.1-2.4 keV) X-ray luminosity respectively.

Table A .
2: Multi-wavelength properties of the 144-MHz detected V-LoTSS × GCNS population sorted by increasing distance d.The symbols G BP − G RP , S V , S V /S I , and log L X correspond to Gaia BP-RP colour, Stokes V flux density, fraction of circularly polarised radio emission, and logarithmic soft (0.1-2.4 keV) X-ray luminosity respectively.

Table A .
3: Multi-wavelength properties of the 3-GHz detected VLASS × GCNS population sorted by increasing distance d.The symbols G BP − G RP , S I , and log L X correspond to distance, Gaia BP-RP colour, Stokes I flux density, and logarithmic soft (0.1-2.4 keV) X-ray luminosity respectively.The last column refers to the VLASS epoch from which the radio data is used for this study, i.e. whether the value of S I is quoted from VLASS Epoch 1 or 2.