© The Authors 2025

1 Introduction

Coronal mass ejections (CMEs) are violent eruptions of plasma from a stellar corona. On the Sun, they occur on average about once every few days (Gopalswamy et al. 2024). The bulk motion of the plasma, when super-Alfvénic, generates a shock that accelerates charges. The accelerated charges eventually cause radio emission via the plasma mechanism (Dulk 1985). The emission occurs at the local plasma frequency (ν_p), given in Gaussian units as $${\nu }_p={\left(\frac{n_e{e}^2}{\pi m_e}\right)}^{1/2}\mathrm{H}\mathrm{z},$$(1)

where n_e is the electron number density, and m_e is the mass of an electron.

As the super-Alfvénic shock propagates away from the star, the ambient plasma density decreases significantly from values typical of a corona (around 10⁸−10¹⁰ cm⁻³). The emission frequency (ν_p ≈ 100–1000 MHz at these densities) decreases, forming the characteristic time-frequency drift of a so-called Type II radio burst (Payne-Scott et al. 1947; Wild & McCready 1950). For the Sun, the connection between drifting radio bursts and CMEs has been confirmed by numerous coincident coronagraph and radio observations (e.g. Wild & Smerd 1972; Kahler 1992; Webb & Howard 2012), with typical drift rates of ~ −0.4 MHz/s at 150 MHz (Aguilar-Rodriguez et al. 2005a) and burst durations of several minutes (Roberts 1959; Nelson & Melrose 1985).

Coronal mass ejections are major contributors to space weather and may influence planetary habitability. They can strip planetary atmospheres, affecting their long-term stability (e.g. Kulikov et al. 2007; Lammer et al. 2013), yet may also stimulate chemical reactions in planetary atmospheres that form prebiotic molecules (Scalo et al. 2007; Lingam et al. 2018). In particular, evidence suggests that the CME event rate or associated energies are elevated in certain classes of stars; if true, this could profoundly impact our understanding of the habitability of their exoplanets (e.g. Glassgold et al. 1997; Turner & Drake 2009).

To date, there has been no unequivocal detection of a stellar CME. Existing methods, such as observing extreme UV/X-ray coronal dimming (e.g. Veronig et al. 2021; Loyd et al. 2022; Namekata et al. 2024), blueshifts in chromospheric lines (e.g. Vida et al. 2019; Namekata et al. 2021; Notsu et al. 2024), blueshifts in transition lines (e.g. Leitzinger et al. 2011; Argiroffi et al. 2019; Inoue et al. 2024), and X-ray/UV absorption (e.g. Haisch et al. 1983; Bond et al. 2001), have revealed signatures consistent with stellar CMEs; however, we lack conclusive evidence that the plasma responsible for these signatures remains confined to the stellar magnetosphere or has dynamically decoupled from it. In contrast, Type II radio bursts provide a more direct signature: they imply the existence of bulk plasma motion at super-Alfvénic speeds. At these speeds the plasma cannot be confined by the stellar magnetic field, implying that the plasma producing the shock was ejected into interplanetary space.

Until recently, observational efforts to detect stellar Type II radio bursts had only yielded non-detections (e.g. Crosley & Osten 2018; Villadsen & Hallinan 2019; Zic et al. 2020; Callingham et al. 2021b). However, Callingham et al. (2025) and Tasse et al. (2025) have now reported the first detection of a stellar Type II burst, made possible by recent advances in low-frequency radio astronomy, particularly through the unprecedented sensitivity, frequency coverage, field of view, and resolution of the LOw-Frequency ARray (LOFAR; van Haarlem et al. 2013), which enables the detection of more bursts and more stringent upper limits to be placed on their occurrence rate. Low-frequency observations (below ~500 MHz) are especially critical, as the anticipated plasma densities in stellar coronae will yield plasma emission at these frequencies (Crosley & Osten 2018; Villadsen & Hallinan 2019).

In this work, we present the largest and most comprehensive search for stellar Type II radio bursts to date, identifying bursts with time-frequency drifts in time-frequency spectra (so-called dynamic spectra) synthesised by the post-processing code DynSpecMS¹ (Tasse et al. 2025). Tasse et al. (2025) present the 233 899 time-frequency spectra for 93 608 unique stars within 100 parsecs, derived from calibrated interferometric data in the LOFAR Two Meter Sky Survey (LoTSS; Shimwell et al. 2017, 2022). The 100-parsec boundary is motivated by detector sensitivity and computational constraints. At this distance, our derived 2.5 mJy 3 σ sensitivity limit for 1-minute bursts corresponds to a radio luminosity of 3 × 10¹⁶ erg s⁻¹ Hz⁻¹, several orders of magnitude brighter than observed solar Type II bursts.

To improve detection sensitivity and reduce false positives, particularly those caused by radio frequency interference (RFI) and sidelobe structures of nearby bright sources, we searched for bursts via their circularly polarised emission (Stokes V). Solar Type II radio bursts are typically weakly circularly polarised (0–40%; e.g. Komesaroff 1958; Stewart 1966; Fomichev & Chertok 1968), although instances of strongly circularly polarised emission are theoretically plausible (Thejappa et al. 2003). While Stokes V data reduce our sensitivity to weak circularly polarised emission, they significantly enhance detection reliability. This is because Stokes V is far less affected by the sidelobe structures of other sources compared to the total intensity, given that Stokes V sources are inherently scarcer (Lenc et al. 2018; Callingham et al. 2023).

In this paper, we address whether stellar Type II radio bursts occur as frequently as their solar counterparts. We begin by detailing our observational setup and data reduction strategy in Sect. 2. We then describe our methodology for identifying and characterising potential burst candidates in Sect. 3. In Sect. 4 we analyse the properties of the detected stellar bursts and establish upper limits on the Type II burst rate. Finally, we interpret our findings in Sect. 5 and summarise our key conclusions in Sect. 6.

2 LoTSS data and post-processing

The methodology used to generate the dynamic spectra used in this paper is described in Tasse et al. (2025); here we summarise the survey characteristics for context: LoTSS is an ongoing survey of the northern sky through eight-hour exposures each of 3168 pointings at 120 to 168 MHz (Shimwell et al. 2017). The dataset encompasses all data processed up to May 2023, comprising 1651 pointings, which totals a combined coverage of 12 500 deg². This includes observations from both LoTSS Data Release 1 and Data Release 2, as well as additional data beyond these releases, forming part of the upcoming Data Release 3 (DR3). The standard LoTSS pipeline² produces high-fidelity, thermal noise-limited images and source-catalogues using direction-dependent self-calibration (Tasse 2014; Smirnov & Tasse 2015; Tasse et al. 2018, 2021). With DynSpecMS, the clean source components are subtracted from the visibilities using their direction-dependent gains. The residual visibilities are then coherently beam-formed towards our stellar targets to produce full-Stokes dynamic spectra at the native time–frequency resolution of the processed visibilities (8 seconds and ~100 kHz). Per pointing, ~100 dynamic spectra are also generated on ‘blank-sky’ locations near the stellar targets, away from any radio sources, to cross-check residual systematic errors from faint interference and calibration issues. A typical LoTSS pointing contains ~100 stars within 100 pc, resulting in 100 × 8 h worth of stellar dynamic spectra data. In total, 93 608 unique stellar objects were observed, each for approximately 8 hours, though some pointings were observed multiple times.

The full observational campaign totals 1 264 977.74 hours, or 144.31 years. However, this total exposure time includes stellar objects like white dwarfs, which are unlikely to host CMEs. To calculate a more representative exposure time, we refined our Gaia sample by applying specific data quality filters: we included all Gaia DR3 stellar objects within 100 pc (Gaia Collaboration 2023) and we adopted a 10% relative precision criterion: parallax_over_error > 10. Similarly, we applied filters to the relative flux error on the G_BP, and G_RP photometry: phot_rp_mean_flux_over_error > 10, and phot_bp_mean_flux_over_error > 10. We included only stars with high-quality astrometric solutions: visibility_periods_used > 6. Finally, as outlined by Gaia Collaboration (2018), we used an additional filter to limit our analysis to sources within the empirically defined locus of the I_BP + I_RP/I_G fluxes ratio as a function of G_BP − G_RP colour: phot_bp_rp_excess_factor > 1.0+0.015 (G_BP − G_RP)² and phot_bp_rp_excess_factor < 1.3+0.06 (G_BP − G_RP)².

We removed white dwarfs from the total sample by iteratively grouping data points into two distinct clusters based on their proximity in the colour–magnitude space and then excluding the cluster corresponding to white dwarfs. The final sample, shown in Fig. 1, focuses specifically on potential CME hosts and has been observed for a total of 107.21 years. Note that we processed and analysed the entire 144.31 years of data, but we used the 100-year duration when calculating total CME event rates.

Fig. 1 Colour-magnitude diagram of Gaia sources observed with LoTSS, applying the quality filters described in Sect. 2. The sources are categorised into two clusters based on their proximity in the colour–magnitude space.

3 Search for and discovery of bursts

We searched for drifting bursts by analysing the Stokes V time series data. We averaged over the frequency axis after correcting for different trial sweep rates. To mitigate the effects of residual RFI, the frequency averaging uses inverse-variance weighting, where each data point is weighted by the inverse of its blank-sky variance.

3.1 Inverse-variance weighted time series

The inverse-variance weighted average (${\hat{y}}_t$) at a given time (t), after drift-correction, is calculated as $${\hat{y}}_t=\frac{\sum _i{y}_i/{\sigma }_i^2}{\sum _{i}1/{\sigma }_i^2},$$(2)

where y_i is the flux density of the target frequency channel (i), and σ_i² is the variance for frequency channel i at time t, which is calculated using the blank-sky spectra as $${\sigma }_i^2=\frac{\sum _j{w}_j{(b_{j,i}-{\overline{b}}_i)}^2}{\sum _j{w}_j},$$(3)

where b_j,i is the flux density of blank-sky spectrum j in frequency channel i, w_j is the weighting factor for blank-sky spectrum j, and ${\overline{b}}_i$ is the weighted mean: $${\overline{b}}_i=\frac{\sum _j{w}_j{b}_{j,i}}{\sum _j{w}_j}.$$(4)

Due to RFI and calibration errors, some blank-sky pixels appeared excessively bright, causing inverse-variance weighting to down-weight the target pixel, even when it was unaffected. To address this, we implemented several strategies to adjust the contribution of each individual blank-sky spectrum so that the combined ${\sigma }_i^2$ best represents the target data: $$w_j=\{\begin{array}{ll}\frac{N_{\text{off-facet}}-N_{\text{in-facet\hspace{0.17em}}}}{N_{\text{in-facet}}}& \text{\hspace{0.17em}if spectrum\hspace{0.17em}}j\text{\hspace{0.17em}is in-facet\hspace{0.17em}}\\ 0& \text{\hspace{0.17em}if spectrum\hspace{0.17em}}j\text{\hspace{0.17em}is in top\hspace{0.17em}}10\mathrm{\%}\text{\hspace{0.17em}outliers\hspace{0.17em}}\\ 1& \text{\hspace{0.17em}otherwise (off-facet, non-outlier).\hspace{0.17em}}\end{array}$$(5)

With this ancillary weighting and flagging, we obtained an inverse-variance weighted time series in which RFI is largely suppressed.

3.2 Drift-corrected search

We constructed the frequency averaged time series by applying a range of trial drift corrections to the dynamic spectrum, effectively compensating for frequency-dependent time delays. To achieve this, we adopted a simplified drift model in which a shock propagates at a constant speed through a medium where the density (n) decreases as $n\propto \frac{1}{r^2}$, where r is the radial distance from the star. Since the plasma frequency follows ν_p ∝ n^1/2, this implies ν_p ∝ 1/r. We can express the radial position of a shock moving at constant radial speed as r = r₀ + vt, where r₀ is the starting position, v the speed of the shock, and t the time. This leads to $$\frac{d{\nu }_p}{{\nu }_p}=-\frac{dr}{r}=-\frac{\nu }{r}dt.$$(6)

This in turn gives $$\frac{dt}{d{\nu }_p}=-\frac{C}{v{\nu }_p^2},$$(7)

where C is a constant, arising from ν_p = C/r. If we integrate from an initial frequency ν to +∞, we obtain the arrival time: $$t\left({\nu }_p\right)=\frac{C}{v{\nu }_p}.$$(8)

The time delay across our observing band is therefore $$\mathrm{\Delta }t=\frac{C}{v{\nu }_{min}}-\frac{C}{v{\nu }_{max}}.$$(9)

Using this time-delay model, similar to pulsar de-dispersion techniques, we could ‘unsweep’ the drift of the burst in the dynamic spectrum, counteracting the delay on a per-channel basis and aligning the burst vertically in the dynamic spectrum. Next, we computed the inverse-variance weighted time series on these drift-corrected spectra, enhancing the S/N for a drift that effectively aligns the burst vertically. The drift value that achieves the best vertical alignment, i.e. the highest S/N, represents the true drift of the burst.

We searched for bursts by setting the range of trial drifts to extend up to 30 minutes, as Type II bursts do not often last longer (Kumari et al. 2023). Using a range of trial drifts up to 30 minutes results in 225 times the original number of time series (30 minutes at an 8-second resolution), each corrected for progressively larger drifts. We applied a boxcar matched-filter search with widths ranging from 8 to 600 seconds to each drift-corrected time series, searching for events with S/N≥8 to minimise false-positive detections that are generated by Gaussian noise. Our data consist of 2.73 × 10¹¹ data points, which we sampled an additional 225 times with our drift-corrected search. Gaussian noise is expected to produce an 8-sigma event approximately once every 1.61 × 10¹⁵ samples, resulting in fewer than one false-positive detection due to Gaussian noise. Since the noise in our data may not be perfectly Gaussian, an 8-sigma threshold should at least significantly reduce false detections.

To test if our new method can reliably detect drifting bursts in the dynamic spectrum, we simulated a drifting burst by modelling a stellar density profile using the Parker Solar wind model (Parker 1965), through which a shock propagates (see Appendix B). As the shockwave passes through the wind and emits via plasma emission, the emission frequency will be directly related to the density profile. Assuming a base wind electron density of 10⁹ cm⁻³, a coronal temperature of 3 MK, and a fractional bandwidth of 0.25 (Aguilar-Rodriguez et al. 2005b), a shockwave moving at a constant velocity of 1000 km s⁻¹ could produce a signal in our dynamic spectra as shown in Fig. 2a. We defined the sign of the Stokes V emission as left-hand circularly polarised light minus right-hand circularly polarised light (Callingham et al. 2021a,b). The drift-corrected time series can be visualised in a ‘drift-time’ spectrum, which shows the inverse-variance weighted frequency-integrated burst brightness as a function of drift and time. For the Parker Solar wind burst, this is illustrated in Fig. 2b. In this parameter space, any broadband signal, whether drifting or not, gives rise to a symmetric ‘bow-tie’-like shape.

To further confirm that the detected bursts are of stellar origin and are not created by RFI, we created Stokes V images with a resolution of 20” using WSClean (Offringa et al. 2014a,b), where we imaged along the identified drift; only using the data within the time periods covered by the drift, flagging the rest. From the image, we can determine if the source of the emission originates from a point source in the sky at the location of the star. However, this imaging process is computationally expensive, especially for a large number of burst candidates. Therefore, we first manually inspected each burst to reject candidates created by RFI or sidelobe patterns of nearby bright sources, which often exhibit features such as negative total intensity burst emission or non-physical, rapidly changing intensities across both time and frequency channels. We also rejected any candidate for which the Stokes I emission is less than half of the emission in Stokes V. Only after passing these checks did we image the burst to identify if the emission originates from a point source.

Fig. 2 Panel a: mock Stokes V dynamic spectrum of a radio burst due to plasma emission modelled using the Parker Solar wind model with solar parameters, assuming a shock speed of 1000 km s−1 and a S/N of 30. Panel b: frequency-integrated burst brightness as a function of the drift delay and time, modelled using the Parker Solar wind model with solar parameters, assuming a shock speed of 1000 km s−1. The defining ‘bow-tie’-shaped feature of a broadband swept astrophysical signal is clearly visible. The red cross indicates the brightest pixel and its corresponding burst drift.

4 Burst properties

With our drift-search method, we can also identify bursts with zero drift, as their S/N will peak at zero drift-correction. Using our approach, we detected 21 stellar radio bursts: 19 without an observable drift, and 2 with drifting features. Table 1 shows the stellar hosts and the time of arrival for the 21 detected stellar bursts. Additionally, to provide an indication of the stellar activity levels, we included literature values of X-ray fluxes converted to X-ray luminosities for each star, as CMEs are often associated with enhanced X-ray emission (Sheeley et al. 1975). Where published X-ray fluxes were unavailable, we estimated upper limits using data from the ROentgen SATellite (ROSAT; Truemper 1982) All-Sky Survey (Voges et al. 1999). Figure 3 shows the dynamic spectra of the stellar bursts without any drifting features, where the spectra were smoothed to optimise clarity. The two drifting bursts we detect are shown in Figs. 4 and 5. Burst B20 corresponds to the event reported by Callingham et al. (2025).

4.1 Drift characterisation

To obtain statistical estimates for the drift parameters, we applied a Markov chain Monte Carlo (MCMC) algorithm to fit a drift model to the data, which can be described by four parameters: the constant burst flux density S₀, the arrival time t₀, the drift rate a, and the burst width Δt. We again modelled the time delay as Eq. (9). A drifting burst model takes the form $$S_k=S_{0}W(t_0+\frac{a}{{\nu }_k},\mathrm{\Delta }t),$$(10)

where k is the frequency channel number and W is given by $$W(t,\mathrm{\Delta }t)=\{\begin{array}{ll}1& ;-\mathrm{\Delta }t/2<t<\mathrm{\Delta }t/2\\ 0& ;\text{\hspace{0.17em}otherwise}.\end{array}$$(11)

The MCMC approach maximises the sum of the log-likelihood function (ℒ), defined by $$\mathcal{L}=C-\sum _k\frac{\mid \text{data}-\text{model}{|}^2}{2{\sigma }_k^2},$$(12)

where C is a constant, k a specific frequency channel, and ${\sigma }_k^2$ the variance in that channel. For our drift model, the likelihood is expressed as $$\text{\hspace{0.17em}Likelihood\hspace{0.17em}}=\prod _k\exp (-{[S_k-S_{0}W(t_k+t_0+\frac{a}{{\nu }_k},\mathrm{\Delta }t)]}^2/2{\sigma }_k^2).$$(13)

Thus, the log-likelihood function becomes $$\mathcal{L}=C-\sum _k\frac{{|S_k-S_{0}W(t_k+t_0+\frac{a}{{\nu }_k},\mathrm{\Delta }t)|}^2}{2{\sigma }_k^2}.$$(14)

We derived the posterior probability distributions using the nested sampling Monte Carlo algorithm implemented in the UltraNest package³ (Buchner 2021). UltraNest combines several advanced methods for efficient nested sampling, including the mlfriends algorithm (Buchner 2016, 2019), which adopts ellipsoidal sampling regions based on the local geometry of the likelihood distribution. It also incorporates strategies for dynamic live point management, robust integration with bootstrapped uncertainties, and adaptive termination criteria to ensure accurate results. Uniform priors are adopted for all parameters except for S₀, which is assigned a logarithmic prior, meaning log(S₀) is uniformly distributed. Bursts B20 and B21 show a significant drift.

The characteristics of these two drifting bursts are shown in Table 2. Errors are given by the 95% confidence interval on the posterior distribution.

4.2 Drifting burst properties

Burst B20, shown in Fig. 4, and burst B21, shown in Fig. 5, both show a time-frequency drift rate: $-{0.801}_{+0.003}^{-0.003}$ and − ${0.06}_{+0.002}^{-0.002}$ MHz s⁻¹, respectively. We did not analyse burst B20 in detail as it has been covered by Callingham et al. (2025).

The modelled burst flux for B21, integrated over the entire burst, is ${4.5}_{-1.3}^{+1.4}$ mJy in Stokes V, which places it at the edge of our detection sensitivity limit. The burst originates from stellar object LP 215-56, (Gaia DR3 ID 769964666964655488), an M2.4V dwarf star (Birky et al. 2020), located at a distance of 62.13 pc (Gaia Collaboration 2020). The burst is significantly circularly polarised, with our measurements indicating a polarisation fraction of greater than 50%. The presence of a nearby (<18′) 1.3 Jy radio source imposes dynamic range limitations in the Stokes I data, preventing us from determining a more precise polarisation fraction beyond this lower limit of 50%. The star is likely part of a binary system, as its Gaia re-normalised unit weight error is 5.25; values larger than 1.4 indicate binary systems (Lindegren et al. 2018).

The star has been observed at optical wavelengths by the Transiting Exoplanet Survey Satellite (TESS). We searched the TESS 2-minute light curve data (TIC 253052423; sectors 22, 48, and 49) and find no flares, which allows us to place an upper limit on the flare rate of <1 flare yr⁻¹ for flare energies ≥10³³ erg, assuming a typical power law slope for the flare frequency distribution of −2 (Ilin et al. 2021). The upper limit energy was determined by defining the smallest detectable flare as three consecutive data points 3σ above the noise level of the de-trended light curve. From the so derived minimum detectable equivalent duration (Gershberg 1972), we derived a bolometric energy following the procedure inShibayama et al. 2013, and assuming a blackbody temperature of 10 000 K. This gives a detection threshold of 3 × 10³² erg in a 77 d long light curve. The Hungarian-made Automated Telescope Network (HATNet; Hartman et al. 2011) Exoplanet Survey reports a rotational period of ~44 days (Bakos 2018). However, we are unable to confirm this with the TESS light curve data due to the limited observation time.

As the binary status of the star is uncertain, we considered two scenarios: (i) the star is not in a binary system and (ii) the more likely scenario is that the star is in a binary system. If the star is not in a binary system, we estimated its mass and radius based on its position in the Gaia colour–magnitude diagram: M = 0.507 ± 0.020M_⊙ and R = 0.509 ± 0.015R_⊙ (Stassun 2019). Alternatively, if the star is in a binary system, we made the rough assumption that our target star has half the mass M = 0.25 M_⊙ and an accompanying radius of R = 0.40 R_⊙. Using these mass and radius values, we could estimate the brightness temperature of the emission.

The brightness temperature T_B of the emission is given by $$T_B={10}^{-26}\frac{c^2}{2k_B{\nu }^2}\frac{S}{\mathrm{\Omega }}\mathrm{K}={10}^{-26}\frac{c^2}{2k_B{\nu }^2}\frac{S{d}^2}{\pi R^2}\mathrm{K},$$(15)

where c is the speed of light in cm/s, k_B the Boltzmann constant in erg/K, v the frequency of the emission in Hz, S the flux density of the burst in mJy, Ω the source solid angle in sr, d the distance to the source in cm, and R the radius of the emitting region in cm. To obtain the second equality, we assumed that the emission region has a projected size comparable to that of the stellar disk. For our ${4.5}_{-1.3}^{+1.4}$ mJy burst, assuming that the star is not in a binary system, the brightness temperature is 6.1 × 10¹³ K. If the star is in a binary system with R = 0.40 R_⊙, the brightness temperature increases to 9.9 × 10¹³ K. Such high brightness temperatures at this frequency require a coherent radio emission process, which can only be produced by either plasma emission or through the electron cyclotron maser instability (ECMI; Treumann 2006).

The linear polarisation fraction of solar Type II bursts coinciding with CMEs is minimal, often approaching zero (Grognard & McLean 1973; Morosan et al. 2022; Dey et al. 2024). However, burst B20 does show weak linear polarisation features, as described by Callingham et al. (2025). For burst B21, we calculated the Faraday rotation measure and confirmed that our burst does not emit any significant linear polarisation components (see Appendix A).

Fig. 3 The 19 stellar radio bursts detected without drifting features. Each thumbnail shows a dynamic spectrum (Δν=0.4 MHz, Δt= the boxcar width). The interferometric images are shown as insets in the top left. The time series of the inverse variance weighted frequency-integrated spectra are shown at the top of each panel, with vertical dotted lines indicating the start and end times of the burst. The time-integrated spectra (shown to the right of each panel) are integrated over the time range indicated by the vertical dotted blue lines in the time series. Denoted in the top right of each thumbnail is the burst ID, which corresponds to Table 1. White areas in the dynamic spectra represent areas masked due to the presence of RFI, and the white ovals in the interferometric image represent the point spread function. For visual purposes, the colour maps of the dynamic spectra range from −20 to 20 mJy. Bursts B02, B04, B10, B11, B12, B14, B16, B17, and B20 are also presented in Tasse et al. (2025).

Fig. 4 Panel a: dynamic spectrum of B20 normalised by the blank-sky variances. The inset in the bottom left displays an interferometric image indicating that the emission originates from a point source. The horizontal white bands are frequency channels masked due to the presence of RFI. Panel b: dynamic spectrum of B20, with the highest likelihood drift rate (i.e. the overall drift rate) shown as a dashed red line. Panel c: frequency-integrated burst brightness as a function of the modelled drift delay. The intersection of the dotted red lines indicates the location of the highest S/N at ~1 minute; the horizontal line represents zero drift and corresponds to the original time series.

Fig. 5 Panel a: dynamic spectrum of B21 convolved with a kernel designed to follow the burst’s drift pattern. The drift of the burst is overplotted as a dotted red line. Panel b: raw dynamic spectrum of B21. The horizontal white bands are frequency channels masked due to the presence of RFI. The inset in the bottom left shows an interferometric image reconstructed along the dotted line in panel a, indicating that the emission originates from a point source. Panel c: frequency-integrated burst brightness as a function of the modelled drift delay. The dotted cross indicates the location of the highest S/N at a drift of 13 minutes, where the burst aligns vertically and thus appears offset from the horizontal zero drift line, which corresponds to the original time series.

4.3 Type II burst rate

With two detections of drifting bursts that exhibit characteristics consistent with Type II bursts, including their drift rates, emission frequencies, durations, and brightness, we can establish constraints on the Type II-like burst rate using two separate but complementary methods. The first is an estimate of the Type II burst rate based on observed detections, which constrains the burst rate by considering only the effective sample of stars from which bursts could be detected. This expands on the method presented by Callingham et al. (2025), applying it to each star rather than a single source. With the second approach, we parametrised the cumulative Type II-like burst luminosity distribution, f(L), as a power law and calculated for which parameters the likelihood is maximised. Both methods are visualised in Fig. 6 and are detailed below.

For our more observational approach, we calculated the burst rate for drifting Type II-like bursts with spectral energies (erg s⁻¹ Hz⁻¹) exceeding a given energy threshold (E) as $$R_{\mathrm{c}\mathrm{u}\mathrm{m}}(E)=\frac{N_{\mathrm{e}\mathrm{f}\mathrm{f}}(\ge E)}{T_{\mathrm{e}\mathrm{f}\mathrm{f}}(E)},$$(16)

where N_eff is the effective cumulative amount of detected bursts with ≥E from the effective sample of stars from which a burst with energy E could be detected, given our sensitivity of ~2.5 mJy for bursts of 1 minute. T_eff represents the total time spent observing this effective star sample. This estimation is slightly conservative as it does not include extremely bright bursts (e.g. ≥10¹⁸ erg s⁻¹ Hz⁻¹) from stars past 100 pc.

We show the Poissonian 3σ upper and lower limits for the burst rate per day as the grey monochromatic region in Fig. 6. Our two detections of drifting bursts significantly increase the lower limit of the burst rate, causing pronounced jumps when the energy surpasses the threshold at which the host stars of the drifting bursts can be observed. The first jump at lower energies is caused by burst B20, whereas the second higher jump is caused by both B20 and B21. The true luminosity distribution of these drifting radio bursts, regardless of its exact form, should lie within this grey monochromatic region. The approach by Callingham et al. (2025) can be used for bursts with radio luminosities similar to that of burst B20, whose Type II burst rate they determined to be ${0.84}_{-0.14}^{+1.10}$ × 10⁻³ per day. This is consistent with our upper limit for bursts with spectral energies of ~10¹⁶ erg s⁻¹ Hz⁻¹ or higher.

In our second approach, we parametrised the cumulative burst luminosity distribution, f(L), as a power law: $$f(L)={\left(\frac{L}{L_0}\right)}^{\alpha },$$(17)

where f(L) is the number of bursts with luminosity L exceeding some value per unit time, L₀ is the normalisation point where we expect on average one burst per unit time, and α is the power law index. We can think of each observation of a star as a single experiment (i), where we have a threshold luminosity $L_i^{\text{th\hspace{0.17em}}}=4\pi d_i^2{S}_i^{\text{th}}$ for detection. The total likelihood across all observations of all stars is given by $$\mathcal{L}(L_0,\alpha )=\prod _i\frac{\exp (-{\lambda }_i){\lambda }_i^{k_i}}{k_i!},$$(18)

where $${\lambda }_i=t_i{\left(\frac{4\pi d_i^2{S}_i^{\mathrm{t}\mathrm{h}}}{L_0}\right)}^{\alpha },$$(19)

where k_i is the number of detected drifting bursts from experiment i. We calculated the likelihood across a range of values, which yielded the two-dimensional posterior distribution for α and L₀, shown in Fig. 7. For the priors, we used a logarithmic prior for L₀, meaning that log(L₀) is uniformly distributed, while α is uniformly distributed. The main panel shows the joint posterior distribution, whereas the right and top panel shows the marginalised probability density functions (PDFs) for α and L₀, respectively. The maximum likelihood solution is attained at α = −0.732, and L₀ = 10^3.528 erg s⁻¹ Hz⁻¹, and is highlighted in Fig. 6 as a white dashed line. For α we calculated the 99.7% confidence interval as α = $-{0.7}_{-0.6}^{+0.9}$, L₀ is weakly constrained. Combinations of α and L₀ with a likelihood above 0.01 collectively form the diverging shaded region in the figure, where the intensity of the shading reflects the relative likelihood of each solution, with darker regions corresponding to higher likelihoods. As the likelihood decreases, the shading becomes more transparent, effectively conveying the gradual transition from the most likely solutions to less probable ones.

Both methods yield consistent estimates for the burst rate, with the more observational approach providing robust bounds and the parametric distribution offering insights into the underlying luminosity distribution. Additionally, we plot the cumulative spectral energy distribution of decametric solar Type II bursts from Mohan et al. (2024) in Fig. 6. The solar Type II data were observed by the Radio and Plasma Wave Investigation (WAVES) instruments on board the Wind, STEREO A, and STEREO B spacecrafts. We assumed that the number of Type II events collated byMohan et al. (2024) corresponds to all Type IIs observed during the quoted period from November 2006 to July 2023. We assumed Poissonian errors on the burst numbers and a fractional error of 20% for the reported peak flux values. We fit a power law to the distribution using the Python powerlaw package (Clauset et al. 2009). In addition to performing power law fits, this package identifies the initial turnover point in the cumulative distribution to determine the appropriate minimum spectral energy from which to start the fit. The fitted power law has a slope of −0.81 ± 0.06 ± 0.02, where the first error is determined as the quadrature sum of the error on the fit using the fractional 20% peak flux errors and the second error is determined by bootstrapping.

By restricting the total observation time to stars of a specific spectral type, we can refine our constraints based on the spectral type of the star. This is shown in Fig. 8. We assigned stars to spectral types based on their Gaia colours, using an approximate binning that is sufficient for our analysis. Stars were classified as M, K, G, F, A, or B types if their colours fell within the ranges 4.25–2.27, 2.27–1.08, 1.08–0.71, 0.71–0.33, 0.33 to −0.04, and −0.04 to −0.49, respectively. Using the more observational approach, we plot the upper limits for stars of spectral types for which no drifting bursts were detected. For M dwarfs, where we have two detections of drifting bursts, we again also parametrised and fitted the cumulative luminosity distribution, f(L), to estimate the burst rate. The maximum likelihood is achieved for α = −0.4, and L₀ = 10^2.1 erg s⁻¹ Hz⁻¹.

Fig. 6 Burst rate of circularly polarised drifting stellar radio bursts that exceed a given spectral energy threshold (in erg s−1 Hz−1). The monochromatic grey area represents the region spanned by the 3σ Poissonian limits; the dark dashed line indicates the upper limit. The two red data points represent detected drifting bursts, which cause sudden jumps in the Poisson lower limit when specific energies reach the detection threshold for that star. For energies exceeding that of a drifting burst, the lower limit drops again due to insufficient data on higher-energy bursts. The dashed white line represents the maximum likelihood of the cumulative luminosity distribution, while the blue diverging shaded region shows other possible solutions, with faded colours indicating a lower likelihood. The cumulative spectral energy distribution of decametric solar Type II radio bursts, to which a power law has been fitted using the Python powerlaw package (Clauset et al. 2009), is plotted in orange.

Fig. 7 Joint posterior probability distribution for the cumulative luminosity distribution parameters α and L0. The maximum likelihood is indicated with a plus-sign, and is reached when α = −0.7, and L0 = 103.5 erg s−1 Hz−1. The marginalised PDFs are shown on the top and right axes, respectively.

5 Discussion

In our ~100 stellar-year search, we detected 21 stellar radio bursts, 20 of which originated from M dwarfs. Given that M dwarfs are the most common star type, and exhibit high stellar activity and high magnetic field strengths, it is unsurprising that the majority of radio bursts were found in these stars. Only B12 originated from a RSCVn close binary. Among all 21 bursts, only B20 and B21 exhibit any measurable form of drift.

B20 and B21 are drifting coherent radio bursts, which can be produced either by plasma emission, where the burst would likely be tracing a CME, or by ECMI emission. Stellar ECMI bursts have been reported previously (e.g. Osten & Bastian 2008; Hallinan et al. 2008; Zarka et al. 2025), but they are typically different with respect to their durations and morphology, and are often modulated periodically with the host’s rotation. Assuming our two drifting bursts are produced by plasma emission, based on their specific drift structure, burst B21 could trace the backbone of a Type II burst, while the remainder of the flux of the burst has fallen below our sensitivity threshold. Both fundamental and harmonic emission could produce a burst with brightness temperatures of ~10¹⁴ K (Stepanov et al. 2001). However, harmonic emission cannot reach polarisation fractions >50% (Vedantham 2021). Burst B21 is significantly circularly polarised, 50–100%, and therefore if the burst is produced by plasma emission, it is likely fundamental emission. Burst B20 is likely produced by plasma emission, as argued by Callingham et al. (2025).

Assuming plasma emission, we calculated the corresponding shock speed of a burst based on its drift rate, enabling a direct comparison with solar CMEs. The shock speed must exceed both the stellar wind speed and the Alfvén speed for the shock to exist, and is therefore dependent on the physical properties of the stellar environment. We modelled the stellar wind using an isothermal Parker wind model with base density and temperature as parameters. The Parker wind model is used as a first-order approximation, since the complex M-dwarf wind structures are not well constrained (Lamers & Cassinelli 1999), and more sophisticated models would not qualitatively change our first-order comparison of CME speeds. We assumed the magnetic field to fall off as radial distance to the power of −3 (i.e. dipolar) out to a distance where the magnetic energy density equals the wind kinetic energy density. Beyond this radial distance we assumed that the field lines open up and fall of as distance to the power of −2. With the wind density and field specified, we can calculate the Alfvén speed at any radial distance. For both the host stars of bursts B20 and B21, we assumed a mass of 0.5 M_⊙ and a radius of 0.5 R_⊙, as variations in these parameters minimally affect the results. Additionally, we assumed a surface magnetic field of 100 Gauss. With various base wind densities and coronal temperatures, we determined the minimum radial speed of a shock required to exceed both the wind speed and the Alfvén speed. Next, by transforming the radial density profile to the emission plasma frequencies, we could calculate the maximum duration of a burst that could pass through our 120–168 MHz observing band. This is shown in Fig. 9a, where the white contour lines represent the durations of the drifting bursts B20 and B21. Solutions with shorter burst durations are also allowed as the CME shock speed can exceed the minimum speed required to surpass the wind and Alfvén speeds. The shock velocity ranges that could exactly produce burst B20 and burst B21 are shown in Figs. 9b and 9c, respectively. This demonstrates that a shockwave within an isothermal Parker wind model can produce the emission as seen in burst B20 with shock speeds upwards of ~1000 km s⁻¹, which agrees with the shock speed of 2400 km s⁻¹ quoted by Callingham et al. (2025). For burst B21, the emission can be produced by shock speeds ranging from ~300 km s⁻¹ upwards to a few thousand km s⁻¹, where the lower end is more likely as lower wind density scenarios are preferred. These speeds broadly agree with solar CME speeds, although burst B21 is notably on the slower side.

The ECMI can also produce coherent radio bursts, although this requires specific magnetic field configurations, making such scenarios quite rare (Zarka et al. 2025). The detection of identically drifting harmonic plasma emission would provide convincing evidence in favour of the CME scenario. Furthermore, observing CME signatures across multiple wavelengths would also amount of definitive evidence (e.g. Zucca et al. 2014; Ying et al. 2020; Namekata et al. 2024). However, given the extremely low burst rate (1.0 year⁻¹ with E > 6.8× 10¹³ erg s⁻¹ Hz⁻¹) and the challenges associated with coordinating multi-wavelength observations of the same source, this remains a difficult avenue to pursue. Nevertheless, with two detected drifting bursts, the burst rate presented in Fig. 6 can be interpreted, at best, as the occurrence rate of stellar CMEs, or, more conservatively, as the occurrence rate of drifting stellar radio bursts. This interpretation is particularly useful for guiding the survey design for future searches for CMEs.

The cumulative spectral energy distribution of solar Type II bursts has been overplotted, and the power law index from our cumulative luminosity distribution, α = $-{0.7}_{-0.6}^{+0.9}$, is in agreement with the solar Type II power law index α = −0.81 ± 0.06 ± 0.02. Notably, the cumulative spectral energy distribution of solar Type II bursts lies within our monotonic 3σ Poisson upper limits and slightly below our maximum likelihood model. This suggests that the low number of detected stellar Type II bursts is simply be due to a lack of sensitivity and time-on-sky.

The solar CME-flare relation is often extrapolated to active stars where CMEs might be more readily detectable (e.g. Moschou et al. 2019). However, this extrapolation remains largely unverified. We explored this by comparing our measured Type II burst rate to the M dwarf flare frequency distribution reported by Yang & Liu (2019). These authors parameterise their distribution as $$\frac{\mathrm{d}N}{\mathrm{d}E}=C{E}^{-\alpha },$$(20)

where C is a normalisation constant and where they determined α = 2.13. We extrapolated this power law and compared it with our observed cumulative Type II burst rate using the point where the variance on our burst rate is the lowest: a rate of one event per 10⁹ s, corresponding to energies ≥10¹⁶ erg s⁻¹ Hz⁻¹.

The cumulative flare rate for M dwarfs can be written as $$N(\ge E)={\int }_E^{\mathrm{\infty }}C{E}^{-2.13}dE$$(21) $$=\frac{C}{1.13}{E}^{-1.13},$$(22)

where we estimate C ≈ 7 × 10⁴⁰ from Fig. 3 of Yang & Liu (2019). Substituting our burst rate into this equation yields a flare energy threshold of E ≈ 2.6 × 10³⁷ ergs. This suggests that the rarity of our observed Type II radio bursts is consistent with stellar flares having energies ≳10³⁷ ergs, supporting a potential connection between high-energy flares and detectable radio emission. For context, the most energetic solar flares typically reach ~10³² ergs, meaning that these radio-associated events on M dwarfs would be orders of magnitude more powerful than the strongest solar flares. This could suggest that only the most energetic flares produce detectable radio emission with LOFAR.

Future observations, for example with the Square Kilometre Array (SKA), offer the potential to directly test the predictions of our maximum likelihood power law. SKA1-Low phase 1 (Macario et al. 2022) is estimated to have roughly an order of magnitude sensitivity increase compared to LOFAR, and an instantaneous field-of-view that is around half of that of LOFAR’s High Band Antennas (40 m versus 30 m diameter stations). These two taken together imply that in a survey with comparable telescope time-on-sky as ours the expected numbers of burst detection is 2 × 0.5 × 10^1.5 = 30.

The minimum detected burst energy at 100 pc for SKA1-LOW is 3 × 10¹⁵ erg s⁻¹ Hz⁻¹. At these energies, our derived 3σ Poissonian burst rate upper limit is around 5 × 10⁻⁸ s⁻¹, while our maximum likelihood power law indicates a burst rate of 1.5 × 10⁻⁹ s⁻¹. To estimate the observation time required for SKA to probe these regions, we assumed a stellar sample similar to this study, including all stars within 100 pc. This is a conservative estimate, as SKA’s greater sensitivity could detect bursts from more distant stars. Given this sample, the smaller field of view would result in ~50 visible stars per pointing. To probe our 3σ Poisson upper limit, SKA1-LOW would need to observe ~14 pointings for 8 hours each. Observing a further 450 pointings would probe the maximum likelihood power law regime, either leading to more detections or placing stronger constraints on the Type II burst rate.

Fig. 8 Burst rate of circularly polarised drifting stellar radio bursts that exceed a given spectral energy threshold (erg s−1 Hz−1) based on the spectral type of the star observed. As in Fig. 6, we plot a grey monochromatic region and diverging shaded region, but now only for M dwarfs, as two drifting bursts have been detected from this spectral type. The dashed white line again represents the maximum likelihood of the cumulative luminosity distribution only fitted to M dwarfs, while the blue diverging shaded region shows other possible solutions, with fading colour indicating a lower likelihood. The Poissonian upper limits for other spectral types observed are shown as dashed lines, with triangle markers emphasising their status as upper limits.

Fig. 9 Panel a: maximum burst duration in minutes through our observing band (120–168 MHz) as a function of coronal temperature and base wind density, modelled with an iso-thermal Parker wind model. White contour lines show the observed durations of bursts B20 and B21. Regions above the contours can produce bursts with the observed durations, with CME shock speeds exceeding the wind- and Alfvén speed. Panel b: CME shock speeds that exactly produce burst B20’s observed duration of 0.988 minutes, as a function of coronal temperature and base wind density, modelled with an isothermal Parker wind model. The black line marks the minimum shock speed exceeding both the Alfvén and wind speeds. Solutions exist for shock speeds upwards of ~1000 km s−1. Panel c: CME shock speeds that exactly produce burst B21’s observed duration of 13.334 minutes, as a function of coronal temperature and base wind density, modelled with an isothermal Parker wind model. The black line marks the minimum shock speed exceeding both the Alfvén and wind speeds. Solutions exist for shock speeds upwards of ~300 km s−1.

6 Conclusions and future work

We have conducted the largest unbiased search for stellar Type II radio bursts, targeting all stars within 100 pc in LoTSS; this resulted in a cumulative total of 107.21 stellar years of observation time. Our search for drifting stellar radio bursts led to the detection of two events: the previously published 2-minute burst from the M dwarf StKM1-1262 and a new 13-minute burst from the M dwarf LP215-56. These bursts exhibit drifts of $-{0.801}_{-0.003}^{+0.003}$ MHz s⁻¹ and $-{0.060}_{-0.002}^{+0.002}$ MHz s⁻¹, respectively, and both are significantly circularly polarised (>50%). Given their characteristics, we assumed that both bursts were produced by plasma emission that was likely tracing stellar CMEs. Additionally, we identified 19 non-drifting radio bursts, each originating from a different star, which warrant further investigation in future studies.

Using the two drifting burst detections, we calculated Poisson upper and lower limits for the burst rate of drifting stellar radio bursts, which can be interpreted as the occurrence rate of stellar CMEs. We also fitted our data to a cumulative luminosity distribution, obtaining a power law index (α) of $-{0.7}_{-0.6}^{+0.9}$. This result is consistent with decametric solar Type II observations, where α = −0.81 ± 0.06 ± 0.02. A natural next step is to compare our findings with LOFAR-detected metric solar radio bursts, which would provide a more robust comparison between stellar and solar Type II bursts.

We find that the rarity of drifting stellar radio bursts is comparable to that of high-energy stellar flares on M dwarfs (~10³⁵ ergs). However, this comparison is only based on event rates, and future work to establish a direct physical link between optical flares and Type II bursts will require substantial co-observing in the two bands. For instance, for the Alpha Centauri A star, SKA-Low phase-1 will be sensitive to bursts with spectral energy ≥10¹² erg s⁻¹ Hz⁻¹, which should occur once every 275 hours. Co-observing with a large survey that simultaneously observes multiple stars would also be fruitful. For example, an SKA survey similar to LoTSS, consisting of 500 pointings each lasting 8 hours, could effectively probe the maximum likelihood regime for energies ≥10¹⁵ erg s⁻¹ Hz⁻¹ and explore the burst rate regime derived in this study.

Data availability

The relevant code and data products for this work will be uploaded on Zenodo at the time of publication.

Acknowledgements

The LOFAR data in this manuscript were (partly) processed by the LOFAR Two-Metre Sky Survey (LoTSS) team. This team made use of the LOFAR direction independent calibration pipeline (https://github.com/lofar-astron/prefactor), which was deployed by the LOFAR e-infragroup on the Dutch National Grid infrastructure with support of the SURF Co-operative through grants e-infra 170194 e-infra 180169 (Mechev et al. 2017). The LoTSS direction dependent calibration and imaging pipeline (http://github.com/mhardcastle/ddf-pipeline/) was run on compute clusters at Leiden Observatory and the University of Hertfordshire, which are supported by a European Research Council Advanced Grant [NEWCLUSTERS-321271] and the UK Science and Technology Funding Council [ST/P000096/1]. The Jülich LOFAR Long Term Archive and the German LOFAR network are both coordinated and operated by the Jülich Supercomputing Centre (JSC), and computing resources on the supercomputer JUWELS at JSC were provided by the Gauss Centre for Supercomputing e.V. (grant CHTB00) through the John von Neumann Institute for Computing (NIC). DCK, HKV and EI acknowledge funding from the ERC starting grant ‘Stormchaser’ (grant number 101042416). HKV and SB acknowledge funding from the Dutch research council (NWO) under the talent programme (Vidi grant VI.Vidi.203.093). JRC acknowledges funding from the European Union via the European Research Council (ERC) grant Epaphus (project number: 101166008). MJH thanks STFC for support [ST/Y001249/1]. AD acknowledges support by the BMBF Verbundforschung under the grant 05A23STA. PZ acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 101020459 – Exoradio).

Appendix A Faraday rotation of the burst

We calculated the linear polarisation flux density across various rotation measure (RM) trials to determine whether burst B21 contains a linearly polarised emission component. The Faraday effect causes the plane of linear polarisation to rotate, with the degree of rotation proportional to the line-of-sight component of the magnetic field and the square of the wavelength. As a result, the linearly polarised signal oscillates sinusoidally as a function of frequency.

By combining the Stokes U and Stokes Q emission, we could compute the total linear polarisation flux density for each RM trial. By testing different RM values, we could find the polarised emission. As shown in Fig. A.1 for burst B21, no significant peak of excess emission is observed at any specific RM, indicating that the burst exhibits little to no linearly polarised emission.

Fig. A.1 Linearly polarised emission of burst B21 as a function of RM trials. The absence of a significant peak indicates little to no linear polarisation in the emission.

Appendix B Parker wind model

We modelled the stellar wind as an isothermal Parker wind model under the assumption of a steady, spherically symmetric outflow. We considered a star of mass M and radius R, and defined the wind properties using a mean molecular weight μ, an initial magnetic field strength B, a temperature T, and an initial number density n₀.

The wind velocity v is derived from the Parker wind solution, which satisfies the equation $${\left(\frac{v}{a}\right)}^2-\ln {\left(\frac{v}{a}\right)}^2=4\ln \left(\frac{r}{r_c}\right)+C,$$(B.1)

where a is the isothermal sound speed and r_c is the critical point given by $$r_c=\frac{GM}{2a^2}.$$(B.2)

We obtained the velocity profile, v(r), numerically and extract the wind speed at the stellar surface, v₀. The number density, n(r), follows from mass conservation: $$n(r)=n_0\frac{R^2{v}_0}{r^{2}v(r)}.$$(B.3)

We assumed a dipolar field that decays as r⁻³ until the Alfvén radius, where the ram pressure equals the magnetic pressure: $$B(r)=B_0{\left(\frac{r}{R}\right)}^{-3},\text{\hspace{0.17em}for\hspace{0.17em}}r<R_A,$$(B.4) $$B(r)=B\left(R_A\right){\left(\frac{r}{R_A}\right)}^{-2},\text{\hspace{0.17em}for\hspace{0.17em}}r\ge R_A.$$(B.5)

With this, we could calculate the Alfvén speed: $$v_A=\frac{B}{\sqrt{4\pi n{m}_p}}.$$(B.6)

The speed of the shock must exceed the wind speed in order to be a shock and must exceed the Alfvén speed to produce the instabilities required for radio emission. Using the density profile, we calculated the plasma frequency and determined the shock speed at the moment it passes through our LOFAR observing band (i.e. 120 MHz and 168 MHz) and at what distance r this happens. With this, we determined the drift delay of the burst: $$\mathrm{\Delta }t=\frac{r_2-r_1}{(v_{\text{shock\hspace{0.17em}},1}+v_{\text{shock\hspace{0.17em}},2})/2}.$$(B.7)

References

All Tables

All Figures

	Fig. 1 Colour-magnitude diagram of Gaia sources observed with LoTSS, applying the quality filters described in Sect. 2. The sources are categorised into two clusters based on their proximity in the colour–magnitude space.
In the text

Fig. 2 Panel a: mock Stokes V dynamic spectrum of a radio burst due to plasma emission modelled using the Parker Solar wind model with solar parameters, assuming a shock speed of 1000 km s−1 and a S/N of 30. Panel b: frequency-integrated burst brightness as a function of the drift delay and time, modelled using the Parker Solar wind model with solar parameters, assuming a shock speed of 1000 km s−1. The defining ‘bow-tie’-shaped feature of a broadband swept astrophysical signal is clearly visible. The red cross indicates the brightest pixel and its corresponding burst drift.

In the text

Fig. 3 The 19 stellar radio bursts detected without drifting features. Each thumbnail shows a dynamic spectrum (Δν=0.4 MHz, Δt= the boxcar width). The interferometric images are shown as insets in the top left. The time series of the inverse variance weighted frequency-integrated spectra are shown at the top of each panel, with vertical dotted lines indicating the start and end times of the burst. The time-integrated spectra (shown to the right of each panel) are integrated over the time range indicated by the vertical dotted blue lines in the time series. Denoted in the top right of each thumbnail is the burst ID, which corresponds to Table 1. White areas in the dynamic spectra represent areas masked due to the presence of RFI, and the white ovals in the interferometric image represent the point spread function. For visual purposes, the colour maps of the dynamic spectra range from −20 to 20 mJy. Bursts B02, B04, B10, B11, B12, B14, B16, B17, and B20 are also presented in Tasse et al. (2025).

In the text

Fig. 4 Panel a: dynamic spectrum of B20 normalised by the blank-sky variances. The inset in the bottom left displays an interferometric image indicating that the emission originates from a point source. The horizontal white bands are frequency channels masked due to the presence of RFI. Panel b: dynamic spectrum of B20, with the highest likelihood drift rate (i.e. the overall drift rate) shown as a dashed red line. Panel c: frequency-integrated burst brightness as a function of the modelled drift delay. The intersection of the dotted red lines indicates the location of the highest S/N at ~1 minute; the horizontal line represents zero drift and corresponds to the original time series.

In the text

Fig. 5 Panel a: dynamic spectrum of B21 convolved with a kernel designed to follow the burst’s drift pattern. The drift of the burst is overplotted as a dotted red line. Panel b: raw dynamic spectrum of B21. The horizontal white bands are frequency channels masked due to the presence of RFI. The inset in the bottom left shows an interferometric image reconstructed along the dotted line in panel a, indicating that the emission originates from a point source. Panel c: frequency-integrated burst brightness as a function of the modelled drift delay. The dotted cross indicates the location of the highest S/N at a drift of 13 minutes, where the burst aligns vertically and thus appears offset from the horizontal zero drift line, which corresponds to the original time series.

In the text

Fig. 6 Burst rate of circularly polarised drifting stellar radio bursts that exceed a given spectral energy threshold (in erg s−1 Hz−1). The monochromatic grey area represents the region spanned by the 3σ Poissonian limits; the dark dashed line indicates the upper limit. The two red data points represent detected drifting bursts, which cause sudden jumps in the Poisson lower limit when specific energies reach the detection threshold for that star. For energies exceeding that of a drifting burst, the lower limit drops again due to insufficient data on higher-energy bursts. The dashed white line represents the maximum likelihood of the cumulative luminosity distribution, while the blue diverging shaded region shows other possible solutions, with faded colours indicating a lower likelihood. The cumulative spectral energy distribution of decametric solar Type II radio bursts, to which a power law has been fitted using the Python powerlaw package (Clauset et al. 2009), is plotted in orange.

In the text

	Fig. 7 Joint posterior probability distribution for the cumulative luminosity distribution parameters α and L0. The maximum likelihood is indicated with a plus-sign, and is reached when α = −0.7, and L0 = 103.5 erg s−1 Hz−1. The marginalised PDFs are shown on the top and right axes, respectively.
In the text

Fig. 8 Burst rate of circularly polarised drifting stellar radio bursts that exceed a given spectral energy threshold (erg s−1 Hz−1) based on the spectral type of the star observed. As in Fig. 6, we plot a grey monochromatic region and diverging shaded region, but now only for M dwarfs, as two drifting bursts have been detected from this spectral type. The dashed white line again represents the maximum likelihood of the cumulative luminosity distribution only fitted to M dwarfs, while the blue diverging shaded region shows other possible solutions, with fading colour indicating a lower likelihood. The Poissonian upper limits for other spectral types observed are shown as dashed lines, with triangle markers emphasising their status as upper limits.

In the text

Fig. 9 Panel a: maximum burst duration in minutes through our observing band (120–168 MHz) as a function of coronal temperature and base wind density, modelled with an iso-thermal Parker wind model. White contour lines show the observed durations of bursts B20 and B21. Regions above the contours can produce bursts with the observed durations, with CME shock speeds exceeding the wind- and Alfvén speed. Panel b: CME shock speeds that exactly produce burst B20’s observed duration of 0.988 minutes, as a function of coronal temperature and base wind density, modelled with an isothermal Parker wind model. The black line marks the minimum shock speed exceeding both the Alfvén and wind speeds. Solutions exist for shock speeds upwards of ~1000 km s−1. Panel c: CME shock speeds that exactly produce burst B21’s observed duration of 13.334 minutes, as a function of coronal temperature and base wind density, modelled with an isothermal Parker wind model. The black line marks the minimum shock speed exceeding both the Alfvén and wind speeds. Solutions exist for shock speeds upwards of ~300 km s−1.

In the text

	Fig. A.1 Linearly polarised emission of burst B21 as a function of RM trials. The absence of a significant peak indicates little to no linear polarisation in the emission.
In the text