Coherent radio bursts from known M-dwarf planet-host YZ Ceti

Stellar characterization

Table 1 details the physical stellar properties for our target YZ Cet and comparison objects Prox Cen and GJ 1151, determined by combining multiple empirical relations to jointly constrain the mass, radius and bolometric luminosity (method in ref. ³³). The estimates employ a precise parallax from Gaia Data Release 2 (ref. ³⁴), with the effective temperature derived from combining the luminosity and radius. The bolometric flux measurements for YZ Cet and GJ 1151 are from ref. ³⁵, and for Prox Cen they are from ref. ³⁶. We further quote additional relevant activity properties, including the stellar rotation period (68.46 ± 1.00), from their respective references. There is a more precise V-band rotation period of 68.4 ± 0.05 measured²¹, but the formal statistical error may not encompass all forms of systematic error. To explore the impact of rotation period uncertainty over the ~90 d separation between epochs 2 and 5, we opted to use the less precise period when phase-wrapping with the synodic periods, finding that it broadens the uncertainty curves but does not impact our conclusions. We also note that optical flaring can be seen in each object’s Transiting Exoplanet Survey Satellite³⁷ light curves.

Radio data analysis

For our VLA observations of YZ Cet we used 3C147 as flux calibrator and J0116−2052 as gain calibrator, and calibrated the data in CASA³⁸ using the VLA pipeline. The VLA was in compact configuration: D in epochs 1–3, DnC in epoch 4 and C in epoch 5. We observed with the phase centre located halfway between YZ Cet and the 150 mJy nearby source PMN J0112−1658 (7.5 arcmin away from YZ Cet), to keep that source within the main lobe of the primary beam, then shifted the phase centre to the location of the star before imaging. We imaged the first three epochs together and each of the two follow-up epochs separately, using CASA’s ‘tclean’ task with W projection, multiscale imaging and multifrequency synthesis with three Taylor terms. We used natural weighting to maximize point source sensitivity.

For each dataset, we performed self-calibration on the target field using CASA’s ‘gaincal’ command: one or two rounds of phase-only self-calibration, followed by one or two rounds of amplitude and phase self-calibration. For amplitude self-calibration, we used gaincal’s solnorm=true parameter, normalizing gain amplitudes to an average value of one to avoid artificially increasing flux density. For example, in the self-calibrated image of epochs 1–3 (Supplementary Fig. 1), the peak flux density of PMN J0112−1658 remained consistent (148.4 mJy before self-calibration, 148.5 mJy after), while the uncertainty substantially improved: the root-mean-square error in a region near YZ Cet that is empty of bright sources was 120 μJy before, and 25 μJy after self-calibration.

After imaging, we masked the star out of the model so that the model contained only background sources, then subtracted the model from the visibility data to obtain residual visibilities containing only the star and noise. Supplementary Figs. 2 and 3 show images of radio bursts in epochs 2 and 5 in the background-subtracted data. In Stokes I, residual sidelobes are visible due to imperfect subtraction of PMN J0112−1658.

With the star at the phase centre, we used the ‘plotms’ task to average the residual visibilities across all baselines and frequencies, yielding a complex-valued time series. The real component is equivalent to the centre pixel of a natural-weighted image, yielding the flux density of the star. The imaginary component should not contain stellar flux, but exhibits comparable levels of noise due to thermal noise, radio-frequency interference (RFI), and sidelobes of imperfectly subtracted background sources. These residual background sidelobes can cause ‘ripples’ in the time series as the sidelobe pattern evolves over time. We calculated the standard deviation of the imaginary component to estimate the effective noise levels in the time series including these factors. We used the imaginary component of the time series to estimate noise levels because, in epochs without clearly detected stellar variability, the real and imaginary components have roughly similar standard deviations. For example, in epoch 1, the Stokes I standard deviation is 69 μJy (real) and 79 μJy (imaginary), and the Stokes V standard deviation is 50.5 μJy (both real and imaginary). The greater standard deviation for Stokes I than Stokes V illustrates the effect of residual sidelobes of background sources.

For 3 min integrations, we measured noise levels of 55–80 μJy in the Stokes I time series and 37–50 μJy in Stokes V. Without source confusion, the VLA’s theoretical sensitivity in 3 min is 22 μJy. As source confusion is not an issue in Stokes V, the elevated Stokes V noise levels are probably due to RFI, with data loss due to RFI flagging and low-level RFI in the remaining data. The Stokes I noise levels are affected by both RFI and imperfect background-source subtraction; both of these effects are enhanced by the VLA’s compact configuration.

Figure 1 shows the resulting time series as a function of orbital phase, where the shaded region shows ±3 times the estimated noise level on flux density in each of the 3 min time bins in that epoch. Supplementary Figs. 4 and 5 show detailed time series of epochs 2 and 5, the two epochs with coherent bursts. To identify bursts, we required a flux increase of >3σ during the burst compared with before or after the burst. For example, in epoch 2 (Supplementary Fig. 4), the events at 2.3 h, 5.1 h and 5–6 h satisfy this criterion, whereas a possible left-polarized event at 3.1 h constitutes only a 2σ flux enhancement.

For epochs without bursts, we measured or placed an upper limit on the quiescent emission levels using an image of the full epoch duration after background-source subtraction. For epoch 1, we measured an intensity of I = −39 μJy per beam in the image at the star’s location (star undetected) and an RMS in the image near the star’s location of σ = 25 μJy per beam, leading to a 3σ upper limit on source flux density of: S < 3σ = 75 μJy. In epoch 3, the star was detected with a peak flux density of 313 ± 20 μJy (Stokes I) and 18.7 ± 5.3 μJy (Stokes V). In epoch 4, the star was undetected with an intensity of 36 ± 21 μJy at its location in the image, yielding a 3σ upper limit of 64 μJy on flux density. Weakly polarized, non-thermal, slowly varying quiescent emission on M dwarfs, such as in epoch 3, is typically attributed to incoherent gyrosynchrotron emission³⁹.

The degree of circular polarization of a signal is r_c = V/I = (RR − LL)/(RR + LL), where RR and LL are the right circular polarization and left circular polarization visibility data, respectively. To assess the bursts’ degree of circular polarization, we used a maximum likelihood approach to estimate r_c and construct a 68% confidence interval. To calculate the likelihood, we assumed RR and LL are Gaussian distributed, obtaining the standard deviation for each from the imaginary component of the time series. We generated a probability distribution function (PDF) for obtaining the data in terms of model parameters S_I (Stokes I flux density) and r_c, then marginalized the distribution across S_I to obtain a PDF for r_c alone. The black line in Supplementary Figs. 4 and 5 shows the value of r_c at which the PDF peaks, and the grey confidence interval shows the range of r_c that lies from 0.16 to 0.84 in the cumulative distribution function.

We produced dynamic spectra (Fig. 2 and Supplementary Fig. 6) for all of epoch 2 and for a short time period surrounding the epoch 5 burst, using the baseline-averaging code described in ref. ¹⁷. The flux density uncertainties quoted in the caption are calculated using the imaginary component of the dynamic spectrum (which does not contain stellar emission), by taking the standard deviation in each frequency channel then calculating the median across all channels. The peak flux densities in Stokes V for the epoch 2 bursts at 2.3 h and 5.1 h both exceed 5σ, as does the epoch 5 burst. The incoherent flare in epoch 2 has weak right polarization, appearing only faintly in Stokes V except for the coincident LCP burst at 5.1 h, and in Stokes I it spans the entire 2–3.7 GHz band, consistent with the broadband nature of gyrosynchrotron emission.

In epoch 2, the right-polarized burst and the left-polarized burst at 5.1 h both are brightest at lowest frequencies (<3 GHz). Like the two clearest features in the epoch 2 dynamic spectrum, the epoch 5 event is also brightest at the lowest frequencies. These three events that appear clearly in the dynamic spectra drop off above 2.5–3 GHz. If the emission process is the cyclotron maser, this indicates that the maximum magnetic field in the source regions is of the order of 1 kG.

Stellar magnetospheric environment

To determine whether SPIs could have powered our observed polarized radio emissions we needed to characterize the likely magnetospheric environment impacting the YZ Cet planetary system. We considered two models: (1) a magnetosphere defined by a radial isothermal stellar wind whose properties are set by the corona and the surface magnetic field strength, and (2) a PFSS extrapolation of typical M-dwarf ZDI measurements including an isothermal stellar wind solution beyond the source surface⁴⁰. As the magnetic field and wind environments of low-mass stars are very uncertain, this approach explores the effects of a range of likely stellar magnetic field strengths experienced by the YZ Cet planets.

The first approach, often employed in the literature, uses a stellar wind originating from the stellar surface, which overestimates the magnetic field at the planet location because it does not take into account the rapid radial decay of closed field lines near the stellar surface (for example, ref. ⁷). The second approach accounts for this effect by using a more realistic stellar magnetic field topology (for example, ref. ⁴⁰); however, the inherent assumptions exclude additional stresses to the magnetic field, and may underestimate the strength of the magnetic field at planetary distances from the star, beyond a poorly constrained source surface.

Model A (radial stellar field, high \(\dot{M}\))

To formulate the radial stellar wind solution we used a Weber wind model⁴¹, which solves the ideal magnetohydrodynamics problem in spherical coordinates for an axisymmetric equatorial wind propelled by a rotating star. Note that there is a typo in their equation (23). In the last term of the denominator in brackets, the factor \({{{\varOmega }}}^{2}{r}^{2}{M}_{{\mathrm{A}}}^{2}\) should be \({{{\varOmega }}}^{2}{r}^{2}{M}_{{\mathrm{A}}}^{4}\), for angular rotation rate Ω, radial coordinate r, and radial Alfvénic Mach number M_A. Their solution is similar to that of ref. ⁴², but self-consistently incorporates the stresses of the ionized wind on the magnetic field anchored to the rotating star. We depart from ref. ⁴¹, by employing an isothermal wind, a subcase of their general polytropic approach. The physical wind solution is the one that smoothly passes through the three critical points (one sonic, two magnetic)⁴¹, constraining the solution and fixing the initial radial wind velocity for our choice of boundary conditions. The coupled solution requires as inputs: stellar mass, radius and rotation rate (Table 1), as well as coronal plasma temperature, an average mass-loss rate and the radial magnetic field strength at the surface (see below). With these assumptions, we numerically solved for the wind radial profile. In practice, the system of equations incorporating the critical points fixes the total energy, a constant of the motion. We then employed the energy equation to numerically solve for the radial speed as a function of distance from the star. The other system properties could then be determined from the radial wind profile⁴¹. In summary, the important variable parameters for the boundary conditions reduce to the coronal temperature, radial magnetic field strength and mass-loss rate.

For this model (A), we assumed a coronal temperature of k_BT = 0.25 keV ≈ 3 × 10⁶ K, a constant mass-loss rate of \(\dot{M}\equiv 4\uppi \rho u{r}^{2}=1{0}^{-13}\,{{M}}_{\odot }\,{\mathrm{yr}}^{-1}= 5\,{\dot{M}}_{\odot }\) (fives times the mass-loss rate of the Sun), and radial field of B_r = 220 G. The assumed constant mass-loss rate sets the relationship between the mass density profile, ρ, and the radial wind speed, u. The coronal temperature is similar to that of other inactive late M dwarfs, based on their X-ray observations (for example, refs. ^43,44). The mass-loss rate is a compromise between the expected rates for similar stars based on their Rossby numbers (0.5 for YZ Cet)⁴⁰, and the low rate of Prox Cen (see within ref. ¹⁵), which has similar physical properties, albeit with a weaker magnetic field strength (Table 1).

For the estimate of the average radial surface magnetic field strength of YZ Cet, we utilized the measured large-scale field topology of Prox Cen, from ref. ⁴⁵, as no such measurements of YZ Cet are yet available. Because of their similar properties (Table 1), Prox Cen is a useful analogue for interpreting the magnetic properties of YZ Cet, and is one of the few slowly rotating late M dwarfs with a measured field topology from ZDI. We scaled the measured magnetic field of Prox Cen, in its spherical harmonic decomposition, based on the measured Zeeman broadening measurement of YZ Cet, <B_ZB> = 2,200 G (ref. ⁴⁶). We defined the scaling to achieve an average field flux strength ratio of ζ ≡ < B_ZDI>/< B_ZB> ≈ 0.1, and thus an average radial surface field of 220 G. We chose ζ = 0.1 as a representative value for the sample of stars with similar properties in ref. ⁴⁷ that have both kinds of Zeeman measurements. The low value of ζ originates from field cancellation in the Stokes V ZDI measurements, as opposed to the Stokes I ZB measurements that include the total field strength. It is worth noting that the YZ Cet Zeeman broadening measurement is a high outlier for its Rossby number of ~0.5 (ref. ⁴⁸), and the source measurements⁴⁶ may be systematically high⁴⁸, especially for slow rotators. However, ZDI has measured M-dwarf ζ values up to ~0.3 (Prox Cen⁴⁵), so our estimated average large-scale field of 220 G may be reasonable even if the current YZ Cet ZB measurement is an overestimate. With an average large-scale field of 220 G, surface variations and small-scale fields in the low stellar corona could still lead to regions with kilogauss field strengths, plausibly allowing ECM emission at 2–3 GHz.

Model B (PFSS stellar field, low \(\dot{M}\))

For model B, we used the same assumptions for the wind properties with two key differences. The first is that we moved the inner boundary of the Weber wind model, where the stellar field is purely radial, to 4.5 stellar radii, consistent with magnetohydrodynamic simulations of M-dwarf winds (Table 2 in ref. ⁴⁹). Shifting this ‘source surface’ outwards accounts for closed magnetic field lines near the surface. We modelled this closed field by filling the space between the stellar surface and the wind source surface with a PFSS extrapolation (for example, ref. ⁵⁰) based on the field topology of Prox Cen, scaled to yield an average large-scale radial field of 220 G at the stellar surface. The PFSS extrapolation thus sets the average radial magnetic field strength at 4.5 stellar radii from the star. Second, we also assumed a 20-times lower mass-loss rate of 0.25 \({\dot{M}}_{\odot }\), comparable to the upper limit on Prox Cen wind from ref. ⁵¹ and consistent with the predicted rate from ref. ¹⁵. This value is also close to the expectation (~0.23 \({\dot{M}}_{\odot }\)) calculated from the relation between X-ray surface flux and mass-loss rate^52,53. We compare these distinct magnetic environments in Supplementary Fig. 7.

We further illustrate some properties of model B in Supplementary Fig. 8, showing the velocities relevant to the wind (top panel), and the total wind pressure throughout the model environment around YZ Cet. We considered model B to correspond to our most realistic estimate of the average magnetized environment pervading this planetary system, whereas model A encapsulates typical assumptions in the literature treatment of these questions. While informed by the literature, the wind parameters are typically uncertain for low-mass stars, but as we employed an analytic model, we can readily change the input assumptions to determine their effect on the potential for the YZ Cet planetary system to power radio emissions (see below). To provide some intuition for the isothermal wind solution and the impacts of these parameter assumptions, we note that changing the temperature is the most impactful parameter determining the wind velocity, changes in the mass-loss rate largely impact the wind density, and the radial field strength scales the overall magnetic field as the azimuthal field component is much weaker for slowly rotating systems. In the absence of a three-dimensional wind simulation (for example, ref. ¹⁵), these simplified isothermal approaches provide a reasonable means to examine the approximate interplanetary environment conditions⁵³.

Planet-induced radio emission

Our detection of polarized radio bursts from YZ Cet prompts the question of whether the coherent radio emission could have been powered by the magnetic interaction of the star with its planets (see within ref. ¹). We used models A and B (described above), to define the magnetized stellar wind filling the environment of the YZ Cet planetary system. When this wind interacts with the planets, the dissipated energy can power auroral radio emissions. We estimated the available power through this interaction using the frameworks of ref. ²⁵ (reconnection) and ref. ⁶ (Alfvén wings), similar to the approach taken by ref. ⁹.

The available power released by magnetic reconnection²⁵ is

in cgs units, where γ is a geometric factor, R_o is the radius of the obstacle, that is, the planetary magnetosphere, υ is the velocity of interaction in the frame of the planet, and B is the star’s magnetic field strength at the planet location. Similarly, the available power transmitted through Alfvén wings⁶, a prediction valid in the low Mach number regime, is

as expressed by ref. ³, where \(\bar{\alpha }\) is an interaction strength, θ indicates the angle between the wind’s relative velocity vector and the magnetic field in the planet’s frame of reference, and ρ is the mass density of the magnetized flow. These two approaches differ by a factor of twice the Alfvénic Mach number⁶ and a geometric factor. We take the wind properties for models A and B, and use these expressions to estimate the expected power available to generate ECM radio emissions. At the planet’s location, the wind and magnetic field are aligned and nearly radial due to the star’s slow rotation, but the planet’s orbital velocity (small relative to wind speed) gives θ a small, non-zero value. We focused on YZ Cet b as the closest-in planet and most likely to power the detected hundreds of microjansky bursts in our radio datasets.

In evaluating equations (1) and (2), we take the middle value of the geometric factor, so γ → 1/2, and consider the interaction strength \(\bar{\alpha }\to 1\) . The former is justified by our ignorance of the exact geometry of the interacting magnetic fields²⁵. For the latter, we justify the assumption of the interaction strength based on the likely conductivity of the planetary obstacle through its magnetosphere or ionosphere, considering the environments of the large close-in rocky planets of the YZ Cet system to have high Pederson conductivities (see Appendix A of ref. ⁷).

The last remaining variable in the power expressions is the planetary obstacle radius, R_o. This is defined by the size of the planetary magnetosphere, or at a minimum the radius of the planet itself assuming a thin ionosphere. We use the pressure balance between the supposed planetary field and the wind to define the radius of the planetary magnetopause:

where B_p is the assumed planetary dipole field strength, μ = 0.5 for a fully ionized hydrogen wind and m_p is the proton mass. If the ratio of R_o/R_p from equation (3) falls below unity, we instead use R_p as the obstacle radius.

The YZ Cet system was characterized with radial velocity measurements and does not exhibit transits, so the planet radii are unknown. The planets are likely to be roughly Earth-sized, and YZ Cet b has a minimum mass of 0.7 M_⊕. For the radii of YZ Cet b we consider a range from R_p = 0.89 R_⊕ to R_p = 1 R_⊕, where the lower bound corresponds to the radius of the minimum mass assuming it also has an Earth-like density. As the planet is roughly Earth-sized, we explore a range of planetary dipole field strengths starting from 1 G (Earth-like), increasing it by an order of magnitude (10 G), and decreasing it to below the stellar field strength at the planet location (effectively unmagnetized). These values set the abscissa range in Fig. 3.

With these assumptions, we can compute the energy available to power auroral radio bursts from YZ Cet, using both the reconnection and Alfvén wing prescriptions, as well as considering both the model A and model B wind environments. To convert the power to a possible burst radio flux density we use

where S comes from equation (1) or equation (2), ϵ = 0.01 is the radio efficiency factor^1,6, Δν = 3 GHz is the emission bandwidth, for which we assume the emission spans from low frequencies to our emission band, d = 3.712 pc is the distance to the star, and we use Ω = 0.16 sr for the beaming angle based on the observed value for Jupiter–Io radio emission⁵⁴. The results of our calculations are shown in Fig. 3, and discussed in the main text.

SPI parameter space

The predicted flux densities for SPI depend on a variety of unknown properties for the magnetized environment, most prominently the assumed stellar mass-loss rate, and the stellar field strength. Above, we chose models A and B to represent a range of values consistent with the literature and the known physical properties of the star. Below, we explore two specific effects related to these assumptions: the range of mass-loss rates consistent with the SPI scenario and the dependence of the SPI power on the assumed stellar magnetic field.

Constraints on mass-loss rate

If our detected bursts are indeed powered by sub-Alfvénic SPI, then the corresponding planet must be within the Alfvén surface of the stellar environment. Using our isothermal wind solution, we explored the impact of the assumed mass-loss rate on the Alfvénic Mach number at the position of the planets around YZ Cet. We show these results in Supplementary Fig. 9. Both fiducial models A and B allow for an increase of an order of magnitude in the assumed mass-loss rate before any planets go beyond the Alfvén surface, and even more before YZ Cet b becomes super-Alfvénic.

Because the star is rotating slowly, the wind speed and stellar field are largely radial, and as the orbital motion of the planets is small compared with the wind speed, the results in Supplementary Fig. 9 are well approximated by

where the velocity and field on the right-hand side correspond to the radial components evaluated at the position of the planets. For our fiducial model A (5 \({\dot{M}}_{\odot }\)), the planets a, b and c, become super-Alfvénic at mass-loss rates of approximately 150, 80 and 50 times \({\dot{M}}_{\odot }\), respectively. For our fiducial model B (0.25 \({\dot{M}}_{\odot }\)), the planets a, b and c, become super-Alfvénic at mass-loss rates of approximately 13.5, 7.5 and 4 times \({\dot{M}}_{\odot }\), respectively. Increasing the assumed stellar magnetic field would increase the distance corresponding to the Alfvén surface, and increase the planets’ corresponding sub/super-Alfvénic transition mass-loss rate. If our radio detections are powered by the interaction of YZ Cet b with its host, it should imply that YZ Cet has a mass-loss rate within these bounds, probably <13.5 \({\dot{M}}_{\odot }\), using our more realistic model B stellar field topology.

Scaling the stellar magnetic field

As discussed when introducing models A and B above, we conservatively assumed a surface radial average field strength of 220 G, but this may be an underestimate. For both our models A and B, scaling the magnetic field at the surface towards higher values linearly scales the field at the location of the planet, and similarly scales the minimum predicted SPI flux density as a function of planet field strength (flat regions in Fig. 3). The turning points in Fig. 3 also shift towards higher planet field strengths, as their position encodes where the stellar and planetary fields balance.

For model A, a stellar field increase of a factor of 2 or 3 places flux density predictions above the measured bursts, for both the reconnection and Alfvén wing mechanisms. If our radio detections are indeed SPIs, this underscores that model A overestimates the stellar field strength at the planet’s location by not accounting for closed field structures. Literature SPI predictions with stellar radial fields, like what we have assumed with model A, may overpredict both SPI intensities and the size of the Alfvén surface and which planets it encompasses.

For model B, an increase in the field strength of a factor of 2 or 3 pushes both the reconnection and Alfvén wing predictions to higher values. While the reconnection prediction would become more discrepant with the measured burst flux densities, it would place the Alfvén wing prediction consistent with smaller (but non-negligible) planetary dipole fields. If our detections are indeed SPIs, and if additional measurements reveal that the average global field of YZ Cet substantially exceeds 220 G, then our result would support model B (accounting for closed field near the surface) over model A (open stellar field), and the Alfvén wing mechanism over the reconnection mechanism.