© The Authors 2025

1 Introduction

Brown dwarfs are known to exhibit rotationally modulated coherent radio bursts (Hallinan et al. 2006, 2007, 2008; Berger et al. 2009; Route & Wolszczan 2012; Williams & Berger 2015; Route & Wolszczan 2016; Kao et al. 2016, 2018; Vedantham et al. 2023; Rose et al. 2023). These radio bursts are generally attributed to electron cyclotron maser instability (ECMI; Wu & Lee 1979; Treumann 2006; Zarka 2007; Hallinan et al. 2008) from non-thermal charges trapped in the magnetosphere, which is the principal mechanism that drives the auroral¹ radio emission in Jupiter and other magnetised planets in the Solar System (Zarka 1998).

Owing to remote and in situ observations of Jovian auroral decametric emission (DAM; see Marques et al. 2017 and references therein), the basic properties of electron cyclotron maser (ECM) radio emission stemming from a magnetopshere are well understood: ECM emission occurs at the local cyclotron frequency and is beamed along the surface of a cone. The opening angle of the cone is wide and its axis parallel to the magnetic field line. Therefore, the aggregate radio emission from different magnetospheric regions typically appears as a broadband signal that is modulated by the rotation rate of the magnetosphere (akin to a lighthouse; Callingham et al. 2024). In the case of Jupiter, its ECM emission is modulated at around 10 hours, which is the rotation period of the Jovian magnetosphere. This measurement technique has provided the most accurate measurement of the rotation period of gas giants in the Solar System (e.g. Higgins et al. 1997; Carbary & Mitchell 2013). It has also been successful in determining the rotation period of brown dwarfs down to spectral class T (i.e. methane dwarfs; Callingham et al. 2024).

Jovian ECM emission is also modulated in a second period corresponding to the orbital period of Io (Marques et al. 2017). This component of the emission is induced by the interaction between the conductive ionosphere of Io and the Jovian magnetic field (Zarka 1998). Because the electrons involved in this component only exist along the flux tube connecting Jupiter and Io, beaming ensures that the emission reaches the observer only when Io is at certain longitudes in its orbit. Identifying a second periodicity in brown dwarfs could analogously lead to the identification of exoplanets. Exoplanets around cold methane dwarfs likely exist but have eluded discovery so far (Kipping 2009; Kipping et al. 2012).

The rotation periods of radio-loud brown dwarfs are usually in the range of a few to several hours measured via their radio light curves (Kao et al. 2016, 2018; Williams & Berger 2015; Route & Wolszczan 2012, 2016). The rotational periodic modulation can therefore be detected with a modest amount of telescope exposure. However, the orbital periods of putative exoplanets around brown dwarfs could be in the range of several days. Detecting the second periodicity therefore requires significant telescope time.

Vedantham et al. (2023) recently discovered a new radio-loud brown dwarf system. Located at 23.3 ± 1.0 pc (Kirkpatrick et al. 2019), WISEP J101905.63+652954.2 (= 2MASS J10190575+6529526 = WISEA J101905.61+652954.0; hereinafter J1019+65) consists of two brown dwarfs in a binary of spectral types T5.5 ± 0.5 and T7.0 ± 0.5 (Vedantham et al. 2023). The two brown dwarfs have a projected separation of 423.0 ± 1.6 mas, which corresponds to 9.9 ± 0.4 AU (Vedantham et al. 2023). J1019+65 displayed circularly polarised radio pulsations at 144 MHz, based on data from the LOFAR Two-metre Sky Survey (LoTSS; Shimwell et al. 2017). This marks the first methane dwarf binary to be detected in the radio wavelength. Using the Lomb-Scargle (LS) periodogram (Lomb 1976; Scargle 1982), Vedantham et al. (2023) also revealed a periodic signature of about three hours from the LOFAR-band light curve of J1019+65. Vedantham et al. (2023) interpreted this radio-derived three-hour period as the rotation period of the radio emitter in J1019+65. However, since LoTSS has an absolute astrometric precision of 0.5″ (defined by the Gaussian-equivalent standard deviation; Shimwell et al. 2022), it is not possible to discern which of the two brown dwarfs in J1019+65 is the radio emitter (also noting that possibly both could be) with the LoTSS data alone due to insufficient angular resolution.

Despite the complication of identifying the radio-loud component in J1019+65, it is a well-suited system for monitoring for the following reasons. Since LoTSS is an untargeted allsky survey, the radio detection of J1019+65 does not suffer from selection bias, which is inherent to targeted observations. Thus, by studying J1019+65 further, we can learn more about the population of radio-loud brown dwarfs using J1019+65 as an archetype. Furthermore, unlike isolated brown dwarfs, the binary nature of J1019+65 will allow for an independent mass constraint using future IR observations (Vedantham et al. 2023).

In this work, we present our observational campaign of J1019+65 with the Low-Frequency Array (LOFAR). The campaign comprises 44 hours of new radio observations plus the original eight-hour LoTSS data, spanning six years (2017-2023), making it one of the most intensive long-term monitoring of a brown dwarf system in radio wavelength to date. The paper is structured as follows. In Sect. 2, we present the details regarding the follow-up observations of J1019+65 carried out by LOFAR, which also include the description of the data reduction procedure. We discuss our results in Sect. 3 and conclude in Sect. 4.

2 Observations and data reduction

Since the discovery of J1019+65's radio emission in the LOFAR Two-Metre Sky Survey data release 2 (LoTSS DR2) exposure, we have acquired ten new follow-up observations with LOFAR as part of projects LC15_019 & LT16_013 (details in Table 1). Except for the first epoch (2021) of our follow-up observation (project code: LC15_019), which has the same exposure time (8 hours) as the LoTSS discovery field, all other LOFAR observations (project code: LT16_013) lasted four hours each. This ensured that we were able to covere the full three-hour rotation (with a 2σ uncertainty range of 2.6-3.8 hours; Vedantham et al. 2023) of J1019+65 in each observation.

We processed the LOFAR dataset for each epoch using the standard LOFAR Initial Calibration (LINC) pipeline² to flag interference and correct for direction-independent errors. We refer to de Gasperin et al. (2019) for more details of these effects and the calibration strategy used in the LINC pipeline. Afterwards, we processed each dataset using the direction-dependent self-calibration and imaging pipeline (DDF-pipeline³). The DDF-pipeline is described by Shimwell et al. (2022) and Tasse et al. (2021). The procedure also involved the extraction and self-calibration to remove other sources in the field before selfcalibrating on the target field, as described by van Weeren et al. (2021). Finally, we synthesised the full-exposure Stokes I and V images of J1019+65 for each LOFAR epoch using wsclean with Briggs weighting (Briggs 1995). We used robustness parameters of -0.5 and 0.0 for Stokes I and Stokes V images, respectively.

The post-processing analysis for all datasets follows the procedure described by Vedantham et al. (2020, 2023). For each calibrated dataset, we first used wsclean (Offringa et al. 2014) to create Stokes V dirty snapshot images synthesised over the full bandwidth (120-168 MHz) at a cadence of four minutes (4 min) after phase-shifting the visibilities to the proper-motion-corrected sky location of J1019+65. We used the ephemeris of Kirkpatrick et al. (2019) to account for proper motion. Furthermore, we split the data into three spectral bins (120-136 MHz, 136-152 MHz, and 152-168 MHz) to check for any spectral evolution in the light curves. We then produced light curves for each epoch and spectral bin by extracting the value of the flux density at J1019+65 from the 4 min images. For each light curve, we computed the LS periodogram using the astropy implementation (Astropy Collaboration 2013, 2018, 2022) to determine whether there were any strong peaks in the LS power spectrum of J1019+65. To quantify the significance of a periodogram peak, we computed its false-alarm probability (FAP) using the approach of Baluev (2008) that was included in the astropy implementation. The FAP represents the probability that a light curve with no periodic component would lead to a peak of a similar height as the one detected (VanderPlas 2018). Figure A.3 in the appendix shows the post-processed light curve and periodogram from the original LoTSS 2017 epoch and serves as a reference for our follow-up radio observations.

2.1 Single-epoch analysis

We individually inspected all the Stokes V light curves of the radio observations. None of the ten LOFAR follow-up epochs shows radio bursts in four-minute cadence data above a threshold of 5σ. The lack of new short bursts in the follow-up observations is surprising, but not unexpected because ECMI, a coherent emission mechanism, is intrinsically highly variable. For example, the occurrence probability of Jovian DAM emission (as a function of central meridian longitude in Jupiter’s System III coordinates or CML for short) only reaches a maximum of around 30% at certain values of CML (Marques et al. 2017). As the occurrence probability of ECM emission from known brown dwarf emitters remains largely unconstrained, it is entirely possible that J1019+65 was simply not producing sufficiently bright and singular ECM bursts to be adequately detected during most of our observations.

In terms of periodicity, only the eight hour (8h) follow-up observation (2021-02-23) shows a significant detection (FAP ≪ 0.05). The rest are only four hours long and cannot reliably measure the 3 h rotation period of J1019+65 in a single exposure. Figure 1 shows the radio light curve of the 2021-02-23 epoch and the corresponding LS periodogram. Although we did not detect a single radio pulsation that is as strong as the original pulse from LoTSS DR2, the periodogram clearly detects the original period of ≈3 h. We hereafter refer to this frequency as f₁, as represented in Fig. 1. However, in contrast to the original 2017-03-08 epoch, the f₁ peak is not significant for the full-bandwidth light curve (blue); the f₁ peak is much more significant (FAP ≪ 0.003) at 120-136 MHz (orange), while the other two spectral bins (136-152 MHz and 152-168 MHz; green and purple) did not contribute to this peak. This suggests that there is some spectral shape evolution with time in the rotationally modulated emission. The other nine LOFAR single-epoch light curves and their corresponding LS periodograms can be found in Figs. A.1 and A.2 in the appendix.

Fig. 1 LOFAR observation of J1019+65 carried out on 2021-02-33. Top panel: Stokes V light curves with a cadence of 4 min at different frequencies. The shaded region represents ±1σ uncertainty. The light curve at 152-168 MHz is omitted for the sake of clarity. The top axis represents the time of observation in Modified Julian Dates (MJD). The dashed black line at 0 mJy is drawn for clarity. Bottom panel: LS periodogram of the radio light curves. The three dashed red lines represent the necessary LS power (i.e. peak height) to achieve a false-alarm probability (FAP) of 5%, 1%, and 0.3%. The grey curve represents the (negative) LS periodogram of the window function, which is a light curve with the same timestamps as the original curve, but the flux density of each data point is replaced with unity (i.e. a flat light curve). The two relevant frequencies, f1 = 0.315 ± 0.013 h−1 and f2 = 1.247 ± 0.013 h−1, are indicated by the black arrows in the LS periodogram. The 120-136 MHz (orange) peak at f1 is consistent with the original frequency discovered by Vedantham et al. (2023), while the relevancy of the f2 is discussed in Sect. 2.2.2. The frequency uncertainties were computed using Eq. (52) by VanderPlas (2018).

Fig. 2 LS periodograms of the LOFAR cross-epoch radio light curve, which includes data from both the original LoTSS DR2 epoch and our ten LOFAR follow-up epochs. The spectrum reveals the peak at f1 ≈ 0.325 h−1 originally discovered by Vedantham et al. (2023). Another significant peak not previously discovered is detected at f2 ≈ 1.271 h−1, with a false alarm probability (FAP) of less than 0.3%. The locations of these two peaks are indicated by the black arrows. The three dashed red lines represent the FAP of different values. The grey curve represents the LS periodogram of the window function, which is a light curve with the same timestamps as the original curve; however, the flux density of each data point is replaced by unity (i.e. a flat light curve). The inset plot in each panel shows a zoomed-in version of the observed J1019+65 power spectrum in the vicinity of the tallest LS peak, with the same three FAP dashed red lines.

2.2 LOFAR cross-epoch analysis

We also concatenated all ten LOFAR follow-up light curves plus the original LoTSS DR2 light curve of J1019+65 to create a cross-epoch light curve that spans six years (from 2017-03-08 to 2023-03-20). We then computed its LS periodogram, as shown in Fig. 2. We note that it is possible to see the peculiar fine-scale structure around the spectral peaks in the zoomed-in inset figure. The structure stems from the observational pattern (i.e. the window function) because the observed spectrum is a convolution of the true spectrum and the window spectrum (grey), which has a complex structure due to a combination of long temporal baseline (>2 year) of our data and the 4 min ‘sampling’ interval of the light curve. We refer to VanderPlas (2018) for more details regarding the various effects of the window function on the observed power spectrum in the LS periodogram.

2.2.1 New periodic signature

Figure 2 shows that the original periodogram peak originally discovered by Vedantham et al. (2023) is present at f₁ ≈ 0.325 h⁻¹, as expected. In addition, we discovered a new periodicity in our cross-epoch light curve that was not seen by Vedantham et al. (2023): there exists a statistically significant peak at f₂ ≈ 1.271 h⁻¹, which corresponds to a period of P₂ ≈ 0.787 h. The FAP of this peak is less than 0.3%.

We note that we do not report the formal uncertainty of the measured frequency in this periodogram analysis due to the reasons discussed by VanderPlas (2018). Briefly, the frequency uncertainty is defined by the full width at half maximum (FWHM) of the periodogram peak, which is equal to the inverse of the temporal baseline. However, such an assignment of uncertainty is fraught for sparse sampling over a large temporal baseline, which amounts to six years in our case. This is because the presence of numerous adjacent peaks in the window function makes it impossible to identify the true peak. Therefore, a more pragmatic choice for the uncertainties of f₁ and f₂ in this crossepoch analysis would be the FWHM of the envelope of the first peak of the window function, which is approximately 0.2 h⁻¹ (as shown in Fig. 2). Using this definition of FWHM and Eq. (52) by VanderPlas (2018), the uncertainties in f₁ and f₂ were computed: f₁ = 0.325 ± 0.006h⁻¹ and f₂ = 1.271 ± 0.006h⁻¹.

Fig. 3 Histogram of peak LS power calculated from light curves at random sky locations in the vicinity of J1019+65.

2.2.2 Fidelity of f₁ and f₂

Despite a sufficiently small value of FAP for the significant frequencies, a periodic signature in an LS periodogram could still be spurious if they are caused by, for example, image artefacts such as sidelobes from nearby bright sources. To test the fidelity of the two detected periodic signatures of the radio light curve of J1019+65 at P₁ ≈ 3.077 ± 0.057 h and P₂ ≈ 0.787 ± 0.006 h, we used the following method. We first picked a random⁴ sky location in the field of J1019+65 and sample the corresponding pixel in the LOFAR Stokes V image, which involves computing the cross-epoch light curve and the LS periodogram of this random pixel. The LS power (i.e. height) of the most significant (i.e. highest) peak in this periodogram is then recorded, before this trial is repeated many more times in order to get an ensemble of the maximum peak height of each random-pixel power spectrum. The idea behind this test is that if the periodicities detected in J1019+65’s light curve are, in fact, caused by sidelobes of nearby bright sources, other sky locations in the field should also show spurious periodic signatures of similar significance in their periodograms. Thus, we performed 2396 trials (each on a different pixel) and determined that < 0.3% of all random-pixel power spectra has a peak that exceeds the height of the f₂ peak, and ≪ 0.3% of the f₁ peak, shown in Fig. 3. We can thus empirically eliminate the possibility that the f₁ and f₂ peaks were caused by image artefacts.

Another test of the fidelity of f₁ and f₂ is to randomly shuffle the data points of the original light curve and compute the LS periodogram of the shuffled light curve. The maximum peak height of the power spectrum was then recorded each time we shuffled the data points. We performed 1000 trials and determined that <0.3% of all flux-density-shuffled power spectra has a peak that exceeds the height of the f₂ peak and ≪ 0.3% of the f₁ peak (as shown in Fig. 4).

As shown in Fig. A.4, we can see that unlike the f₁ signature, which is only present in the 120-136 MHz, the f₂ signature shows up still at higher frequencies of 136-168 MHz. This implies a difference between the spectral shape of the radio emission that is responsible for the f₁ and f₂ signatures.

Fig. 4 Histogram of peak LS power calculated by randomly shuffling the measured light curve values, while keeping the same time sampling window.

Fig. 5 Stokes V cross-epoch (2017-2023) light curve of J1019+65 phase-folded at a period of P1 = 3.0742 h, corresponding to the f1 peak seen in Fig. 2. The light curve is binned at intervals of 0.04 in phase and one full period corresponds to a phase interval of unity. The choice of starting phase is arbitrary. The green line represents the best sinusoidal fit using the computed LS model. The J1019+65 radio burst shown in Fig. A.3 is also plotted as a red line (with its width representing the duration of the burst) to indicate at which phase relative to the phase-folded light curves does it correspond to. Note: the burst shown in red only occurs in one epoch out of the whole dataset. The dashed black line at 0 mJy is drawn for clarity.

2.2.3 Phase-folded light curve and stacked image

Figures 5 and 6 show the cross-epoch light curve when phase-folded with each of two significant periods: P₁ = 3.07 h and P₂ = 0.79 h, respectively, derived from the f₁ and f2 peaks discussed in Sect. 2.2.1. From the figures, we can see that both phase-folded light curves are mostly below the zero line, implying that J1019+65 was producing a significant amount of phase-averaged (negative) Stokes V flux density in sub-mJy scales. Since a single normal LOFAR observation can only reach a 1σ sensitivity of ≈0.2 mJy in an 8h exposure (Callingham et al. 2023), a single epoch of observation would not be able to reveal the periodic emission. However, the intrinsic variability of the periodic emission means that it can occasionally be detected in single-epoch observations. Thus, the detection of this persistent emission (derived from LOFAR observations that comprise 52 h worth of data) is entirely consistent with the non-detections in the 4h or 8h observation of individual LOFAR epoch, but with single-epoch detections in the 2017-03-08 and 2021-02-03 epochs.

More importantly, it appears that this emission was detected mainly within a certain phase φ range. For example, as shown in Fig. 5, J1019+65 has significantly more negative Stokes V emission both on average (with a weighted mean S̄_V = -0.136 ± 0.029 mJy) and when φ₁ takes a value between 0.3 and 0.8 approximately, reaching its peak at φ₁ ≈ 0.55. We see the same phenomena in Fig. 6. However, we note that the quasi-sinusoidal profile of the data in Figs. 5 and 6 should be interpreted with caution, as the LS algorithm specifically finds periods that yield a sinusoidal phase-folded light curve. In other words, it may be possible to find a period very close to the LS peak that yields a different (non-sinusoidal) period signature.

As confirmation that the emission detected in our LS analysis originates from J1019+65, Figure 7 shows the inverse variation-weighted average image of J1019+65 made by stacking Stokes V images of individual LOFAR epochs from 2017 to 2023. Here, we detected a source with flux density of -0.161 ± 0.021 mJy near the expected location of J1019+65 in the stacked image. The proper motion-corrected location of J1019+65 derived from IR observations (cross marker in lime colour; Kirkpatrick et al. 2019) is approximately 0.5" away from the radio centroid, which has roughly the same value as the Gaussian-equivalent standard deviation (σ ≈ 0.5") in the absolute radio astrometry of LOFAR (Shimwell et al. 2022). As such, we are confident of the association between J1019+65 and the radio source in the stacked image. In comparing the stacked image in Stokes V (Fig. 7) with that in Stokes I (Fig. 8), we can see that the source is ≈100% circularly polarised. This result is consistent with ECM emission.

Fig. 6 Same as Fig. 5, but instead with a shorter period of P2 = 0.7866 h, corresponding to the f2 peak seen in Fig. 2.

Fig. 7 Stacked Stokes V image of J1019+65 with inverse-variance weighting from all LOFAR data (2017-2023), centred at the radio centroid (marked in magenta). The flux density of the radio centroid (peak pixel) is around −0.16 mJy. The 1σ noise of the image is around −0.02 mJy The proper-motion-corrected location of J1019+65 derived from IR observation is marked by a lime cross, approximately 0.5" away from the radio centroid.

Fig. 8 Same as Fig. 7, but in Stokes I.

3 Results and discussion

3.1 Handedness and emission geometry

As seen in Figures 5 and 6, the emission of J1019+65 occurs in specific phases and has a preferred Stokes V handedness. We note that we define the sign of circularly polarised emission as left-handed circularly polarised emission minus right-handed circularly polarised emission (Callingham et al. 2023). The bright burst arrives at a phase corresponding to the peak in the rotationally modulated emission. A preference in both the emission polarisation and the rotational phase of the radio source was also observed in the case of Jupiter (Marques et al. 2017). In Jupiter, right-hand circular polarised emission (stemming from the northern hemisphere) is more frequent than its left-hand southern counterpart due to the stronger magnetic field amplitude in the northern hemisphere (Zarka et al. 1996; Bonfond et al. 2013). The preferential detection of Jovian radio emission at certain rotational phase is due to ECM beaming, as discussed in Sect. 2.1. Therefore, the phenomenology of rotational modulated emission in J1019+65 shares several similarities with the Jovian case, suggesting a Jupiter-like radio engine operating in J1019+65.

The origin of the persistent emission of J1019+65 can therefore be interpreted as follows: the magnetic field lines of J1019+65 host ECM emission in a large region covering all azimuthal sectors. As the magnetosphere of J1019+65 rotates, the ECM emission beams from different azimuthal sectors come into and out of our line of sight. The rotationally modulated emission of J1019+65 detected within a particular phase is simply the collective ECM emission from all the visible beams during said time frame and its maximum magnitude occurs when the maximum range of sectors is visible to the observer. Typically, magnetic obliquity restricts the visible period to a small fraction of the rotation period and could make the emission from both magnetic poles visible at some fraction of the rotation period (Williams et al. 2017). In contrast, J1019+65 shows a persistent emission of the same handedness. This could indicate that the magnetic obliquity is smaller than the angular thickness of the emission cone, thereby keeping the source from a single magnetic pole in view for most of its rotation.

The detected radio burst of J1019+65 occurred near the peak of the phase-folded light curves, as seen in Figs. 5 and 6. The location of the burst is also consistent with a Jupiter-like engine because we have the highest chance of seeing anomalously bright ECM-induced radio bursts when the persistent emission is detected at the highest significance, which should correspond to the rotational phase at which the maximum amount of ECM conical beams are towards our line of sight. We also note that the flux density of the burst of J1019+65 is around several mJy, which is ~ 1 order of magnitude brighter than its persistent emission. This is consistent with what was seen in the case of Jupiter as well (Zarka 2007).

3.2 The true nature of f₂

The presence of two distinct frequencies in the periodogram is unusual and requires an explanation. There are four possibilities for the secondary peak at f₂ : (i) a spurious frequency caused by aliasing; (ii) some harmonic of f₁ ; (iii) a true frequency corresponding to the rotation period of the other brown dwarf; or (iv) a true frequency due to satellite interaction. We tackled each possible explanation individually in this order.

Aliasing. As with any other method based on Fourier transform, the LS periodograms can have aliased power. In general, for an observation whose window function has a strong power at δf, we can expect an observed peak, f_peak, to be an alias if there are peaks at |f_peak ± nδf|, where n ∈ ℕ⁺ is small (n ~ 1). With that in mind, we inspected the window function for our cross-epoch analysis (plotted in A.5) and identified the different prominent peaks. For example, our ‘sampling’ rate of 4 min corresponds to δf = 15 hour⁻¹, so clearly this cannot be responsible for creating aliases of f₁ at f₂. We were unable to identify a peak that could explain the value of f₂ and f₁ based on aliasing.

Harmonic. If the signal is periodic but not sinusoidal, its power in the periodogram will be present at the fundamental frequency f also at higher harmonics nf, where n ∈ {2,3,4,...}. This is especially relevant for Jupiter-like radio-loud objects such as J1019+65, since (i) we anticipate a light curve that resembles an on-off function and (ii) there are certain magnetic configurations that allow a brown dwarf to pulsate more than once per rotation; Jupiter typically emits two pulses for every rotation. As f₂/f₁ = 3.91, the closest harmonic number is n = 4. This results in a residual of 4f₁ - f₂ = 0.029 ± 0.025 hour⁻¹. In comparison, the frequency uncertainty of f₁ in the original LoTSS DR2 data shown in Fig. A.3 has a value that is about three times smaller than this residual. Moreover, there are no signs of a significant harmonic peak at 2f₁ or 3 f₁ in the LS periodogram (see Fig. 2), which is contrary to what was expected if f₂ is indeed the fourth harmonic. Furthermore, the LS power at f₂ can exceed that at f₁ in some spectral bands as seen in Fig. A.4. If f₂ were to be a harmonic of f₁, f₂ should not have any LS power in the 136 -168 MHz band - and yet it does. For the above reasons, f₂ is unlikely to be a higher harmonic of f₁.

Rotation of the second brown dwarf. Since J1019+65 is a brown dwarf binary, f₂ may correspond to the rotation rate of the ‘second’ brown dwarf, given that f₁ corresponds to the ‘first’ brown dwarf. In this scenario, both brown dwarfs in J1019+65 emit persistent ECM emissions in the LOFAR band; since LOFAR does not have sufficient angular resolution, the emissions of both brown dwarfs were superposed, thus creating a light curve with two periodic signatures of P₁ = 3.07 h and P₂ = 0.79 h. Using the estimation of J1019+65's physical parameters by (Vedantham et al. 2023), we calculated the rotation break-up period to be 0.45 h for the more massive T5.5 brown dwarf and 0.51 h for the lighter T7.0 brown dwarf in J1019+65, so a 0.79 h period is physically possible. If true, the second brown dwarf would be the fastest rotating brown dwarf ever detected (cf. Tannock et al. 2021). However, interpreting P₂ as the rotation period can be fraught (see e.g. Williams et al. 2017) without any independent confirmation from IR photometry. For instance, it is also plausible that the radio period is a harmonic of the true rotation period. Regardless, for this scenario to work, both radio emitting objects must have comparable radio luminosities, which is somewhat contrived. Follow-up IR observations on J1019+65 may be the key to determining the nature of f₂ ; a detection of these (P₁ and P₂) periodicities in the IR light curve of J1019+65 due to the anticipated rotational modulation would definitively confirm the authenticity of f₂ . Moreover, using adaptive optics, the two components of J1019+65 were easily resolved (Vedantham et al. 2023). Therefore, further IR observations could also help determine which brown dwarf has the longer or shorter period.

Satellite interaction. Based on the Jovian analogy, an intriguing possibility is that the two periodicities come from two emission components: one modulated by the rotation of the emitter (0.79 h) and another modulated by the orbital period of a companion satellite (3.07h). At this juncture, we must consider whether such short rotation periods are even feasible. To answer this, we first checked the Roche limit by computing the net force on a test particle at the surface of the satellite. To avoid disintegration, the centrifugal force plus the gravitational force due to the satellite’s mass must exceed the gravitational force of the brown dwarf at the surface of the satellite. The calculations show that if the brown dwarf mass is 30M_J, then disintegration of the satellite can only be prevented if the mean density of the satellite exceeds ≈4 g cm⁻³. Such a high density is unlikely for a rocky planet, but possible for a massive gas giant planet or another brown dwarf. Next, we considered the tidal evolution of the orbit. The putative satellite is beyond the co-rotation radius, which means its orbital distance will increase due to tidal forces. Equations for the approximate timescale for such an evolution have been derived by Bolmont et al. (2011). We used their Eq. (2) with an eccentricity of 0, primary and secondary masses of 30 M_J and 5 M_J respectively, primary and secondary radii equal to the Jovian radius and an orbital period of 3 h as characteristic values. We find that the orbital evolution timescale is under 10⁶ yr, even if we assume that the putative planet is perfectly tidally locked to the brown dwarf. Because we have no indication that J1019+65 would be classed as a very young system, the short tidal evolution timescale is a serious problem for the existence of a planet in such a short orbit.

3.3 Radio-loud fraction and duty ratio of brown dwarfs

J1019+65 was the only brown dwarf system within 25 pc detected in LoTSS DR2 among hundreds of catalogued brown dwarfs and within the survey footprint. Therefore, prior to our follow-up radio observations, it was reasonable to assume that J1019+65 has some atypical property (e.g. binarity, spectral type, and rotation rate) that led to a detectable level of radio emission. Our follow-up exposures combined with the 8h long LoTSS exposure amount to 52 hours over two years, where we only detected one radio burst lasting for several minutes. Because J1019+65 was identified due to its bright 4 min burst, we are led to an interesting question regarding whether J1019+65 had been discovered because it is atypical. Alternatively, J1019+65 could, in fact, be typical, but also happened to be in a high-flux-density state during the LoTSS DR2 observations. In other words, here we consider whether radio-emitting brown dwarfs are rare or whether radio emission is actually common among brown dwarfs, but flux-limited surveys end up picking up only the small numbers of brown dwarfs that happen to be in a high-activity state.

3.3.1 Rarity versus activity

We start this subsection with an explanation of the ‘rarity’ and ‘activity’ of J1019+65 within a mathematical context: we let F_radio = N_radio/N_total and 〈D〉_BD be the fraction of brown dwarfs that are ‘radio-loud’ and the expected (average) duty ratio of a ‘radio-loud’ brown dwarf, respectively. In this context, we define a ‘radio-loud’ brown dwarf as one that is capable of emitting ECM emission detectable in a typical 8h LoTSS DR2 exposure, which has a median continuum noise of 95 μJy in the 120168 MHz band. For our 4h observations, this corresponds to a median continuum noise of around 134 μJy. Therefore, the radio detectability of a brown dwarf mainly hinges on whether said brown dwarf is capable of generating such ECM emissions in the first place and, if so, how often we expect the brown dwarf to generate them. F_radio encapsulates the former, while 〈D〉_BD the latter.

Intuitively, we can see that the two quantities of F_radio and 〈D〉_BD correspond to the rarity and activity of a radio-loud brown dwarf, respectively: a small value of F_radio implies that radio-loud brown dwarfs are rare, while a small value of 〈D〉_BD implies that radio-loud brown dwarfs are mostly inactive, thus leading to a small number of detections in a survey.

Mathematically, the duty ratio of a radio-loud brown dwarf is defined as D = P_on/P_BD, where P_on is the length of time for which the brown dwarf’s radio signal is ‘on’ in a single rotation period P_BD ≡ P_on + P_off of said brown dwarf. In the context of this analysis, the radio signal of J1019+65 at time t can be considered as ‘on’ when its ECM emission exceeds a certain detection threshold during said time, t. This definition is equivalent to the following statement: D_J1019 represents the probability of detecting significant (i.e. above the detection threshold) radio emission from J1019+65 at any given time (ignoring beaming for now). We shall revisit this statement in a follow-up analysis.

3.3.2 Constraints on 〈D〉_BD and F_radio

With the concept of rarity and activity mathematically defined, we can proceed to explaining how constraining the duty ratio of J1019+65 could ultimately lead to constraints on F_radio and 〈D〉_BD. We recall from Sect. 1 that by the virtue of an untargeted radio sky survey, the radio detection of J1019+65 in LoTSS does not suffer from selection bias that is inherent to targeted observations (where most other radio-loud methane dwarfs were found). Moreover, as highlighted in Sect. 2.2, the radio emission from J1019+65 matches the characteristics of the canonical Jovian periodic pulsed emission profile previously seen from brown dwarfs and ultracool dwarfs (e.g. Hallinan et al. 2007, 2008, 2015), meaning J1019+65 should be a typical radio-loud brown dwarf in this regard. Both reasons are crucial if we were to justify any estimation of F_radio and 〈D〉_BD using the statistics of J1019+65. Therefore, we assumed that the magnetic field strength of J1019+65 is not atypical for the underlying brown dwarf population.

Now, because J1019+65 is at d = 23.3 pc, we took all known brown dwarfs within 25 pc as our parent sample. This gave us a total of 106 brown dwarfs (spectral types L0-Y2) within the LoTSS DR2 footprint based on the volume-complete sample by Best et al. (2024). Therefore, this naturally leads to the following question: given that the total number of brown dwarfs in LoTSS DR2 is N_LoTSS = 106, we need to determine the probability of detecting only one radio-loud brown dwarf (i.e. J1019+65). In such a case, even if all brown dwarfs were to be radio-loud, only a fraction would be visible from the Earth due to beaming. And since we cannot differentiate between a brown dwarf that is not radio-loud and a radio-loud brown dwarf that is not beaming towards us, the observed fraction of radio-loud brown dwarfs is therefore F^′_radio ≡ F_radio × F_beam, where F_beam is the fraction of (radio-loud) brown dwarfs that beam towards us (and thus visible to us).

Finally, under these assumptions, we can proceed with determining the probability of detecting only k = 1 brown dwarf in LoTSS DR2. This should follow a Poisson process with the rate parameter, λ, expressed as $$P(\lambda ,k=1)=\lambda e^{-\lambda },\lambda =N_{\mathrm{L}\mathrm{o}\mathrm{T}\mathrm{S}\mathrm{S}}\times {⟨D⟩}_{\mathrm{B}\mathrm{D}}\times F_{\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{i}\mathrm{o}}^{{}^{\prime }}.$$(1)

Already, computing P(λ) gives us the joint probability distribution function (PDF) of (D)_BD and F^′_radio. The singular LoTSS DR2 detection can only constrain the product of (D)_BD and F′_radio, meaning J1019+65 could have been the singular radio-loud object, or it could be just a typical brown dwarf that happened to be in a bright state during the LoTSS DR2 exposure. This is clearly illustrated in red in Fig. 9, where we compute the covariance of the posterior distribution of 〈D〉_BD and F'_radio based on Eq. (1).

To break this degeneracy between 〈D〉_BD and F′_radio, we now consider the non-detections of radio pulsation in our follow-up observational campaign of J1019+65. As explained in Sect. 3.3, we observed J1019+65 with LOFAR for a total of 52 hours, which corresponds to n ≈ 17 rotations of the radio emitter in J1019+65 based on the radio-derived period P₁ = 3.0742 hour. We note that we ignored P₂ in this analysis⁵. We recall that the phenomenon of radio pulsations is the consequence of rotational modulation due to beaming; intuitively, we expect to detect 1 pulse per rotation⁶ for a brown dwarf with D = 1. Therefore, the probability of detecting only x = 1 bursts given that the LOFAR observations covered n rotations of J1019+65 is $$P(p,n=17,x=1)=(\frac{n}{x})p^x(1-p{)}^{n-x}=17p(1-p{)}^{16},p=D_{\mathrm{J}1019},$$(2)

where $(\frac{n}{x})\equiv \frac{n!}{x!(n-x)!}$, is the binomial coefficient. As explained earlier in this section, we can interpret J1019+65 as an accurate representation of a radio-loud brown dwarf in this sample owing to the virtue of an untargeted survey. This is a powerful statement as it implies $D_{\mathrm{J}1019}\simeq {⟨D⟩}_{\mathrm{B}\mathrm{D}}$; hence, we can constrain 〈D〉_BD simply by estimating the duty ratio of J1019+65. As we can see from Eq. (2), the binomial distribution depends solely on the duty ratio (in contrast to Eq. (1)).

The likelihood of our LoTSS DR2 and follow-up observations is thus given by the product of Eqs. (1) and (2). The result is shown in blue in Fig. 9; compared to the data in red that only contains the LoTSS data, the extra restriction provided by Eq. (2) greatly improve the constraints on 〈D〉_BD and consequently F′_radio. Figure 10 shows a detailed version of the marginal PDFs presented in Fig. 9. Both curves in Fig. 10 are unimodal. Therefore, using the method of maximum likelihood estimation (MLE), we obtained ${⟨D⟩}_{\mathrm{B}\mathrm{D}}={0.030}_{-0.030}^{+0.034}$ and $F_{\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{i}\mathrm{o}}^{{}^{\prime }}={0.088}_{-0.088}^{+0.168}$. Moreover, we can reject any value of 〈D〉_BD > 0.177 and F′_radio > 0.959 with ≥ 95% confidence. The rejection of a large 〈D〉_BD implies that in order to capture the moment of a radio burst from a brown dwarf, we must either (i) observe a radio-loud brown dwarf for many rotations or (ii) observe many brown dwarfs at the same time (e.g. a sky survey). Therefore, if LoTSS were to be repeated, we would expect to find new radio-loud brown dwarfs with detected burst. Furthermore, we can estimate F_beam = Ω/4π ≈ 12.7% based on the instantaneous solid angle of the measured Jovian emission beam Ω = 1.6 sr (Zarka et al. 2004). We note that this value serves as a lower limit of F_beam; for brown dwarfs with magnetic obliquity, their beaming pattern would thus sweep across the sky, as the magnetic pole wobbles around when observing over a whole rotation period. This estimate gives us F_radio = F′_radio/F_beam ≲ 0.691, meaning that brown dwarfs such as J1019+65 cannot be particularly rare. In comparison, our estimation of F_radio is consistent with the work by Kao & Pineda (2024), which showed that the fraction of ultracool dwarf binary systems that emit persistent radio emission is ${0.56}_{-0.11}^{+0.11}$. We note that the focus of their work is on the incoherent (gyro)synchrotron emission stemming from radiation belts of brown dwarfs, rather than coherent ECM emission.

Fig. 9 Corner plot showing the covariance of the joint probability distribution function (PDF) and marginal PDFs of the two parameters: the observed fraction of radio-loud brown dwarfs, F′radio, in our sample and its expected duty ratio, 〈D〉BD. Information plotted in red only considers the Poisson probability (Eq. (1)) derived from LoTSS DR2 observation, whereas information in blue also considers the binomial distribution (Eq. (2) multiplied by Eq. (1)) derived from the follow-up LOFAR data. The 2D distribution have contours indicating the regions containing 68% and 95%.

Fig. 10 Marginal probability distribution functions of 〈D〉BD and F′radio. In each plot, the relative likelihood is defined such that it is unity at the dashed black line, corresponding to the most probable value of each parameter, i.e. ${⟨D⟩}_{\mathrm{B}\mathrm{D}}={0.030}_{-0.030}^{+0.034}$ and $F_{\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{i}\mathrm{o}}^{{}^{\prime }}={0.088}_{-0.088}^{+0.168}$. The error bars of the parameters correspond to 68% confidence interval (CI), as represented by the blue-shaded region. The parameter value within orange shaded region can be rejected with ≥95%.

4 Conclusions and outlook

Our follow-up observational campaign of J1019+65 with LOFAR led to 52 hours of observations spanning close to six years. We only detected a periodic signal at 120-136 MHz in one out of ten individual LOFAR epochs using LS periodogram, while the rest of the epochs were all non-detections. Moreover, we did not detect any radio burst from J1019+65 similar to the original pulse by Vedantham et al. (2023) in any of our follow-up observations.

Despite the lack of radio bursts, we found that J1019+65 shows persistent circularly polarised emission that is periodically modulated. We find two distinct periodicities: a new periodic signature of P₂ = 0.79 h and the 3h period initially identified by Vedantham et al. (2023). We find that the periodogram peak corresponding to the new period in the LS periodogram is not a spurious signal due to aliasing, harmonic structure, or the window function. If f₂ is indeed the rotation rate of the second brown dwarf in J1019+65, it would have the shortest rotation period of any brown dwarf to date; however, confirming such an association between the two peaks and the two objects of the binary would require an independent conformation, possibly from IR band photometry, in the future. Moreover, the ratio between the two inferred periods is about 4; we do not have systematic samples of T dwarf binary periods for comparison, however, in low-mass stellar binaries, the two components tend to have similar rotation periods (see e.g. Harding et al. 2013; Stauffer et al. 2016). If the same trend continues for cold brown dwarfs, this period ratio would be unusual.

Based on the non-detection of bursts in most of our LOFAR observations, we computed the posterior distribution for the fraction of radio-loud brown dwarfs and the radio duty cycle of a typical brown dwarf, allowing us to constrain the rarity of a radio-loud brown dwarf and the average activity of a radio-loud brown dwarf. As such, we were able to statistically determine that radio-loud brown dwarfs in the LOFAR band would not be expected to be particularly rare (F_radio ≲ 0.691); we simply need to stare at brown dwarfs for long enough to counteract to their low average duty ratio ($⟨D⟩={0.030}_{-0.030}^{+0.034}$). In comparison, Kao & Pineda (2024) showed that the fraction of ultracool dwarf binary systems that emit persistent radio emission is ${0.56}_{-0.11}^{+0.11}$ which is consistent with our estimation, although it should be noted that the persistent emission in their work stems from synchrotron radiation - and not ECM emission.

Acknowledgements

T.W.H.Y. and H.K.V. acknowledge funding from EOSC Future (Grant Agreement no. 101017536) projects funded by the European Union’s Horizon 2020 research and innovation programme. HKV acknowledges funding from the European Research Council via the starting grant ‘STORMCHASER’ (grant number 101042416) and from the Dutch research council under the talent programme (Vidi grant VI.Vidi.203.093). JRC acknowledges funding from the European Union via the European Research Council (ERC) grant Epa-phus (project number 101166008). LOFAR is the Low Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage facilities in several countries, which are owned by various parties (each with their own funding sources), and which are collectively operated by the ILT foundation under a joint scientific policy. The ILT resources have benefited from the following recent major funding sources: CNRS-INSU, Observatoire de Paris and Université d’Orléans, France; BMBF, MIWF-NRW, MPG, Germany; Science Foundation Ireland (SFI), Department of Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The Science and Technology Facilities Council, UK; Ministry of Science and Higher Education, Poland; The Istituto Nazionale di Astrofisica (INAF), Italy.

References

Appendix A Supplementary figures

Fig. A.1 Follow-up LOFAR observations of J1019+65 carried out on the first three dates (excluding 2021-02-23).

Fig. A.2 Follow-up LOFAR observations of J1019+65 carried out on the remaining six dates.

Fig. A.3 LOFAR observation of J1019+65 carried out on 2017-03-08 (i.e. LoTSS DR2 epoch). Top panel: Stokes V light curves with a cadence of 4 min at different frequencies. The shaded region represents ±1σ uncertainty. The negative ≈ 10mJy burst at 120-136 MHz (enclosed in red) has a signal-to-noise ratio of > 7. The light curve at 152-168 MHz is omitted for the sake of clarity. The top axis represents the time of observation in Modified Julian Dates (MJD). The dashed black line at 0mJy is drawn for clarity. Bottom panel: LS periodogram of the radio light curves. The three dashed red lines represent the necessary LS power (i.e. peak height) to achieve a false-alarm probability (FAP) of 5%, 1%, and 0.3%. The grey curve represents the (negative) LS periodogram of the window function, which is a light curve with the same timestamps as the original curve, but the flux density of each data point is replaced with unity (i.e. a flat light curve).

Fig. A.4 LS periodograms of the LOFAR cross-epoch radio light curve, which includes data from both the original LoTSS DR2 epoch and our ten LOFAR follow-up epochs.

Fig. A.5 Power spectrum of the sampling window function for the full LOFAR dataset (cross-epoch analysis of Sect. 2.2). The light curve sampling cadence of 4 min has an identifiable peak.

All Tables

All Figures

Fig. 1 LOFAR observation of J1019+65 carried out on 2021-02-33. Top panel: Stokes V light curves with a cadence of 4 min at different frequencies. The shaded region represents ±1σ uncertainty. The light curve at 152-168 MHz is omitted for the sake of clarity. The top axis represents the time of observation in Modified Julian Dates (MJD). The dashed black line at 0 mJy is drawn for clarity. Bottom panel: LS periodogram of the radio light curves. The three dashed red lines represent the necessary LS power (i.e. peak height) to achieve a false-alarm probability (FAP) of 5%, 1%, and 0.3%. The grey curve represents the (negative) LS periodogram of the window function, which is a light curve with the same timestamps as the original curve, but the flux density of each data point is replaced with unity (i.e. a flat light curve). The two relevant frequencies, f1 = 0.315 ± 0.013 h−1 and f2 = 1.247 ± 0.013 h−1, are indicated by the black arrows in the LS periodogram. The 120-136 MHz (orange) peak at f1 is consistent with the original frequency discovered by Vedantham et al. (2023), while the relevancy of the f2 is discussed in Sect. 2.2.2. The frequency uncertainties were computed using Eq. (52) by VanderPlas (2018).

In the text

Fig. 2 LS periodograms of the LOFAR cross-epoch radio light curve, which includes data from both the original LoTSS DR2 epoch and our ten LOFAR follow-up epochs. The spectrum reveals the peak at f1 ≈ 0.325 h−1 originally discovered by Vedantham et al. (2023). Another significant peak not previously discovered is detected at f2 ≈ 1.271 h−1, with a false alarm probability (FAP) of less than 0.3%. The locations of these two peaks are indicated by the black arrows. The three dashed red lines represent the FAP of different values. The grey curve represents the LS periodogram of the window function, which is a light curve with the same timestamps as the original curve; however, the flux density of each data point is replaced by unity (i.e. a flat light curve). The inset plot in each panel shows a zoomed-in version of the observed J1019+65 power spectrum in the vicinity of the tallest LS peak, with the same three FAP dashed red lines.

In the text

	Fig. 3 Histogram of peak LS power calculated from light curves at random sky locations in the vicinity of J1019+65.
In the text

	Fig. 4 Histogram of peak LS power calculated by randomly shuffling the measured light curve values, while keeping the same time sampling window.
In the text

Fig. 5 Stokes V cross-epoch (2017-2023) light curve of J1019+65 phase-folded at a period of P1 = 3.0742 h, corresponding to the f1 peak seen in Fig. 2. The light curve is binned at intervals of 0.04 in phase and one full period corresponds to a phase interval of unity. The choice of starting phase is arbitrary. The green line represents the best sinusoidal fit using the computed LS model. The J1019+65 radio burst shown in Fig. A.3 is also plotted as a red line (with its width representing the duration of the burst) to indicate at which phase relative to the phase-folded light curves does it correspond to. Note: the burst shown in red only occurs in one epoch out of the whole dataset. The dashed black line at 0 mJy is drawn for clarity.

In the text

	Fig. 6 Same as Fig. 5, but instead with a shorter period of P2 = 0.7866 h, corresponding to the f2 peak seen in Fig. 2.
In the text

Fig. 7 Stacked Stokes V image of J1019+65 with inverse-variance weighting from all LOFAR data (2017-2023), centred at the radio centroid (marked in magenta). The flux density of the radio centroid (peak pixel) is around −0.16 mJy. The 1σ noise of the image is around −0.02 mJy The proper-motion-corrected location of J1019+65 derived from IR observation is marked by a lime cross, approximately 0.5" away from the radio centroid.

In the text

	Fig. 8 Same as Fig. 7, but in Stokes I.
In the text

Fig. 9 Corner plot showing the covariance of the joint probability distribution function (PDF) and marginal PDFs of the two parameters: the observed fraction of radio-loud brown dwarfs, F′radio, in our sample and its expected duty ratio, 〈D〉BD. Information plotted in red only considers the Poisson probability (Eq. (1)) derived from LoTSS DR2 observation, whereas information in blue also considers the binomial distribution (Eq. (2) multiplied by Eq. (1)) derived from the follow-up LOFAR data. The 2D distribution have contours indicating the regions containing 68% and 95%.

In the text

Fig. 10 Marginal probability distribution functions of 〈D〉BD and F′radio. In each plot, the relative likelihood is defined such that it is unity at the dashed black line, corresponding to the most probable value of each parameter, i.e. ${⟨D⟩}_{\mathrm{B}\mathrm{D}}={0.030}_{-0.030}^{+0.034}$ and $F_{\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{i}\mathrm{o}}^{{}^{\prime }}={0.088}_{-0.088}^{+0.168}$. The error bars of the parameters correspond to 68% confidence interval (CI), as represented by the blue-shaded region. The parameter value within orange shaded region can be rejected with ≥95%.

In the text

	Fig. A.1 Follow-up LOFAR observations of J1019+65 carried out on the first three dates (excluding 2021-02-23).
In the text

	Fig. A.2 Follow-up LOFAR observations of J1019+65 carried out on the remaining six dates.
In the text

Fig. A.3 LOFAR observation of J1019+65 carried out on 2017-03-08 (i.e. LoTSS DR2 epoch). Top panel: Stokes V light curves with a cadence of 4 min at different frequencies. The shaded region represents ±1σ uncertainty. The negative ≈ 10mJy burst at 120-136 MHz (enclosed in red) has a signal-to-noise ratio of > 7. The light curve at 152-168 MHz is omitted for the sake of clarity. The top axis represents the time of observation in Modified Julian Dates (MJD). The dashed black line at 0mJy is drawn for clarity. Bottom panel: LS periodogram of the radio light curves. The three dashed red lines represent the necessary LS power (i.e. peak height) to achieve a false-alarm probability (FAP) of 5%, 1%, and 0.3%. The grey curve represents the (negative) LS periodogram of the window function, which is a light curve with the same timestamps as the original curve, but the flux density of each data point is replaced with unity (i.e. a flat light curve).

In the text

	Fig. A.4 LS periodograms of the LOFAR cross-epoch radio light curve, which includes data from both the original LoTSS DR2 epoch and our ten LOFAR follow-up epochs.
In the text

	Fig. A.5 Power spectrum of the sampling window function for the full LOFAR dataset (cross-epoch analysis of Sect. 2.2). The light curve sampling cadence of 4 min has an identifiable peak.
In the text