© The Authors 2025

1. Introduction

Fast radio bursts (FRBs) are microsecond- to millisecond-long bursts of radio emission of extragalactic origin (see Petroff et al. 2019, 2022, for review). Over the almost two decades since the discovery of the first FRB (Lorimer et al. 2007), multiple FRB emission theories have been proposed (for the live catalog see Platts et al. 2019). To date no consensus about their exact origin has emerged, although neutron star progenitors currently appear to be favored.

It is even possible that several classes of progenitors emit FRBs. At the moment, the most obvious empirical demarcation between these types is the dichotomy between one-off sources and FRB repeaters. Only the former are potentially cataclysmic. These two classes show statistical differences in the spectro-temporal properties of the bursts (Pleunis et al. 2021a). In comparison to one-off sources, repeaters offer much more information about their circum-burst environment and the burst emission mechanism: they can be localized with much better precision and the pulses provide dynamical estimates of the electron content and magnetic field in the vicinity of the emitting plasma. Any emission theory must explain the entire distribution of burst fluences, and the various spectro-morphological properties of pulses that originate in the same local environment.

While some FRBs are only seen to repeat a handful of times, FRB 20201124A is a veritable FRB factory: it is capable of emitting prolifically, and datasets covering it contain hundreds of pulses. It is thought to be located in a dynamically evolving magnetized environment, as suggested by irregular rotation measure (RM) variations on short timescales and the presence of Faraday conversion (Xu et al. 2022). FRB 20201124A is notorious for its high but exceedingly variable pulse emission rate. The source was first detected at the end of 2020 by the Canadian Hydrogen Intensity Mapping Experiment (CHIME), but only after about 40 hours of earlier observations of the same field contained no detections (Lanman et al. 2022). By 2021 March-May FRB 20201124A had entered a high-activity phase, hereafter called the Spring 2021 (S21) epoch¹, reaching rates almost 50 bursts per hour, as observed by the Five-hundred-meter Aperture Spherical radio Telescope (FAST) at 20 cm wavelength (Xu et al. 2022). The S21 FAST burst sample was complemented by observations performed with CHIME, the Effelsberg 100 m and the Parkes 64 m dishes, the upgraded Giant Metrewave Radio Telescope (uGRMT), and the Australian Square Kilometre Array Pathfinder (ASKAP) (Lanman et al. 2022; Hilmarsson et al. 2021; Marthi et al. 2022; Kumar et al. 2022).

Bursts from FRB 20201124A exhibit a high degree of circular and linear polarizations, with predominantly flat position angle (PA) curves. These polarization properties hint at a magnetospheric origin (Jiang et al. 2022). Such an origin agrees with a rotating neutron-star progenitor hypothesis. However, despite extensive searches, no periodicity has been found in the high number of bursts, over a broad range of trial periods spanning from milliseconds to days (Niu et al. 2022; Du et al. 2024).

The FAST observations indicate that the S21 activity epoch ended abruptly between 2021 May 26 and 29 (Xu et al. 2022, observations at 1250 MHz). On May 27 CHIME/FRB recorded one more burst from this FRB at the lower frequencies of 400–800 MHz (Lanman et al. 2022), followed by a bright burst at 1350 MHz detected by the Stockert telescope on May 28, during a three-day gap in FAST coverage (Kirsten et al. 2024). The gaps in the observing schedules of these three telescopes did not allow further constraints on any possible frequency-dependent boundaries of the activity phase, an effect seen in one other repeating FRB, FRB 20180916B (Pastor-Marazuela et al. 2021). The end of the activity window was, however, also covered by GMRT observations at 300–750 MHz, which are presented in this work. Strikingly, GMRT detected several low-frequency pulses on June 1, several hours after the FAST non-detections showed the higher-frequency bursts had already turned off.

Despite continued monitoring (Mao et al. 2022; Trudu et al. 2022), these June 1 pulses were the only signs of activity detected from FRB 20201124A until the pulses reappeared at some time before 2021 September 21 (Main et al. 2021, also CHIME²). An extensive FAST monitoring campaign then observed an exponentially increasing FRB rate, reaching activity levels an order of magnitude higher than those Xu et al. reported earlier. The activity in this Fall 2021 (F21) epoch abruptly quenched again between September 28 and 29 (Zhou et al. 2022). No detections were published until late 2022 January (the “Winter 2022”, W22 epoch; Ould-Boukattine et al. 2022).

The repeater field was regularly scheduled in the Apertif-LOFAR³ Exploration of the Radio Transient sky (ALERT) survey, that started in 2019 at the Westerbork Synthesis Radio Telescope (WSRT; Maan & van Leeuwen 2017; van Leeuwen et al. 2023; Pastor-Marazuela et al. 2025). During the last observing run of this survey, in the first week of 2022 February, the three observing sessions toward FRB 20201124A yielded ten FRB detections.

This work describes the properties of those bursts detected within the ALERT survey in W22, and of GMRT bursts from the S21 and W22 activity windows. We include the analysis of dispersion measures (DMs), fluence distributions, (quasi-)periodicities, and the scintillation of these bursts.

The evidence we find for a frequency-dependent activity window, together with the energetics and morphological properties of the recorded pulses – both in our sample and in the large FAST sample – provide a unique test of the hypothesis that FRBs originate from low-twist ultra-long-period (ULP) magnetars (Wadiasingh & Timokhin 2019; Wadiasingh et al. 2020; Beniamini et al. 2020; Caleb et al. 2022). That is interesting because a variety of scenarios have been proposed to explain how neutron stars might form FRBs, but testing and distinguishing these observationally is challenging.

In the ULP hypothesis, FRBs are generated via a pulsar-like emission mechanism in the magnetospheres of very slowly rotating magnetars. In these sources, the non-potential magnetosphere (i.e., one in which currents flow) is characterized by an unusually weak large-scale magnetic field twist, much weaker than that of commonly observed galactic magnetars. This results in a low plasma density on the closed magnetic field lines. In such charge-starved conditions, deformations in the crust that dislocate the footpoints of the magnetic field lines generate strong transient electric fields. This, in turn, leads to avalanche pair production and, ultimately, the emission of FRBs via a pulsar-like coherent mechanism.

The low-twist magnetar hypothesis of FRB generation makes specific predictions for the times of arrival (TOAs) of individual bursts, and for the quasi-periodicity of sub-bursts. In this work, we compare these predictions to existing observational evidence. Finally, we put constraints on the location and shape of active regions on the surface of the star by using two classical phenomenological models of radio pulsar emission: radius-to-frequency mapping and rotating vector models.

2. GMRT observations and analysis

2.1. Observations

We observed FRB 20201124A with the uGMRT (Gupta et al. 2017) on 2021 June 01, July 03 & 04, and 2022 February 05, 06 & 07, under Director’s Discretionary Time. Except for one session, FRB 20201124A was observed simultaneously in the band 3 (300–500 MHz) and band 4 (550–750 MHz) of GMRT, by combining the antennae in two different sub-arrays. These dual-band observations thus provided a frequency coverage of 300–750 MHz. The observing setup utilized the phased-array beam (PAB) sensitivity of typically 11 antennae in each of the two sub-arrays. At band 4, data were recorded at 4096 channels across 200 MHz bandwidth centered at 650 MHz, with a sampling time of 0.328 ms. At band 3, the data were coherently dedispersed in real-time at a DM of 410.9 pc cm⁻³, and then recorded to disk with 1024 sub-bands across 200 MHz bandwidth centered at 400 MHz, using a sampling time of 0.164 ms. The last session on 2022 February 07 was conducted only in band 4 utilizing a higher sensitivity obtained by combining 16 antennae in a single sub-array, and the PAB data were coherently dedispersed at a DM of 410.9 pc cm⁻³ in real-time and then recorded to the disk with 1024 sub-bands and a sampling time of 40.96 μs. During all the observations conducted in 2022, data from the incoherent-array beam (IAB) formed using the same antennae in the individual sub-arrays were also recorded simultaneously with the PAB. The availability of the IAB data facilitate a better excision of the radio frequency interference (RFI).

2.2. Preprocessing and search for radio bursts

The individual band data from all the sessions were processed through the following series of data reduction steps. For the band 4 observations, which did not employ the real-time coherent dedispersion, we used SIGPROC’s dedisperse to create 1024 sub-bands, by dedispersing sets of 4 channels to a single sub-band, using a DM of 410.9 pc cm⁻³. As mentioned above, all the band 3 observations were already recorded with 1024 coherently dedispersed sub-bands (the real-time coherent dedispersion removes only the intra-sub-band smearing). The individual band data from each of the sessions were then subjected to down-sampling from 16 bits to 8 bits per sample using digifil, and RFI excision using RFIClean⁴ (Maan et al. 2021) and rfifind from the pulsar search and analysis toolkit PRESTO (Ransom et al. 2002).

For the observations conducted in 2022, we formed the post-correlation beam (Roy et al. 2018) using the PABs and IABs (i.e., we subtracted the IAB from the PAB after appropriate scaling⁵). This post-correlation beam contains less red noise, mitigates some RFI and thus reduces false candidates while searching for bright pulses. This beam was processed through the same steps described above.

The sub-banded data were then incoherently dedispersed in steps of 0.05 and 0.2 pc cm⁻³, covering the DM range 400–425 pc cm⁻³ in 500 and 125 trial DMs, for band 3 and band 4 respectively, using prepdata from PRESTO. The above configurations limit the dispersive smearing to a maximum of 1 ms throughout the explored DM range for both bands. Each dedispersed time-series was searched for the presence of bright pulses above a signal-to-noise ratio (S/N) threshold of 8 and a maximum pulse-width of 50 ms using single_pulse_search.py. All the single pulse candidates were next grouped by looking for the brightest candidates within roughly 100 or 200 ms wide time windows and across all the trial DMs. Waterfall plots for all the grouped candidates were scrutinized by eye to identify the genuine bursts.

Each dedispersed time-series was also searched for periodic signals using accelsearch from PRESTO, with zmax set to 256. For each observation, all the periodic signal candidates were sifted for harmonics and duplicates at different DM trials, and the final candidates were folded using prepfold and the diagnostic plots were examined by human experts.

Overall, 16 bursts were detected, six during the 2021 June 01 session, one on 2022 February 6 and the rest on February 07. The bursts are labeled G01–G16 following their TOA order. Figure 1 shows spectra and band-integrated profiles for the obtained burst sample.

Fig. 1. Spectra and band-integrated sample intensity profiles for the 16 bursts from GMRT observations. Bursts G01–G06 were detected on 2021 June 01, the rest on 2022 February 06–07. All bursts are dedispersed to DM = 414.73 pc cm−3 determined from burst G05. For plotting, the spectra were normalized by the mean and standard deviation in the off-burst region in each sub-band and the colors were saturated at ±8σ for G05 and at ±4σ for the other bursts. All bursts from 2021 except for G05 are plotted with a time resolution of 5.24 ms and a frequency resolution of 22–25 MHz. Burst G05 is shown on a finer 2.62 ms/391 kHz spectro-temporal grid. For the 2022 bursts, the time resolution is 3.3–5.9 ms, and the frequency resolution is 24 MHz.

Table 1 lists the burst properties. Burst fluences were determined by integrating the emission over the 100 ms window centered at the burst peak. The equivalent width was calculated by dividing the fluence by the peak flux density.

All bursts except G05 are faint and little more can be inferred from their spectro-temporal shapes than a hint of downward-drifting trombone features, when dedispersed to the DM obtained from G05 (see Section 4.1). Overall, the peak flux densities, fluences and equivalent widths are comparable to the FRB 20201124A bursts recorded earlier with GMRT, in 2021 April (Marthi et al. 2022).

3. ALERT observations and analysis

Apertif is a phased-array front-end system installed on 12 of the 14 WSRT dishes⁶ (Adams & van Leeuwen 2019; van Cappellen et al. 2022). Apertif observed FRB 20201124A as part of the scheduled visits of repeater fields in the ALERT survey (Oostrum et al. 2020; van Leeuwen et al. 2023), on 2021 July 03 & 04, and on 2022 February 01, 05, & 06. The last two observations were coordinated to overlap with GMRT (see Sect. 3.2). All sessions lasted for three hours except for a 2.4-h-long session on February 05.

Apertif consists of phase array feeds on a multi-element interferometer, that form a hierarchical system of beams. During the FRB 20201124A observations, the central compound beam CB00 was pointed at the J2000 sky coordinates RA = 05^h08^m03.5^s, Dec = +26° 03′38.4″ (37.8″ for July 2021 sessions). These coordinates are close to the source position at RA = 05^h08^m03.5073^s ± 4.7 mas, Dec = +26° 03′38.5032″ ± 3.9 mas as reported by Nimmo et al. (2022). The offset lies well within the Apertif localization limits (Oostrum 2020).

Total intensity samples were recorded at a time resolution of 81.92 μs, and with 1536 channels of 195.312 kHz for a sample bandwidth of 300 MHz centered at 1369.6 MHz. As part of the standard real-time FRB search, the data from all 40 compound beams was independently searched for FRBs, using the AMBER⁷ search code (the Apertif Monitor for Bursts Encountered in Real-Time; Sclocco et al. 2016) in the DARC⁸ pipeline (the Data Analysis of Real-time Candidates; Oostrum 2021) on the Apertif Radio Transient System (ARTS; van Leeuwen 2014). Real-time candidate selection was carried out by a neural network, as described in Connor & van Leeuwen (2018).

The DARC pipeline did not find any candidates in the July sessions. This is in line with the reports by Main et al. (2021) on targeted observations with both uGMRT and the Effelsberg Telescope over 2021 June–Aug. Mao et al. (2022) also did not find any bursts down to 4 Jy ms during a 104-h observing run with Nanshan 26 m Radio Telescope in 2021 June–July.

During the 2022 February ALERT run, two bright FRBs were detected on the 1st, at DM = 410 pc cm⁻³. Initially, the bursts were found in CB17, which partially overlaps CB00 (Fig. 3 in van Leeuwen et al. 2023). Fluences and dispersion measures for these bursts were previously reported in Atri et al. (2022). Subsequent reanalysis showed the pulses were not detected in CB00 because of the residual RFI. The next session, on the February 05, yielded one more burst. No bursts were detected on the February 06, despite comparable system parameters and RFI.

3.1. Deep search

We subsequently performed an offline search over a finer grid of trial DMs and matched-filter boxcar widths than is possible in the real-time search. First, the filterbank files from CB00 were thoroughly cleaned of RFI using the iqrm⁹ software implementation of the Inter-Quartile Range Mitigation outlier detection algorithm (Morello et al. 2022). After cleaning with iqrm, we applied rficlean¹⁰, which operates in the Fourier domain (Maan et al. 2021). We employed the default threshold values, but omitted the excision of spiky RFI in the band-averaged time series.

Pulse candidates were selected using TransientX¹¹ (Men & Barr 2024). We searched for FRBs over DMs ranging from 360 to 460 pc cm⁻³ with a step of 0.25 pc cm⁻³. The search was performed on time series at the original time resolution. Boxcar filter widths ranged from 0.18 ms (slightly over two samples) to 200 ms. Each session yielded about 500 potential FRBs with integrated S/N ≥ 7, which were inspected visually. Most of the candidates were residual RFI, displaying sharp, narrow-band positive and negative jumps in the frequency-resolved signal. Only 24 pulses possessed the FRB-like properties we defined, being relatively broadband and displaying smooth variation of signal strength with frequency. Three brighter candidates matched earlier DARC detections.

Table 2 lists the times of arrival, best DM, integrated S/N and boxcar filter width for all 24 candidates. The majority (83%) have DMs larger than 411 pc cm⁻³, the DM of the brightest pulses from FRB 20201124A. However, transientX optimizes DM to maximize the S/N and for some bright pulses this clearly compromises the intrinsic spectro-temporal structure.

About 70% of all candidates come from one observing session, 2022 February 01. The 2021 July sessions and session 2022 February 06 yield two candidates each, and session 2022 February 05 resulted in only one (relatively bright) burst. It is possible that some of our faint candidates are due to chance noise fluctuations. In order to estimate the rate of occurrence of such noise candidates, we performed the same search in DM range between 660 and 760 pc cm⁻³, leaving all other parameters intact. In this manner three candidates were visually filtered from about 500 candidates per session. Two candidates were detected in session 2021 July 04 and one in 2022 February 05, their integrated S/N values were ≤8.1, widths ranged from 0.3 to 6.3 ms, and on the diagnostic plots the spectra looked indistinguishable from the spectra of faint FRB candidates from Table 2.

Table 2 lists 13 pulses with 7 < S/N ≤ 8.1, significantly more than the ∼3 that would have been expected by chance detections following the test described above. This likely means that some of those candidates were emitted by FRB 20201124A, although we cannot tell which ones exactly. Their faintness precludes any meaningful analysis, thus we do not further include them in the sample. The spectra of the remaining ten bursts were next computed using dspsr and are shown in Fig. 2.

Fig. 2. Spectra (bottom subpanels) and intensity profiles, integrated over the frequency band (top subpanels), of the ten bursts discovered within the ALERT survey. For plotting, the spectra were normalized by the mean and the standard deviation in the off-burst region in each sub-band, and the colors were saturated at ±10σ for W01, W07, and W10, and at ±4σ for the other bursts. Bursts W01 and W10 are plotted with the original tres = 81.9 μs, burst W07 with 0.33 ms and the rest with 2.62 ms. The frequency resolution is 25 MHz for all bursts except W07 and W10 (6.25 MHz), and W01 (0.48 MHz). The time in ms is zeroed on the burst peak. Spectra at frequencies above 1450 MHz are not shown, since the signal there is mostly corrupted by RFI.

Calibration of the ALERT FRBs is performed with the help of drift scan observations of the bright quasars 3C147, 3C286, and 3C48 at the beginning and the end of each observing run (cf. Connor et al. 2020; Pastor-Marazuela et al. 2023). We used the calibrator observation closest to a given FRB 20201124A session and estimated the System Equivalent Flux Density (SEFD) using the known quasar flux (Perley & Butler 2017). For 2022 February 01, the SEFD was 94 Jy and for the last two sessions it was 82 Jy. There is a slight (10%) variation within the band which was ignored. However, we took into account the excised parts of the band. Scaling the SEFD according to the radiometer equation resulted in about 170 or 200 times smaller SEFD for the band-integrated signal at the original time resolution (the denominator in Eq. (A1.16) of Lorimer & Kramer 2005), depending on the number of excised channels. Following Pastor-Marazuela et al. (2021) we assume 20% errors on the flux density values.

The quasar driftscan observations in the end of the 2022 June–July observing run failed, but test observations of pulsars immediately before and after FRB 20201124A observations did not indicate any malfunction. Taking a typical SEFD of 85 Jy as derived in van Leeuwen et al. (2023), the fluence $F=\mathrm{S}/\mathrm{N}\times \mathrm{SEFD}\sqrt{w_{\sec }}/\sqrt{n_{\mathrm{pol}}\mathrm{BW}}$ of the faint (i.e., below the adopted S/N threshold) bursts from Table 2 ranges from 0.7 to 4 Jy ms, which is comparable to the limits by Mao et al. (2022) and larger than the 0.02 Jy ms limits for 5 ms pulses after the emission quenching as reported by Xu et al. (2022). Still, we believe that the small excess of burst candidates detected near the plausible source DM (four around 410 pc cm⁻³ versus one candidate around the incorrect DM of 610 pc cm⁻³) does not provide compelling evidence for the detection of faint FRBs between S21 and F21 activity windows. More robust estimates of the chance probabilities of such detection are beyond the scope of this work.

Among the ALERT bursts, W07 stands out because of its complex structure, appearing to consist of two groups of pulses with separations comparable to the duration of the groups themselves, clearly visible in Fig. 2. The burst was actually detected as two separate events by transientX, but in what follows we analyze it as one cluster-burst, following the convention of Zhou et al. (2022), who, based on the waiting time distribution of the emission peaks from Xu et al. (2022), define such a “cluster-burst” as a collection of emission peaks with a separation less than 400 ms, without signs of bridge emission between them.

The ALERT rate of 3 bursts per hour is seemingly smaller than the 5.6–45.8 h⁻¹ reported by Xu et al. (2022). However, taking into account only those FAST pulses which satisfy the width-dependent fluence threshold based on an integrated ALERT S/N of 8.1, $F=0.94\sqrt{w_{\mathrm{ms}}}$ Jy ms, we find that the FAST rate was close to the ALERT values in the beginning of the FAST observing campaign and extrapolates to 10–14 bursts per hour for the FAST sensitivity limits.

On 2022 February 01, no pulses other than W07 appear clustered. Session 2022 February 05 yielded only one relatively bright burst, 11 minutes into the observation, but no other bursts, even faint ones, were detected later. Xu et al. (2022) report at least one instance of a change in rate by a factor 3 (5 on their full sample) between daily sessions. The absence of emission after 2022 February 5 is unlikely to be the end of this activity window, since Takefuji et al. (2022) observed a burst at 2.3 GHz on 2022 February 18. After that no other detections were reported.

3.2. Simultaneous observations by GMRT and ALERT

Apertif and the GMRT were co-pointing on 2021 February 05 from 18:20–19:30 UT, and on February 06 from 16:30–19:30 UT. During this time there is one burst detection, G07 (Table 1). There is no evidence for this same burst in the Apertif data. We conclude the burst emission is band limited, and does not extend from 650 MHz up to 1.4 GHz, similar to the behavior we have found earlier in FRB 20180916B (Pastor-Marazuela et al. 2021).

4. Burst analysis

4.1. Dispersion measure

Dispersion is an important characteristic of FRBs as it measures the integrated electron density on the line of sight (LOS) provided that the effect of dispersion can be separated from any intrinsic spectral shape of the burst. Variations of the DM may indicate a complex and dynamic circum-burst environment, since rapid DM changes are not expected at the galactic and inter-galactic level. When combined with the RM, the DM places limits on the magnetic field strengths encountered by the bursts (e.g., Xu et al. 2022; Lu et al. 2023).

As a quantity, the DM characterizes the magnitude of the pulse delay – a delay that is inversely dependent on the square of the observing frequency. In the absence of any a priori knowledge about the intrinsic spectral shape of a burst, the DM can be estimated by maximizing the coherent power in the burst, across the observing bandwidth (Seymour et al. 2019). This technique is employed by DM_phase¹², which operates as follows. During dedispersion, the time series from each particular frequency channel is shifted to counteract the expected dispersive delay for this frequency. For a Fourier transform along the time axis, this operation corresponds to a multiplication by exp(iϕ), with ϕ = ω C DM/ν², where ω signifies the Fourier frequency. The phase of the Fourier transform ϕ is integrated over ω and ν, and the resulting dependence of integrated coherent power on the trial DM values is examined for peaks. The error of the DM measurement is assumed to be the error of peak position determination.

We have utilized DM_phase with an automatic cutoff along the ω axis to measure the DM for all our bursts that have an integrated S/N > 50. The GMRT sample yielded one such burst, with a single component and slightly asymmetrical pulse shape. For this burst, the DM was measured to be DM = 414.73 ± 0.48 pc cm⁻³.

The WSRT sample contained three sufficiently bright bursts, all composed of a few distinct components. The DM measured here was lower, 410 pc cm⁻³, with a characteristic error of 0.3 pc cm⁻³ (Table 3). The initially reported measurement of DM = 410.9 ± 0.2 pc cm⁻³ for bursts W01 and W07b (Atri et al. 2022) was based on the less precise method of straightening the pulse structures visually. Both the GMRT and WSRT measurements agree with the distribution of DMs within the respective activity windows (Xu et al. 2022; Kirsten et al. 2024).

The structure-maximizing method of DM estimation is a de-facto standard in the FRB field at the moment. However, it does not provide unambiguous measurements for all bursts. For some FRBs, optimal DMs determined from the upper/lower halves of the observing band are inconsistent with each other, and burst spectral features cannot be aligned across the whole band (see, e.g., Platts et al. 2021; Kumar et al. 2022; Zhou et al. 2022).

Below we demonstrate that the existence of a DM value that maximizes burst coherent structure does not mean that this DM can be readily used to measure the integrated electron density along the path the emission traveled. Depending on the burst spectral shape, the structure-maximizing DM may be biased by fine structure buried in noise. We illustrate this using the comprehensive study of FRB 20201124A bursts by Xu et al. (2022). The authors compiled a set of DM values measured with DM_phase for bursts with integrated S/N of 20–3000 and a variety of time-frequency profiles. The authors rule out secular DM trends on the level of 2.9 pc cm⁻³ per two months, but record a large spread of DMs (with variations on the order of a few pc cm⁻³) within individual 2 h sessions, sometimes on a timescale of less than a minute or, remarkably, even less than a second.

There is a correspondence between the reported DMs and the integrated S/Ns of the bursts in the Xu et al. dataset: the DMs of the brightest bursts are listed as lower (around 411 pc cm⁻³), whereas the reported DMs of fainter bursts are spread between 409 and 424 pc cm⁻³ (Fig. 3). Kumar et al. (2022) and Lanman et al. (2022) report a similar spread of DM values for comparable S/N values¹³.

Fig. 3. Compilation of DM measurements from previous studies: Xu et al. (2022, X21), Kumar et al. (2022, K22), Lanman et al. (2022, L22), and this work, plotted as a function of integrated S/N.

The questions are whether this correlation signifies a causative relationship with the intrinsic brightness of the burst and/or whether it is related to the method or to instrument noise. To investigate the effect of noise on the performance of DM_phase, we utilized a nine-burst sample of burst spectra that were made public by Xu et al. at the time of manuscript preparation¹⁴. We take burst #377 from Xu et al. (2022) as an example, and discuss other bursts afterward.

The observed behavior of burst #377 is shown in the top left panel of Fig. 4. It has three partly merged components and an overall drift in frequency. The spectrum of the ω-resolved power versus the trial DM shows different zones corresponding to these scales. If we now add Gaussian noise to each frequency channel of the data (pre-normalized by the standard deviation of noise in the off-pulse region), the fine temporal structure is washed out ever more, and above a certain amount of added noise only the low-ω feature remains. Above this edge in additional noise (i.e., below this boundary in terms of S/N), DM_phase aligns the entire spectrum, converging on a DM about 10 pc cm⁻³ higher that the original value. The step is clearly visible when comparing the middle column in Fig. 4 to the right column.

Fig. 4. Influence of noise on the DM estimate from DM_phase. Burst #377 from Fig. 6 in the extended data in Xu et al. (2022) was normalized in the off-burst window in each channel (such that the standard deviation of the noise std = 1), and Gaussian noise was added in each channel. Each simulation was then run through DM_phase with an automatic spectrum cut. Top row: Example of spectra before and after DM correction for three different values of noise std. Middle row: Corresponding DM_phase spectra. Bottom row:ΔDM vs. a set of noise realizations, six simulation runs per each noise std. For clarity, some jitter along the horizontal axis was added to each set of points.

We note that even without such additional noise the DM derived with DM_phase depends on the choice of the frequency and time averaging, as well as the channel normalization method. The resulting spread of DM measurements is a few times larger than the estimated DM_phase uncertainties. The same level of discrepancy is observed between DMs from Xu et al. (2022) and our measurements. Also, on our normalized data we measure integrated S/N of about 1.5 times larger than reported by Xu et al. (2022). This could be partly attributed to normalization or, partly, be due to a finer grid of trial widths we used for the integrated S/N calculation. We also note that the some of the information in the header of burst spectra that accompany Xu et al. (2022) (e.g., DM and cardinal burst number) deviates from the corresponding entries in their data table.

How much the DM varies with pulse S/N, as well as the character of this variation (e.g., with or without a step), is determined by the temporal-spectral shape of the burst. The least amount of variation, less than a few times the DM_phase-reported error, was recorded for bursts with distinct, widely separated components (e.g., bursts #779, #1377, and #1398 from Extended Data Figure 6 in Xu et al. 2022). Mean while, bursts with an overall drift in frequency (e.g., #377, #460) exhibited steps of 1 − 10 pc cm⁻³ even at integrated S/Ns as large as 100. Thus, for the majority of the burst population in Fig. 3 the DMs are likely overestimated if unresolved drift in frequency was present.

FAST provides by far the largest and brightest sample of bursts, since for other telescopes the S/N is usually smaller, meaning that DM overestimation is widespread. As this is a matter of S/N, the same effect will occur for very bright bursts observed with less sensitive telescopes. It is worth noting that a DM measurement bias introduced by drifting of separate burst components (that can go unrecognized if the burst does not have sharp sub-bursts) was investigated on a large sample of bursts from another prolific FRB repeater, FRB 20121102A, by Jahns et al. (2023). In that work, the measured tilt of each individual sub-burst was directly converted to the resulting change in the fitted DM. This equivalent DM can be used as an upper limit on structure-maximized DM bias, since the algorithm for the latter is typically using information from a collection of sub-bursts, which can have different drift rates. The upper limits on this DM bias have the same range as in our work, namely ≤12 pc cm⁻³. While the authors speculate that the actual bias in structure-maximized DM estimates is likely to be less than 1 pc cm⁻³, no direct studies have been performed.

Assuming that the DM excess of Xu et al. is indeed due to absorption of the trombone effect, we can calculate the absorbed drift rate from an extra dispersion delay between the edges of the band. For a DM ≥ 412 pc cm⁻³, the frequency drifts are between 25 and 225 MHz/ms, comparable to the values found by Zhou et al. (2022), who determine their DM values from a sample of bursts with sharp, separate components.

If the true DM at the time of the GMRT observations was 411 pc cm⁻³, then the DM measured from burst G05 would imply a drift rate of 8 MHz/ms (Fig. 5). It is useful to compare this value to the study of Marthi et al. (2022), who calculated drift rates for a sample of 48 FRBs recorded during a single session with the GMRT. In that study, a DM of 410.78 ± 0.54 pc cm⁻³ was determined using a singular-value decomposition of the spectrum of the brightest burst, exhibiting a sharp burst rise and a few partially merged components. All bursts in the study of Marthi et al. (2022) displayed trombone drift, with rates between 0.75 and 20 MHz/ms. We thus conclude that the DM determined for G05 has some drift absorbed in it.

Fig. 5. Burst G05 dedispersed with structure-maximizing DM = 414.73 pc cm−3 (left) and with DM = 411 pc cm−3 (right), suggested by the brightest bursts observed by FAST.

For the F21 activity window, Zhou et al. (2022) measure an average DM of 412.4(3)–411.6(3) pc cm⁻³, using only those pulses that show sharp edges or well-separated components. Individual measurements are spread within ±2 pc cm⁻³, larger than the reported errors. Kirsten et al. (2024) measure DMs consistent with Zhou et al. (2022) for that day. Statistically, their DMs for the W22 activity window are no different from our measurements. Overall we conclude that reliably detecting any secular DM trend on the level of 1 pc cm⁻³ between activity windows requires a more robust method of dealing with the influence of burst structure than is currently used.

4.2. (Quasi-)periodicity

Some FRBs consist of multiple components arranged in a seemingly regular fashion. A reliable detection of such (quasi-)periodicity may have interesting implications for FRB emission theories. For example, such periodicity may be a direct manifestation of the relatively fast spin period of an emitting compact object. Or, it may reflect the features of spark generation or non-stationary plasma flow in the neutron star magnetosphere (e.g., Mitra et al. 2015). Finally, quasi-periodicity of subpulse components appears naturally when FRB generation is driven by crust motion (Wadiasingh & Chirenti 2020).

So far, the most reliable detection of periodicity comes from a one-off FRB. Chime/FRB Collaboration (2022) detected a 217 ms periodicity in the nine components of FRB 20191221A. This 6.5σ periodicity is consistent with beamed emission from a neutron star rotating at ∼5 Hz. Other detections are close to the 3σ level and exhibit quasi-periods ranging from submilliseconds to tens of milliseconds (Pastor-Marazuela et al. 2023; Chime/FRB Collaboration 2022).

FRB 20201124A presents an interesting example of a repeating FRB source characterized by apparent periodicity in some of its sub-bursts. Nevertheless, timing analysis of the sub-burst TOAs derived from 53 bursts within the F21 activity window, as conducted by Niu et al. (2022), did not reveal a single period with possible harmonics, but unveiled a broad distribution of periods spanning the range of 1 to 10 ms, with an extended tail reaching up to 50 ms. All identified periods exhibited significances not exceeding 3.9σ of the normal distribution.

In our sample, only bursts W01, W07, and W10 exhibit multiple components with potentially periodic spacing. To test for this subsecond periodicity we computed the auto-correlation function (ACF) of W01, W10 and the two parts of W07 (Fig. 6), but no prominent periodicity was found. The quasi-periods corresponding to the peaks on ACF are in good agreement with the distribution obtained by Niu et al. (2022) on a larger sample of bursts.

Fig. 6. Left column: Bursts W01, W07a, W07b, and W10 with potential quasiperiodic components. Right: Autocorrelation function from all four bursts.

We complemented the ACFs with a timing analysis, which may be more sensitive to short periodicities potentially buried in the zero-lag peak of the ACF. In the timing analysis, periodicity is searched for in the sample of sub-burst TOAs, and the accuracy of the timing analysis is greatly influenced by the precision of the TOA measurement. If sub-bursts have complex shape and are closely spaced, it becomes difficult to determine which parts of the time series belong to different sub-bursts, and which reflect the intrinsic shape of an individual sub-burst. An example of this can be seen in the sub-bursts of FRB W07 on the upper subplot of Fig. 7 around t = −10 ms.

Fig. 7. Top: TOAs of individual sub-bursts of FRB W07 (red dots) identified as peaks of the signal smoothed with Gaussian kernel (thick line). The vertical lines show periodicity with best quasi-period P = 1.44 ms. Middle: Phase residuals with respect to model based on the best-fit quasi-period P. Bottom:Z2 statistic on the range of trial periods.

For the subsequent analysis, we defined each TOA as the time stamp of the peak on the signal smoothed with a Gaussian kernel with a five-sample standard deviation. The peaks were located using scipy.signal.find_peaks.

For bursts with prominent, ostensible quasi-periodicity it is logical to assume that the period is close to sub-burst separation or its integer multiplicative (e.g., FRB 20191221A in Chime/FRB Collaboration 2022, where sub-bursts arrive every one or two periods). The situation becomes more complicated when bursts are less frequent, as in FRB 20201124A. Searching for periods around the minimum sub-burst component separation then leads to a large variation of these periods from one burst to another, as was demonstrated by Niu et al. (2022).

In what follows we thus do not assume that the minimum separation between sub-bursts strictly constrains the possible periodicity. Instead, we performed a uniform search over a grid of trial periods using the Z² statistic:

$$\begin{array}{c}\hfill Z^2=\frac{2}{N}[{({\displaystyle \sum _{i=0}^{N-1}}cos{\varphi }_i)}^2+{({\displaystyle \sum _{i=0}^{N-1}}sin{\varphi }_i)}^2].\end{array}$$(1)

Here ϕ_i = 2πt_i/P is the phase of ith sub-burst component computed with trial period P. The statistic Z² was maximized over the range of P between 0.5 and 6 ms with an increment of 0.01 ms. We estimate the significance by repeating the analysis on 10⁵ samples of random TOAs, keeping the first and last t_i fixed, and drawing the remaining TOAs from a uniform distribution spanning (t₀,  t_N − 1).

Applying this analysis to the 11-component burst from Niu et al. (2022), we find an optimal P of 3.07 ms, Z² = 11.48. In the batch simulated of sets containing 11 random TOAs, 99.75% of sets had a maximum Z² < 11.48 for P = 3.07 ms. This equates to a significance for the observed burst equivalent to 2.84σ for a normal distribution. The results are close to those obtained by Niu et al.: P = 3.06 ms with significance of 3.3σ¹⁵. However, since we searched over a grid of P, the number of trials should be taken into account. Comparing the maximum Z² score of the real data to the pool of maximum Z² scores lowers the significance to 0.3σ. Thus, for these weak signals the presence of a priori constraints on P is crucial for obtaining a significant result.

Among the bursts in our sample, W07 exhibited the strongest evidence for quasi-periodicity. This candidate had a period of 1.44 ms, and a significance of 3.3σ directly, and 0.68σ after correction for the number of trials (Fig. 7). For Bursts W01, W10, and the two sub-burst groups of W07, no discernible periodicity with single-trial significances greater than 2σ was identified. It appears the sub-bursts lack prominent quasi-periodicity, at least when not considering sub-burst shape appropriately or in the absence of motivated constraints on P.

Under the assumption that sub-bursts appear almost every period, and that P can change from one burst to another, we examined the distribution of P from Niu et al. (2022) in order to investigate whether the frequencies 1/P were clustered around specific harmonics, as predicted by the crust motion and low-twist theory (Wadiasingh & Chirenti 2020). We did this by obtaining the eigenmode l-number from the formula $\nu =0.5{\nu }_0\sqrt{(l-1)(l+2)}$ over a range of ν₀ trials. For each trial ν₀ we examined the deviations of the obtained values for l from their respective nearest integer counterparts. These distributions appeared to be uniform and indistinguishable from the same distributions computed on sets of uniformly distributed random period values. Although no clustering was hence found, we cannot disprove crust motion theory this way. Since the eigenfrequency depends also on the strength of magnetic field at the location and time of the crust motion, it may vary from one burst to another (Wadiasingh & Chirenti 2020).

4.3. Burst fluences

Tables 1 and 3 list the peak flux densities and fluences of the bursts recorded in this study. Prior to measuring the peak flux densities, the band-integrated time series were averaged by several time bins (see the caption of Fig. 2). For fainter bursts, the peak flux density depends on how much of such averaging was performed.

The fluence was estimated by summing the signal in a fixed-size window around the burst peak. For all bursts except W07, this window size comprised 60 ms. For W07 we used a larger window, from −120 to 20 ms, to encompass all components. The fluence of W01 is in agreement with the CB17 detection we reported in Atri et al. (2022). For W07, we measure a 3 times larger fluence in the larger window, as we are now including the component group W07a, in contrast to the earlier reported results from the standard pipeline.

4.3.1. Fluence distributions

Statistical distributions of burst fluences provide valuable information about the FRB emission mechanism, and about propagation in the circum-burst environment (Cui et al. 2021; Xiao et al. 2024). Owing to the large number of bursts observed, FRB 20201124A has the potential to offer one of the best-measured distributions of burst fluences among repeating FRB sources. In practice, however, the distributions from different observations exhibit little agreement with each other (Fig. 8). The apparent mismatch is not entirely caused by the difference between the mean event rates at different epochs, but the shape of the distributions themselves seem to have intrinsic changes.

Fig. 8. Upper row: FRB fluence vs. MJD of arrival for three activity windows. The markers and colors correspond to the legend in the lower panel (see inset), which shows rate survival function vs. burst fluence. FRB information was collected from Lanman et al. (2022, Lan22), Xu et al. (2022, Xu21), Nimmo et al. (2022, Nim22), Hilmarsson et al. (2021, Hil21), Kumar et al. (2022, Kum22), Marthi et al. (2022, Mar22), Zhang et al. (2022, Zha22), and Kirsten et al. (2024, Kir24). Where possible, fluence distributions obtained from observations on different days and with different telescopes were computed separately.

A series of dedicated observations by the FAST telescope demonstrated that the FRB 20201124A burst rate is highly variable (by up to two orders of magnitude), both within and between activity windows (Xu et al. 2022; Zhang et al. 2022). The fluence distribution is bimodal (Zhang et al. 2022), but unfortunately, the pulse rate is insufficient for exploring any changes in the distribution shape on the timescale of the rate change, although some studies have been performed (Sang & Lin 2023).

The fluences of bursts from WSRT observations fall within the range of fluences reported by other authors. Fainter bursts likely belong to the fainter component of the bimodal fluence distribution seen in the FAST observations Zhang et al. (2022). No contemporaneous observations were conducted during W22 activity window, however, and the position of the fainter component is known to shift from one activity window to another. Fluence distributions for W22 activity window recorded with WSRT and S21 activity window recorded with GMRT are flatter than normal, although similarly flat distributions have been recorded previously.

4.3.2. Fitting methods

Generally, a cumulative distribution of burst fluences can be approximated with a power-law (PL) distribution, with a possible flattening at lower fluences, either intrinsic or due to sample incompleteness close to observational sensitivity limit. PL fits are easy to perform, have only few free parameters, and allow for quick comparisons with theoretical models (e.g., Wadiasingh et al. 2020). There are, however, limitations that should be kept in mind as these can bias the physical interpretation of the results.

One of the caveats concerns the fitting method. Historically, fluence distributions of individual pulses of pulsar radio emission have been approximated with PL functions by performing a least-squares linear fit to the survival function on a log–log scale (hereafter the “graphical method”), and this practice is still sometimes used for the FRB fluence distribution (Popov & Stappers 2007; Bilous et al. 2022; Pastor-Marazuela et al. 2021; Kirsten et al. 2024). Despite the ostensible transparency of the graphical fitting method, the least-squares minimization does not provide an accurate and unbiased estimate of the PL parameters – even if the number of the sample is quite large by FRB standards (∼100 pulses, see Goldstein et al. 2004; Hoogenboom et al. 2006). A Maximum Likelihood Estimator (MLE) provides a much more accurate estimate of the PL index. In Appendix A we provide a detailed comparison of the two methods.

As an illustration, we refit the fluence distributions recently obtained by Kirsten et al. (2024) using both methods: graphical and MLE. The outcome of the original fits prompted the authors to speculate that there may be two separate populations of pulses emitted by FRB 20201124A, with the more energetic population having substantially flatter distribution. Unlike the authors, we fit fluences instead of spectral energy densities, however this does not affect the shape of the distribution. While we use all available Onsala measurements, we necessarily exclude the Stockert observations from fitting with both methods, as the large difference in exposure time and sensitivity makes combining the datasets for subsequent MLE fitting very difficult. As in Kirsten et al. (2024), the optimal minimum fluence value for both the FAST and Onsala distributions was determined with powerlaw¹⁶. This software finds the optimal minimum fluence by minimizing the Kolmogorov-Smirnov distance between the data and the fit on a sample of fits that start from each unique value in the dataset. In the same way as Kirsten et al. (2024), the obtained minimum fluence value was used to filter fluences for fitting with the graphical method.

Despite the above-mentioned differences in approach, our implementation of the graphical method reproduces the values of the measured power-law index (PLI) α_m and its corresponding bootstrapping errors as reported by Kirsten et al.: α_m(Onsala) = − 0.50 ± 0.06 and α_m(FAST) = − 2.11 ± 0.09¹⁷.

The MLE method, however, results in α_m(Onsala) = − 0.63 ± 0.26 and α_m(FAST) = − 2.15 ± 0.19. Clearly, the MLE determines the errors to the fit to be larger than the graphical method indicates. In order to determine the statistical significance of the difference between the Onsala and FAST measurements we compared this difference to the respective measurement errors added in quadrature. This method assumes that both measured values are drawn from normal distributions, with a mean equal to the true value and a variance defined by measurement error. Therefore the difference between measurements is also normally distributed, with a mean of 0 and a variance equal to the sum of variances of individual measurements. In our case, for the graphical method the ratio of the measurement difference to its supposed standard deviation is 13.7, corresponding to practically zero chance probability. For MLE, the ratio is 4.7, corresponding to a chance probability of 1.3 × 10⁻⁶. This indicates that the difference in PLI between Onsala and FAST is much less significant when measured with the MLE method.

However, the distribution of α_m around true value α_t is not Gaussian: α_m|α_t obeys a gamma distribution, which deviates from a normal distribution for small sample sizes, N ≲ 100. For a uniform prior on α_t this results in the posterior distribution for α_t|α_m being skewed toward steeper PLI for small N. In Fig. 9 we show the posterior distributions for both the FAST and Onsala samples constructed from the fitted α_m on a grid of trial α_t. The posterior α_t is well described by the gamma distribution from James et al. (2019) modified for the unbiased estimate α′ (in their notation) by taking M = N − 1:

Fig. 9. Left: Survival functions of the fluence distributions from the FAST and Onsala observations during the S21 activity window, together with MLE PL fits. The vertical lines mark the optimal minimum fluence as found by the powerlaw package, after modification according to Eqs. (A.2)–(A.3) in Appendix A. Top right: Posterior distributions of αt|αm. The gray lines indicate gamma distributions from Eq. (2). Bottom right: Probability of αt(Onsala)−αt(FAST) being smaller than the threshold value.

$$\begin{array}{c}\hfill p({\alpha }_{\mathrm{t}}|{\alpha }_{\mathrm{m}})\sim {(-{\alpha }_{\mathrm{t}})}^{N-1}{\left(\frac{N-1}{-{\alpha }_{\mathrm{m}}}\right)}^{N}exp(-\frac{{\alpha }_{\mathrm{m}}(N-1)}{{\alpha }_{\mathrm{t}}}).\end{array}$$(2)

Having two posterior distributions, one may calculate the probability that the difference between the α_t values for the FAST and Onsala samples is smaller than some threshold value (Fig. 9). These probabilities remain low for differences in PLI less than ∼0.5, indicating that a significant difference exists between α_t(FAST) and α_t(Onsala). This significance is, however, appreciably smaller than implied by graphical method. For example, p(Δα_t ≤ 0.5) = 5.3 × 10⁻⁴ graphically, but 3.2 × 10⁻² according to the MLE.

To summarize our findings up to here, accurate estimates of the PLI for small burst samples critically depend on two factors: first, using an unbiased PLI fitting method, and second, taking into account the skewness of the posterior probability distribution.

There is, however, one more caveat connected to sampling from flattened distributions. If the observed fluence distribution flattens at the low end, either due to instrument sensitivity, or intrinsically, then this flattening cannot always be recognized in small samples, and the measured power-law index (PLI) can be biased toward shallower values. The amount of this bias depends on the extent of flattening. As an example, we show the distribution of 10-pulse PLIs from the FAST S21 sample in Fig. 10. The sample was truncated at 5 Jy ms and the flattening signifies an intrinsic property of the pulses. The PLI distribution is skewed toward shallower values, which is expected because low-fluence pulses with shallower survival functions are more abundant. However, about 17% of these 10-pulse samples do not exhibit apparent flattening at the lower-fluence edge of survival function. Their survival functions are well-fitted with a single PL for the whole extent of the 10-sample distribution. A chance observation resulting in one such 10-burst realization may well lead to the erroneous conclusions that the distribution of burst fluences obeys a PL with a shallow index, and to a subsequent unfounded scientific interpretation.

Fig. 10. Survival functions of the fluence distribution from the S21 FAST bursts (Xu et al. 2022) with fluences larger than 0.5 Jy ms, together with an example survival function from a random ten-pulse subsample. The vertical lines indicate the optimal minimum fluence values from the PL fit. The example survival function does not show apparent flattening. Inset: Distribution of αm for all 104 ten-pulse subsamples (green) and those that do not show flattening (orange).

The exact amount of PLI bias toward shallower values depends on the sample size and on the underlying fluence distribution as observed with a specific instrument under a specific setup. Assumptions about intrinsic fluence distribution as well as models of instrumental limitations (see, e.g., Gardenier et al. 2019; Wang & van Leeuwen 2024, for such models) should both be taken into account while interpreting the measured PLIs.

4.3.3. Fitting results

In this section, we review the PLI measurements for the fluence survival functions of bursts recorded in our observations and reported in the literature. Table 4 summarizes the results of the MLE fitting procedure. Unlike Kirsten et al. (2024), we did not combine bursts recorded by different telescopes, and we did not correct for the unknown possible underestimation of CHIME fluences (Lanman et al. 2022). For the FAST sample and our WSRT observations, we combined bursts with separation smaller than 100 ms to facilitate comparison with Kirsten et al. (2024). Such proximity threshold lies within the short-recurrence component of the bimodal waiting time distribution and is smaller than the 400 ms sub-burst separation threshold in Zhou et al. (2022), which was derived from the trough between the two components of the waiting time distribution. In practice this means that for the small fraction of bursts the fluences of individual sub-burst clusters, as defined by Zhou et al. (2022) were counted separately. This, however, had only minor influence on the shape of fluence distributions.

In Table 4, all measured α_m values are within the range of −0.5 to −4, but with the caveat that nominal PLI values, even with their respective MLE fit uncertainties, should be treated with caution due to the apparent flattening of distributions at the lower-fluence end, which can be either instrumental or intrinsic. To illustrate this, we show the powerlaw fit for four individual sessions of the F21 FAST burst sample (Fig. 11). During the observing run, the rate of pulses increased more than tenfold, from 31 to 412 combined bursts per one-hour session. The overall shape of the fluence distribution did not change dramatically, yet for the first session, the PL fit yielded a flat α_m = −0.49 ± 0.09 without a low-fluence cutoff. For the subsequent sessions, the optimal fit excluded the low-fluence region, resulting in a much steeper PLI (e.g., α_m = −1.8 ± 0.2 for MJD 59485). It is possible that the shallow α_m for the first session is solely due to the aforementioned bias present in small samples drawn from a flattened distribution; however, intrinsic variability cannot be ruled out.

Fig. 11. Un-normalized survival function of the fluence distributions (red line) for the bursts recorded at the end of F21 activity window (Zhang et al. 2022). The gray rectangles mark the regions in the fluence–burst number parameter space that contain bursts excluded from the PL fit (see Section 4.3.2 for details). The colored rectangles contain bursts included in the fit; the color represents the steepness of the best-fit power-law index (PLI) ranging from light yellow (αm = −0.5) to dark blue (αm = −2.5). Light blue bands of progressively lighter shades in the background correspond to regions containing 10, 100, and 1000 bursts.

The fits for per-session fluence distributions of the FAST sample from the S21 activity window also exhibited a dependence on the number of bursts recorded per session. The number of bursts varied from 11 to 99 per session. The measured PLIs were more diverse for smaller burst samples, ranging from approximately −4 ± 2 to −0.7 ± 0.3¹⁸. Occasionaly, α_m reached much steeper values of about −10 ± 9. The fraction of power-law-obeying bursts ranged from 5% to 90%, with steeper indices corresponding to larger minimum fluences and a smaller number of bursts in the PL tail. The two most prolific sessions, with more than 80 FRBs each, yielded α_m values of −1.2 ± 0.3 and −1.7 ± 0.3, with 25% of bursts belonging to the PL tail.

All three burst subsamples in our study are small, with ten or fewer bursts each. The bursts from two GMRT sessions in the S21 and W22 activity windows resulted in dramatically different α_m values, −0.80 ± 0.40 and −3.59 ± 1.47. The steep PLI resembles similar steep values of the Spring FAST per-session fit when the burst sample size was small, and a relatively large minimum fluence was found for the PL tail. The WSRT burst sample is described by a shallow distribution with α_m = −0.71 ± 0.25 and no apparent low-fluence flattening. This PLI is close to the ones reported by Kirsten et al. (2024) and Hilmarsson et al. (2021), but extends to smaller fluences. Our α_m for the CHIME data is shallower than the MLE fit from Lanman et al. (2022), as they excluded the brightest burst due to uncertainties in fluence determination.

Fig. 12 provides a graphical representation of the PLI measurements from Table 4. Burst samples are sorted by the number of bursts in the PL tail. Overall, PLIs determined from burst samples smaller than about a hundred tend to have diverse values that are inconsistent between different studies. For the largest sets of pulses (696 and 1715 FRBs from the two FAST studies), only a small fraction (about 10%) of the brightest pulses follow a power-law distribution, with α_m close to −2 and minimum fluence differing by a factor of 3. Since these observations were performed under the same observing setup, the difference signifies intrinsic variability between activity windows.

Fig. 12. Un-normalized survival function of the fluence distributions for the bursts recorded in the previous studies and current work (see caption of Fig. 11 for the explanation of the legend and caption of Fig. 8 for the literature abbreviations). Here the burst samples are ordered by the number of bursts in the best-fitting PL tail. MJD refers to the start of respective observing campaigns.

Given the inevitable instrumental bias, possible intrinsic frequency evolution and temporal variability, constraining the PLI is a daunting task. One can also well ask the question whether PL approximation should be used in this case at all, since for both of the two samples of bursts that were largest and most sensitive, only 10% of bursts have fluences that obey a PL distribution.

4.4. Scintillation

In order to measure the decorrelation bandwidth, we obtained the spectra of the brightest bursts in our sample, G05, W01, W07A, W07b, and W10, by averaging their emission in time over the time bins where the signal is four times larger than the noise standard deviation. Note that for W07a, only the last, brightest component satisfies this criterion. For each ALERT burst, we removed the data above ∼1423 MHz, where RFI becomes strong in the observations. For the uGMRT burst G05, we fitted the burst spectrum to a Gaussian and only took the frequencies within two standard deviations of the center in order to have enough signal for the analysis; this is between 580 MHz and 709 MHz. Next we compute the ACFs of the spectra, removing the zero-lag frequency value, and fit the central peak of the ACF to a Lorentzian. The decorrelation bandwidth is often defined as the half-width at half-maximum of the ACF’s fitted Lorentzian (see, e.g., Lorimer & Kramer 2005, Section 4.2.2). For the ALERT bursts, we obtain an average decorrelation bandwidth (weighted by the inverse standard deviation of each measurement) Δν_sc = 1.1 ± 0.2 MHz at the central frequency 1336.5 MHz. For G05, we obtain Δν_sc = 0.15 ± 0.05 MHz at the central burst frequency, 644.5 MHz. Both Δν_sc values are much larger than the channel width of 195 and 48 kHz, respectively.

The frequency-dependent intensity variations produced by scintillation are expected to follow a power law evolution of the form Δν_sc = Aν^γ, with ν the frequency in GHz, A a constant that gives the decorrelation bandwidth in MHz at 1 GHz, and γ the scintillation index, expected to be γ = 4 for scintillation produced by a thin screen and γ = 4.4 for scintillation produced in an extended, turbulent medium. To measure the power- law index of the decorrelation bandwidth, we divide the bandwidth into several sub-bands, and we measure the half width at half maximum (HWHM) in each sub-band as described above. We divide the Apertif bursts into five sub-bands, and the uGMRT one into four. For the Apertif bursts, we compute the HWHM weighted average in each frequency sub-band, and then we fit all decorrelation bandwidths as a function of frequency to a power-law spectrum. We obtain A = 0.44 ± 0.03 MHz, and γ = 2.4 ± 0.2. The results from the scintillation analysis are presented in Fig. 13.

Fig. 13. Decorrelation bandwidth analysis. The top panel shows the ACF from 1220 MHz to 1453 MHz of the Apertif bursts W01 (orange), W07a (green), W07b (red), W10 (purple), and G05 (blue). The middle panel shows the ACF from 580 to 709 MHz of the uGMRT burst G05 in blue. The bottom panel shows how the decorrelation bandwidth evolves with frequency. The Apertif bandwidths were divided into five sub-bands, and the uGMRT bandwidth into four sub-bands. The gray dashed line with a shaded area shows the best-fitting power law, with an index α = 2.4 ± 0.2. As a reference, the black solid line shows an index α = 4 expected from a thin screen. The gray shaded region indicates decorrelation bandwidths that cannot be resolved by the uGMRT resolution.

Our measurements at both frequencies are based on a small number of bursts. Previous studies have shown that both the decorrelation bandwidth and γ vary substantially when measured on individual bursts within a relatively narrow bandwidth of 500 MHz centered around 1250 MHz: Xu et al. (2022) reported a mean γ of 4.9 with a standard deviation of 6.4. This variation may be at least partially influenced by the intrinsic spectral structure of the bursts. Inferring γ from the combined spectra of a few dozen bursts detected within an hour-long observing session resulted in a shallower γ of 3.0 ± 0.2 (Zhou et al. 2022).

Main et al. (2022) increased the robustness of γ determination by comparing decorrelation bandwidths at 700 and 1400 MHz. Their setup is similar to ours in terms of frequency coverage, but their sample of bursts is larger. They reported a power-law index of γ = 3.5 ± 0.1, which is also shallower than the one for the thin screen model. However, both in our measurements and those from Main et al. (2022), the decorrelation bandwidth at the lowest radio frequencies may be biased by insufficient frequency resolution, resulting in γ being biased toward shallower values.

The study of Main et al. (2022) utilized bursts from the S21 activity window, with low- and high-frequency observations conducted several days apart. Our measurements at 1.4 GHz come from the W22 activity window, 245 days after the S21 observations at 650 MHz. Wu et al. (2024) reported small annual variations in the observed decorrelation bandwidth attributed to the Earth’s movement with respect to a moderately anisotropic scattering screen located in the Milky Way (see also Main et al. 2022). This variation is much smaller than the scatter associated with measurements from individual sessions and therefore cannot be the main source of systematic uncertainty in our γ measurement. The available bulk of observational data does not show evidence for a secular trend in the decorrelation bandwidth or γ across three activity windows.

Main et al. (2022) detected enough nearby bursts from FRB 20201124A to measure its scintillation timescale of 13.3 ± 0.8 min in the L-band. Unfortunately, the most closely spaced bursts with sufficient S/N in our sample (W01 and W07) are separated by ∼110 min, which does not allow us to probe the scintillation timescale of the source. The correlation coefficient between W01 and W07a/b is close to zero, consistent with the previous measurement.

4.5. TOA statistics

The similarity in burst arrival statistics between FRBs and the high-energy short bursts from magnetars was one of the key observational facts put forward by Wadiasingh & Timokhin (2019, hereafter: WT19) in support of their crust motion and low-twist model of FRB generation. Both observed burst types feature a log-normal distribution of waiting times, and a log-uniform distribution of TOAs as measured from the outburst start time. WT19 studied the distribution of burst arrival times using a sample of 93 FRBs from FRB 20121102A recorded over the span of five hours. Some of the FAST observations of FRB 20201124A provide much higher a pulse rate (up to a few hundred pulses per hour), but the duration of the observing session is smaller, typically less than two hours.

Using the publicly available information about burst fluences and arrival times from Xu et al. (2022) and Zhang et al. (2022), we constructed the distribution of waiting times between successive bursts. As mentioned before, this distribution consists of two non-overlapping log-normal components. Based on a common definition of sub-bursts, we combined bursts which were separated by less than 400 ms into one, using an iterative routine. For each resulting cluster of sub-bursts, the TOA was taken to be the time of arrival of the first burst, and the cluster fluence to be the sum of the individual component fluences. Combined bursts from each observing day were analyzed separately.

Figure 14 shows the distribution of burst TOAs measured with respect to the session start for the most prolific observing session, on September 29 2021 (MJD 59486; see Table 4). The distribution is uniform with linear time bins and, unlike for FRB 20121102A (WT19), there is no linear correlation between the logarithms of the TOAs and the corresponding waiting times (the brown points in the top-right subplot). This discrepancy can, however, be explained if the observation took place some time after the outburst started, when the burst rate no longer changed significantly during a session.

Fig. 14. TOA and waiting time distribution for the most prolific FRB 20201124A observing session from Zhang et al. (2022) on MJD 59486 (upper row), and a series of log-uniform random variables (lower row) that simulate an FRB outburst. In our model only a small fraction of the outburst duration is covered by observations, long after the burst onset. The columns are: a) distribution of TOAs measured from the first burst arrival (brown) and from the hypothetical outburst start (violet), with the pale shade of violet showing the whole underlying mostly unobserved outburst; b) same distribution, but binned with logarithmic time bins; c) distribution of TOAs vs. waiting time between two bursts; d) distribution of waiting times. Bursts separated by less than 400 ms combined, producing a single TOA for the start of the combined burst.

To illustrate this, we simulated burst arrival times by generating 6 × 10⁴ log-uniformly distributed random variables t, using scipy.stats.loguniform. These random variables fell within the range of 1 < t < 10⁴ s with a probability density function of p(t)∼1/t. Approximately 1.5 hours after the start of the simulated outburst, the rate of occurrence for t remains relatively constant over the duration of the FAST session.

If TOAs are measured from the start of the simulated outburst, TOA = t − 5400 s, their statistics closely resemble the statistics of real bursts (Fig. 14). This suggests the possibility that the observed bursts are part of an outburst that began several hours before the observations. In this model, the variation in the session-to-session burst rate, which shows significant growth toward the end of the F21 activity window, could be attributed to different outbursts overlapping, or occurring in close succession. Unfortunately, the scope of this work does not permit a quantitative assessment of the probability of FAST missing the start of an outburst after observing FRB 20201124A for almost 60 days with < 10% duty cycle. Such an assessment would require an estimation of the tentative duration of an outburst and the frequency of their occurrence.

4.6. Frequency-dependent activity window

FRB 20201124A was extensively observed, by multiple telescopes, around the end of the S21 activity window. Figure 15 provides a summary of the time/frequency coverage of these observations, including the TOAs and fluences of the detected FRBs, plus fluence upper limits for any non-detections.

Fig. 15. Upper panel: Overview of daily detection statistics at the end of S21 activity window from all telescopes combined. A filled rectangle indicates that the source was detected, a hatched rectangle for non-detected, and a white rectangle signifies there were no observations that day. Lower frequencies are shown at the top, following Sect. 5.1. Middle panel: FRB fluences, together with upper limits for non-detections vs. the integer part of the observation MJD. The time region around the end of the S21 activity window is shown. The red triangles mark GMRT observations at 400 or 650 MHz (bursts G01–G06 from Fig. 1). The darker and lighter red shades show detections (pentagons) or upper limits for bursts in similar frequency regions from the observations by CHIME/FRB (Lanman et al. 2022) and NC (Trudu et al. 2022). The results of higher-frequency L-band observations by FAST (Xu et al. 2022) and NSRT (Mao et al. 2022) are shown in darker and lighter blue, with dots marking FAST detections. Detections by Torun/Onsala/Stockert (T/O/S) dishes from Kirsten et al. (2024) are shown as blue triangles (see text for details on fluence upper limits). Lower panel: Observation time coverage.

The majority of the FRB detections in the L-band were provided by sensitive FAST observations. These observations occurred at intervals of 1 or 3 days, with each 2-hour session yielding the detection of dozens of bursts with fluences above 0.053 Jy ms. Notably, there was an abrupt cessation of emission between MJD 59360 and 59363, after which no FRBs were recorded despite the unchanged observing setup and cadence.

These FAST observations were complemented by Kirsten et al. (2024), who conducted near-daily observations using several smaller telescopes. Their observations necessarily featured a significantly higher detection threshold, approximately 10 Jy ms. The last known pulse from the S21 session, detected by Kirsten et al., occurred on May 28, one day before the first non-detection by FAST, on a day when FAST was not conducting observations.

Shortly after this then yet unknown quenching, Mao et al. (2022) executed an extensive targeted search for bursts from FRB 20201124A in the L-band using the Nanshan 26 m radio telescope (NSRT). The authors determined a minimum detectable fluence of 4 Jy ms. On June 02 and 03, when FAST was not observing, the NSRT sessions were significantly longer than those generally used at FAST. If the source had persisted as active as before, several bright pulses should likely have been detected on these dates – but none were. On June 07, both NSRT and FAST observed FRB 20201124A, a few hours apart, with neither telescope detecting any bursts.

At the lower frequencies centered around 400–600 MHz, the majority of observations are provided by CHIME/FRB. That transit instrument records FRB 20201124A for 3.13 minutes virtually every day in the frequency range of 400–800 MHz. Lanman et al. (2022) estimate the burst rate after March 20 to be between 0.9–2 × 100 day⁻¹, and the bursts exhibit Poissonian repetition. No bursts were detected in the five sessions between May 27 and our GMRT detections described below; this absence has a Poissonian probability of between 0.11 and 0.37, assuming constant observing and instrument conditions. Beyond purely Poissonian variations, it is noteworthy that the rate may also intrinsically vary on timescales shorter than a month, as indicated by the apparent inconsistency between the rate derived from CHIME/FRB observations and that from a 3-hour April session at 550–750 MHz at GMRT (Marthi et al. 2022). For the CHIME/FRB fluence limit, we take the lowest-fluence burst detection from Lanman et al. (2022), namely 3 Jy ms.

Trudu et al. (2022) observed FRB 20201124A at similar frequencies, 400–416 MHz, with the Northern Cross (NC) radio telescope, for 68 hours in April and June 2021. The authors estimate a minimum detectable fluence of 44 Jy ms and expect 1 ± 1 bursts to be detectable over this whole observing campaign. That expectation is based on the rates and the power-law fluence distribution slope from CHIME/FRB monitoring (Lanman et al. 2022), and assuming there is no burst rate variability. Similarly, Kirsten et al. (2024) observed FRB 20201124A at 350 MHz with Wb-RT1 during several days, right before the L-band quenching, with an estimated minimum detectable fluence of 42 Jy ms, and not detecting any bursts.

Our observations with GMRT yielded six bursts on 2021 June 01, all of them with fluences above 40 Jy ms and one burst reaching fluence of 680 Jy ms. FAST observed FRB 20201124A a few hours before these GMRT observations and placed a stringent 0.05 Jy ms upper limit on the FRB fluences in L-band.

Taken together, these detections and strict upper limits show that after it stopped emitting at L-band, FRB 20201124A continued to produce bursts at 400–600 MHz. The evidence is shown in Fig. 15. We estimate that the low-frequency radio emission may have lasted for 3–6 days after bursts stopped in L-band. The lower limit of 3 days comes from a scenario where the L-band quenching happened right before FAST observations on May 29, and lower frequency emission ended right after the GMRT observations on June 01. The upper limit of 6 days follows from assuming the high frequencies cease right after the last detected Stockert burst on May 28 and low-frequency emission persist up to the non-detection session at the Northern Cross telescope on June 03.

Fig. 16. Maps of footpoint loops spanning all geometries. Shown are the coordinates of the footpoints of the magnetic field lines which are potentially visible for a LOS characterized by viewing angle θob, as it passes through a dipole field inclined to spin axis by an angle α. The line colors reflect the observing frequencies shown in Fig. 15. The blue loop corresponds to emission seen at WSRT frequencies νhi = 1300 MHz, coming from chosen emission altitude rhi = 5RNS; while the red loop contains the footpoints visible in the lower GMRT band νlow = 430 MHz, from rlow = 15RNS. The dots mark the quarters of the magnetar rotation. The coordinate system is centered on the magnetic dipole axis. The letters indicate the sets of geometry angles with emission regions plotted in Fig. 17.

Fig. 17. Field lines that are potentially visible, for a number of neutron-star geometries. The geometry determines the LOS, which is characterized by the viewing angle θob, and passes through a dipole field that is inclined to the spin axis by an angle α. These spin and magnetic axes are shown with thick black lines, and the NS coordinate lines are drawn in its rest frame. The frequency-color mapping follows the previous figures. Field lines with radio emission observable at 430 MHz are plotted in red, those observable at 1300 MHz in blue. Dots mark the emission points, at an altitude of 15RNS for 430 MHz and at 5RNS for 1300 MHz. The corresponding footpoint loops on the stellar surface are shown in Fig. 16. Left column: Side view with spin axis pointing vertically. Right column: Top view with spin axis pointing at the observer. Upper row: Almost aligned rotator, with our LOS passing close to the spin–magnetic axes. The emission points and field line footpoints form two almost concentric circles. Middle row: Inclined rotator with a LOS passing away from the magnetic axis, but close to the spin axis. The LOS samples a restricted range of longitudes and latitudes on the stellar surface, and the emission points and footpoints loops are completely separate for the two frequencies, meaning that the lower-frequency loop is not encircled by the higher-frequency one. Lower row: Inclined rotator with LOS passing away from the spin axis, but close to magnetic axis. The emission points and footpoints form concentric loops that are flattened on the side closer to the magnetic axis.

So far, only FRB 20180916B is known to display a similar frequency-dependent activity window in a repeating FRB (Pastor-Marazuela et al. 2021; Pleunis et al. 2021b). This source has a well-determined activity period of 16.3 days (CHIME/FRB 2020). From comparing the activity phase-resolved burst detection rate in simultaneous observations using WSRT/Apertif and LOFAR, and in earlier CHIME/FRB observations, the authors conclude that higher frequencies appear to arrive earlier in phase. This trend was confirmed to extend to frequencies as high as 5 GHz. Between 150 MHz and 5 GHz the center of activity window shifts by about 6 days and the width of the window shrinks from 3.6 to 1.0 days (Bethapudi et al. 2023).

For FRB 20201124A, the activity window is challenging to determine precisely because of the uneven observational coverage, but evidence suggests that its duration is highly variable, with the duration of the inactive sessions ranging from a few months to at least two years (Lanman et al. 2022), and the active stages spanning months. Targeted searches have ruled out a periodicity of up to 10 days (Xu et al. 2022; Niu et al. 2022; Du et al. 2024), but longer periods are not yet disproved.

5. Geometric constraints on emission regions

In the low-twist magnetar model of FRB generation, bursts are produced in the magnetosphere of a neutron star via a pulsar-like emission mechanism (WT19). The exact nature of that invoked mechanism remains a long-standing mystery despite a plethora of observational pulsar facts and decades of ongoing modeling efforts. Nonetheless, the key observational properties of pulsar radio emission can be explained by a pair of phenomenological models, known as the radius-to-frequency mapping (RFM; Cordes 1978) and the rotating vector model (RVM; Radhakrishnan & Cooke 1969). These models produce constraints on the overall magnetospheric geometry and the location of emission regions, based on a set of plausible assumptions. Below we apply RFM/RVM techniques to bursts from FRB 20201124A.

5.1. Radius-to-frequency mapping (RFM)

The RFM model postulates that radio waves decouple from magnetic field lines at a certain altitude, and propagate tangentially to the local magnetic field line at the decoupling point. We call this the emission point, although radio waves may actually be produced elsewhere (Philippov et al. 2020). Lower radio frequencies are related to emission points situated higher in the magnetosphere, where the dipolar field has diverged further, offering a natural explanation for the broadening of the on-pulse window at lower frequencies that is observed in radio pulsars. In pulsars, the on-pulse activity window is thus frequency dependent, and wider at lower frequencies.

Now, the presence of similar frequency-dependent activity window behavior is evident in both FRB 20201124A and FRB 20180916B. Nevertheless, a clear distinction arises in terms of timescale: in pulsars, the widening amounts to fractions of a second, whereas for FRB sources, the frequency-dependent edge of the activity window extends over the course of days. However, the difference in spin phase remains consistent between radio pulsars and FRB repeaters if the latter are ULPs with spin periods (P) of weeks or longer.

Qualitatively, the frequency-dependent edge of the activity window can be explained as follows: while the rotation of the magnetar slowly moves our LOS through the magnetosphere, we detect radio waves from field lines with footpoints in the active regions on the stellar surface. Plasma propagates along the active field line and generates radio emission of continuously decreasing frequency as it moves outward. Since this emission is highly beamed and its frequency is altitude-dependent, at any given moment of time the observed broadband spectrum is composed of radio waves coming from a range of altitudes and originating from different field lines. Consider, for clarity, only two radio frequencies: high and low. We, the observer, detect emission along our single LOS. But the two different frequencies originate from field lines with footpoints that form two separate, different paths on the surface of the star. These two paths may cross the edge of the active region at different times. Fortuitously, the low-frequency path may stay longer in the active region, resulting in detection of low-frequency FRBs after the cessation of higher-frequency emission. Since the effect is purely geometrical and there are no constraints on the shape of the active region for FRB 20201124A, many exact lags between high- and low-frequency quenching are possible, including negative lags, where low-frequency ceases first. The lags may also vary from one episode of activity to another. Once accumulated, the statistical distribution of lag magnitudes and signs may provide clues for the magnetospheric geometry configuration and the spread of the active regions on the stellar surface.

That is the qualitative description to introduce the underlying concepts. In what follows we constrain the location of observable field line footpoints quantitatively. We assume the external magnetic field to be dipolar, designating the angle between spin and magnetic axis as α. Since no external constraints are known for α we take it to belong to a grid of trial values ranging from 5° to 90° and spaced by 5°.

In a spherical coordinate system aligned with the spin axis, the observer LOS is defined by a pair of longitude/latitude angles. The latitude θ_ob (also called viewing angle) does not change with time while the longitude corresponds to the spin phase at the moment of time t: ϕ_ob ≡ 2πt/P. In the absence of any external constraints, we review a grid of θ_ob ranging from 5° to 90° with the spacing of 5°.

Mathematical expressions for the dipolar magnetic field have the simplest form in the rest frame of the pulsar with magnetic moment directed along z. In this coordinate system the location of the emission region is (r_em, θ_em, ϕ_em). The relationship between (θ_ob, ϕ_ob) and (θ_em, ϕ_em) is set by the requirement of the tangent to the magnetic field line at the emission point to be aligned with the LOS at the spin phase ϕ_ob. Following Lyutikov (2016):

$$\begin{array}{c}\hfill tan{\varphi }_{\mathrm{em}}=\frac{sin{\theta }_{\mathrm{ob}}sin{\varphi }_{\mathrm{ob}}}{cos\alpha sin{\theta }_{\mathrm{ob}}cos{\varphi }_{\mathrm{ob}}-sin\alpha cos{\theta }_{\mathrm{ob}}},\end{array}$$(3)

and

$$\begin{array}{c}\hfill \frac{3cos2{\theta }_{\mathrm{em}}+1}{\sqrt{6cos2{\theta }_{\mathrm{em}}+10}}=cos\alpha cos{\theta }_{\mathrm{ob}}+sin\alpha sin{\theta }_{\mathrm{ob}}cos{\varphi }_{\mathrm{ob}}.\end{array}$$(4)

The coordinates of the footpoint of an active field line (θ₀, ϕ₀) can be obtained using the equation for that dipolar field line in the rest frame:

$$\begin{array}{c}\hfill \frac{R_{\mathrm{NS}}}{{sin}^2{\theta }_0}=\frac{r_{\mathrm{em}}}{{sin}^2{\theta }_{\mathrm{em}}},\end{array}$$(5)

and

$$\begin{array}{c}\hfill {\varphi }_0={\varphi }_{\mathrm{em}}.\end{array}$$(6)

Without constraints on the emission altitude r_em, multiple field lines can contribute to the emission at any given ϕ_ob. To maximize the stellar surface area under consideration in our exploration, we set the altitude of the observed higher-frequency radio emission to the smallest possible value. For the frequency at the center of WSRT band, this value corresponding to the minimum altitude from which radio waves can escape the low-twist magnetosphere. Assuming characteristic values of the limiting twist and crust oscillation frequency, and taking surface magnetic field to be 10¹⁴ G, we find this minimum escape altitude is 5R_NS (WT19, Beniamini et al. 2020).

While our numerical calculations are grounded on the assumption that the observed radio waves originate from curvature radiation (Wang et al. 2019, and references therein), the outcome is qualitatively similar when employing a different relationship between the radio wave frequency and the altitude of the emission region. For curvature radiation, the relationship between the radio frequency and the plasma/magnetic field parameters is

$$\begin{array}{c}\hfill \nu =\frac{3\pi {\gamma }^{3}c}{4\kappa },\end{array}$$(7)

where γ is the Lorentz factor of the emitting particles, c is the speed of light, and κ is the curvature radius of the field lines. For a dipole magnetic field, κ can be calculated from the coordinates of the emission point (r_em, θ_em, ϕ_em) in the rest frame of the star:

$$\begin{array}{c}\hfill \kappa =\frac{r_{\mathrm{em}}{(1+3{cos}^2{\theta }_{\mathrm{em}})}^{3/2}}{3sin{\theta }_{\mathrm{em}}(1+{cos}^2{\theta }_{\mathrm{em}})}·\end{array}$$(8)

Combining Eqs. (3)–(8), we obtain (r_em, θ_em, ϕ_em) for both the WSRT and GMRT radio frequencies ν_hi and ν_lo, for a range of chosen inclination and viewing angles. If γ does not change between two frequencies (Wang et al. 2019), r(ν_hi)≈0.3r(ν_lo) over the entire range of input parameters. This next determines the low-frequency emission to take place at 16.7R_NS, above the minimum escape altitude of 15R_NS for GMRT frequencies. If γ_lo < γ_hi, then r_lo is closer to r_hi, and the lower-frequency θ₀(t) moves closer to its higher frequency counterpart.

The absolute value of γ (if taken constant between two frequencies) varies from ∼50 to ∼200, with the higher values corresponding to field lines closer to the magnetic pole. These values are similar to the γ = 300 adopted by Wang et al. (2019). We note that our values of γ are lower limits, since we adopted the smallest possible r_em.

In the absence of any other constraints, the footpoints of the field lines with emission visible to the observer form a closed loop on the star surface (Fig. 16). For the chosen frequencies ν_hi and ν_lo, loops subtend ≲0.2% of star surface for small θ_ob, regardless of inclination angle. The fraction of the stellar surface encircled by this footpoint loop increases with growing θ_ob, reaching about 20% for large viewing and small inclination angles.

Since the footpoint loops are closed, there exists a field line that can potentially provide radio emission at any spin phase ϕ_ob. Whether this line will actually emit is determined by physical requirements for pair production and radio wave generation. For radio pulsars, such a requirement is that radio emission is produced by the open field lines with θ₀ ≤ θ_PC, where θ_PC is the radius of the polar cap. Figure 18 provides an illustration of this, showing an example of the footpoint loops at two radio frequencies for an inclination angle of α = 30° and θ_ob = 60°. For plotting convenience, r_lo was set to be 2R_NS and ν_hi = 1.5ν_lo. The polar cap with θ_PC = 15° cuts segments from the footpoint loops, with the lower-frequency segment spanning a larger amount of spin phase bins. This implies that radio emission at lower frequencies is observed at earlier and later ϕ_ob than the high-frequency emission; the average pulse profile is wider at lower frequencies.

Fig. 18. Radius-to-frequency mapping of radio emission observed at two radio frequencies together with the footpoints of visible emission lines on the surface of the star for radio pulsars (left column) and FRB 20201124A (right column). Spin and magnetic axis are marked with vertical and inclined lines, respectively. The path of the LOS and corresponding line footpoints are shown with circle markers, gray for no emission, and red (resp. blue) for low (resp. high) frequency. For radio pulsars, only open field lines originating within the polar cap (green circle) can produce radio emission. At lower frequencies the LOS spends more time within the polar cap, resulting in a widening of observed radio profile. For FRB 20201124A, emission is not restricted to the open field lines and comes from an active region of unknown shape. At the end of the S21 activity window the LOS crossed the edge of an active region (right edge on the figure) in two places, with the edge point corresponding to the lower radio frequency extending farther along the LOS path.

In the WT19 model the FRB 20201124A emission can originate on closed field lines, eliminating the polar cap requirement. The minimum emission altitude requirements imposed by the plasma transparency set limits on the magnetic colatitude θ₀, defined by r_min, alt = R_NSsin²θ₀, indicating that θ₀ ≲ 15° and 26° for lower and higher frequencies, respectively. This broad region may host several active patches during the activity window.

A frequency-dependent activity edge arises when the LOS crosses the edge of the active patch. The exact lag between high- and low-frequency quenching is determined by the skewness of the active region shape in the magnetic longitude direction and the unknown value of the spin period. For instance, assuming the FRB 20201124A spin period to be 2 months (the duration of the S21 activity window), then δϕ_ob = 20°–40° for a 3–6 day lag, which translates to a similar upper limit on δϕ_em, with the latter being smaller for non-concentric footpoints loops. Thus, a small variation in active region shape may lead to a large variation in observed cessation lag.

A frequency-dependent activity window places constraints on the overall size of an active surface area. Individual crustal motions, on the other hand, are limited in size by the trombone frequency drift (Wang et al. 2019; Bilous et al. 2022; Lyutikov 2020). Such frequency drift by δν = 300 MHz at 1400 MHz or by 125 MHz at 650 MHz would translate to δθ₀ ≈ 0.5θ₀δν/ν ≈ 0.1θ₀ ≲ 3°. If the size of the active “spark” which produces an individual FRB is comparable to the amplitude of footpoint dislocation ξ, then the latter is less than 500 m, below ξ_max = 3 km from Beniamini et al. (2020).

5.2. The rotating vector model (RVM)

In the rotating vector model, the plane of the linearly polarized part of the radio emission is defined by the field line curvature plane at the emission point. For an inclined dipolar field, the position angle (PA) has a characteristic S-shape dependence on the spin phase ϕ_ob and the angles introduced previously, as defined by the following analytical relation:

$$\begin{array}{c}\hfill tan\mathrm{PA}=\frac{sin\alpha sin{\varphi }_{\mathrm{ob}}}{sin{\theta }_{\mathrm{ob}}cos\alpha -cos{\theta }_{\mathrm{ob}}sin\alpha cos{\varphi }_{\mathrm{ob}}}·\end{array}$$(9)

For both normal radio pulsars and magnetars, individual pulses exhibit a large diversity of PA behavior (Mitra et al. 2016; Johnston et al. 2024; Kramer et al. 2007). While numerous deviations are known (Johnston et al. 2023), the measured single-pulse PAs average to an S-curve defined by Eq. (9) in a considerable number of cases.

If we refocus now on FRB 20201124A, we see that the behavior of the PA curve is known only for a four-day stretch at the end of the F21 activity window (Jiang et al. 2022). There, the collection of the mean PAs of individual bursts had a relatively wide distribution around some most probable value, which in turn remained constant during four observing days. The authors do not specify the limits of the constancy, but the distribution of mean PAs in Jiang et al. (2022) has a standard deviation of about 20°. For the purpose of the arguments following below, we adopt a PA-gradient upper limit of 10° per four days as best describing the observations.

We proceed by selecting (α, θ_ob) pairs for which there exists at least one spin longitude window that can produce the observed |max(PA)−min(PA) |< 10°. The size of the spin window depends on the unknown spin period, so we explore three trial periods, of one, two, or four months. Figure 19 shows the allowed ranges of inclination/viewing angles. For the smallest trial period of one month α ≲ θ_ob unless the inclination angle is close to 90°. This excludes the non-concentric types of footpoint loops that are visible in Fig. 16. The inclination angle becomes less constrained as the potential spin period grows. Combinations of large α and small θ_ob (e.g., α = 80° and θ_ob = 10°) still, however, produce a PA gradient that is too steep, even if P = 4 months.

Fig. 19. Top: Constraints on the magnetospheric geometry placed by constant median single-burst PA during the four-day observing epoch from the F21 activity window (Jiang et al. 2022). The shaded areas encompass inclination–viewing angle pairs related to three trial values of spin period P (see legend). For each of these areas, the maximum PA variation was less than 10° in at least one of the spin phase regions spanning 4 days/P. The stars mark three representative geometries whose PA curves are shown in the bottom plot. Geometries that approach aligned rotators and have small viewing angles predominantly show larger PA gradients than allowed by the observations, unless the spin period exceeds several months.

6. Summary and conclusions

FRB 20201124A is a prolific source of FRBs. Its bursts exhibit rich and complex structure and are undoubtedly capable of providing valuable insights into the underlying emission physics and propagation mechanisms. However, the current body of observational evidence contains gaps in areas of high scientific interest, and is inevitably affected by both instrumental and processing limitations.

6.1. Determining DM and fluence changes requires careful analysis

One example of these limitations is our current inability to disentangle the dispersive delay imposed by the interstellar medium from the intrinsic burst properties. We have demonstrated that the DM measured using the de facto standard structure-maximizing method absorbs the trombone drift for fainter pulses, leading to an overestimation bias in DM measurements that can be as large as 10 pc cm⁻³, even for bursts with an integrated S/N on the order of 100. For brighter bursts, the complex structure often cannot be maximized with a single DM value (e.g., Zhou et al. 2022). This bias should be taken into account while analyzing DM distributions aggregated from several studies with different sensitivities since it can lead to apparently multi-modal DM distributions, which can subsequently be interpreted as plasma lensing or frequency-dependent DM (Wang et al. 2023b). On the other hand, directly incorporating external constraints on spectro-temporal FRB behavior would allow for testing interesting physical effects. One example is investigating the filamentation of FRBs in the relativistic winds of magnetars, which implies a dependence of the measured DM on pulse luminosity due to propagation in the near-source environment (Sobacchi et al. 2023).

We confirm the absence of secular DM trends between different activity windows of FRB 20201124A down to the level of a few pc cm⁻³. This places a limit on any electron density changes along the LOS that could have been caused by an active circum-burst environment (Metzger et al. 2017). Combining WSRT bursts from the W22 activity window and burst G05 from the S21 activity window, we obtain a broadband PL dependence of the decorrelation bandwidth, shallower than expected for Kolmogorov turbulence and close to the measurements of Main et al. (2022), who had a similar observing setup.

Despite the relative abundance of available data, the FRB 20201124A burst-fluence distribution remains poorly constrained. The FAST observations demonstrate that there is a significant amount of variability in the fluence distribution over a range of timescales: between subsequent days to between distinct activity windows (Xu et al. 2022; Zhang et al. 2022). Corroborating the findings of Kirsten et al. (2024), we observed a shallower PLI for bursts in the W22 activity window than in the FAST measurements during previous windows (all at frequencies around 1300 MHz). However, as we demonstrate, PLI measurements tend to be imprecise and biased for small burst samples, especially when drawn from distributions with flattening at the low-fluence end, even if the observed distribution appears to be well fit with a single PL.

6.2. A ULP magnetar model for the burst behavior

We explored the possibility that FRB 20201124A is an ultra-long-period magnetar with a spin period on the scale of months, a period that may remain undetectable for current searches. For such ULP sources, FRBs could be triggered by motion in the star crust, if these cause dislocation of magnetic field-line footpoints that leads to the production of plasma, which subsequently emits radio waves via some pulsar-like emission mechanism (WT19). Similar to the short X-ray bursts thought to be triggered by similar crustal motion events, the arrival times of individual FRBs would then follow a log-uniform distribution if measured from the trigger moment. Such a distribution has been reported for the repeating source FRB 20121102A (WT19). We demonstrated that TOAs of individual bursts from FRB 20201124A in both the S21 and F21 FAST samples exhibit a distribution that is uniform, not log-uniform. The sample could, however, still actually be part of the tail of a log-normal distribution, from a trigger event that occurred before these FAST observing sessions. As the burst rate, in a crescendo, increases exponentially before the emission quenches at the end of the F21 activity window (Zhang et al. 2022), there must be multiple trigger events during a single activity window. If these occurred several hours before the observation, our data is consistent with this model.

Another prediction from the crust motion and low-twist theory links crustal oscillations to quasi-periodic sub-burst spacing within individual FRBs. Our WSRT sample contains three bursts with seemingly periodic components; however, we were unable to find any significant periodicity. We show (and caution) that detecting periodicity for closely spaced sub-bursts of arbitrary shape is difficult, especially when the number of period trials is unconstrained.

If radio emission originates from only one frequency-dependent range of heights, then at any given spin phase, the LOS samples only a small fraction of the stellar surface. In contrast, X-ray bursts can be detected from the entire stellar surface at any moment. This difference may explain the intriguing result obtained by Tsuzuki et al. (2024), who found that regular radio pulses from the occasionally FRB-emitting magnetar SGR 1935+2154 exhibit correlation properties similar to those of extragalactic FRBs (Totani & Tsuzuki 2023) and earthquakes. However, X-ray bursts from the same source did not show such correlation properties. We hypothesize that the lack of correlation for X-ray bursts stems from the superposition of several starquake events occurring simultaneously at different locations. In radio, individual starquake events are observed sequentially.

6.3. The FRB 20201124A activity window ends chromatically

At the end of the S21 activity window, our GMRT observations recorded several strong bursts at 300–600 MHz. This is surprising because only a few hours before, FAST had placed upper limits on the absence of emission at 1250 MHz that were 1000 times more stringent. By combining burst detections and telescope scheduling information, we were able to paint a complete multi-frequency picture of the end of this activity window. We showed that the edge of the S21 activity window depends on the observing frequency: low-frequency emission is present for 3–6 days after the higher-frequency quenching. This marks the second detection of a chromatic activity window for a repeating FRB source, following FRB 20180916B (Pastor-Marazuela et al. 2021; Pleunis et al. 2021b). However, for FRB 20180916B, a 16.3-day periodicity is known, and the chromatic window was established by combining bursts from several periods, whereas for FRB 20201124A, we have so far recorded only a one-time event. Unfortunately, the end of the F21 and W22 activity windows was not covered by multi-frequency observations.

The fact that lower-frequency emission was recorded a few days after the disappearance of the higher-frequency emission would be hard to explain within the plasma lensing theory of burst quenching. In this theory, the non-uniform distribution of free electrons along the LOS acts as a diverging lens, causing a significant increase in burst rate followed by a sudden drop at the end of the F21 activity window (Chen et al. 2024). For the one-dimensional Gaussian lens model reviewed by Chen et al., low-frequency emission disappears before the high-frequency emission. This is the opposite of what we observed.

6.4. Two well-known pulsar models can explain the chromatic activity end, and the flat polarization, assuming a rotation period of order one month

The chromatic activity window we observed resembles a long-known effect observed for radio pulsar emission, namely the widening of the on-pulse window toward lower frequencies. This similarity can actually have a physical interpretation if FRBs are indeed generated via a pulsar-like emission mechanism. We applied the classical phenomenological radius-to-frequency model of average profile widening to the frequency-dependent edge of the S21 activity window. Our goal was to constrain the magnetospheric geometry, and the location of active regions on the star surface. Assuming that the active regions are not confined to the open field lines, and given the fact that the spin period is unconstrained, we find a multitude of possible active regions for any combination of dipole inclination angle and the viewing angle of the LOS. The possible extent of these regions in magnetic colatitude is only limited by the requirement of plasma transparency and is on the order of ∼20°. The observed chromatic edge of the activity window is defined by the exact shape of the active regions, and small variations in this shape may cause different amounts of chromatic lag. If active regions differ from one activity window to another, we expect the amount of frequency-dependent quenching timescale to vary in subsequent activity windows as well. In this model an inverted dependency (i.e., low-frequency emission disappearing before the high frequencies) is also allowed. Similarly, a chromaticity at the onset of the activity might also exist, again depending on the exact shape of the active regions. In the absence of a physical constraint on the active region shapes, all possibilities, such as symmetric, asymmetric, or one-sided behavior of the activity window would be possible. It is worth noting that FRB 20180916B exhibits such asymmetric chromaticity, with the emission window at 5 GHz preceding the one at 150 MHz (Bethapudi et al. 2023).

The magnetospheric geometry of FRB 20201124A can be further constrained using a second classical aspect of our understanding of pulsars, the rotating vector model. For this, we determined which geometries are allowed under the observed PA constancy during the four-day observing stretch at the end of the F21 activity window (Jiang et al. 2022). For potential spin periods of a month and above, we can exclude aligned rotators with a LOS close to the spin axis. The constraints are less strict for larger potential spin periods.

6.5. Forward look

To the best of our knowledge, as of July 2024 there have been no indications that FRB 20201124A entered a new activity window, with the last pulses detected in the end of March 2022 (Wu et al. 2024). Continuous monitoring of this source remains critically important by, among others, CHIME on the low-frequency end and a collection of small telescopes operating at higher frequencies (Ould-Boukattine et al. 2024).

We showed that catching a new activity window early is essential for testing the crust motion and low-twist FRB emission theory, as future observations may unveil a prolific cluster of bursts that could mark the onset of a single trigger event. Additionally, such observations could provide valuable insights on the potentially diverse shape of the active regions through the recording of another chromatic edge of the activity window, after the first one presented here. Finally, further constraints on the secular behavior of the PA would be helpful in refining our understanding of magnetospheric geometry.

Data availability

To facilitate the future use in the community of the correct, MLE algorithms for fitting burst fluences, we have made available an ipython notebook that implements the various types of PL fits. It includes instructions on how to adapt the powerlaw package as described in Appendix A. The notebook is hosted on Zenodo (https://doi.org/10.5281/zenodo.12644702) and GitHub (https://github.com/TRASAL/FRB_powerlaw).

Acknowledgments

AB thanks Matteo Trudu, Adam Lanman, Jian-Ping Yuan, and Franz Kirsten for providing extra information about FRB 20201124A observations, Kejia Lee for help with downloading FAST data, Zorawar Wadiasingh and Andrey Timokhin for discussions on the low-twist ULP magnetar model, and Heng Xu for providing critical feedback on fluence distributions. This work was supported by the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement No. 617199 (‘ALERT’) and Vici research programme ‘ARGO’ with project number 639.043.815, financed by the Dutch Research Council (NWO). JVL further acknowledges support from ‘CORTEX’ (NWA.1160.18.316), under the research programme NWA-ORC by NWO. YM further acknowledges support from the Department of Science and Technology, India, via the Science and Engineering Research Board (SERB) Start-up Research Grant (SRG/2023/002657). IPM further acknowledges funding from an NWO Rubicon Fellowship, project number 019.221EN.019. We would like to thank the Centre Director, NCRA-TIFR, and the GMRT observatory for the prompt director’s discretionary time allocation and scheduling of our observations. We thank the staff of the GMRT who have made the GMRT observations possible. The GMRT is run by the National Centre for Radio Astrophysics of the Tata Institute of Fundamental Research. This work makes use of data from the Apertif system installed at the Westerbork Synthesis Radio Telescope owned by ASTRON. ASTRON, the Netherlands Institute for Radio Astronomy, is an institute of NWO.

References

Appendix A: Comparing least-squares and MLE methods of PLI estimates

In this section we investigate accuracy and bias of PLI estimates provided by the graphical and MLE methods. We start with demonstrating the sub-optimal performance of the least-squares estimates. To do this, we generated three realizations of PL-distributed random variables using a Pareto distribution from the python library scipy, with survival function defined as:

$$\begin{array}{c}\hfill N(X\ge X_0)=N{\left(\frac{X_0}{X_{\min }}\right)}^{\alpha }.\end{array}$$(A.1)

We took N = 100, X_min = 1, and α = −2. PLIs were estimated using graphical method on the unbinned survival functions. In order to investigate the dependence of the PLI produced by the graphical method on the sample size, we performed fits on reduced samples, selecting X ≥ X_min, where X_min ranged from the smallest to the fourth largest values of X. For each fit both the least-squares and bootstrapping errors were estimated. For the latter we followed the procedure from Kirsten et al. (2024), removing 10% of the sample without replacement (or one element if N(X_min) < 10), calculating the PLI, repeating this 100 times and then taking standard deviation of the acquired distribution as ϵ_α, the error on α. Figure A.1 displays the survival functions and the resulting α(N). Even for N = 100, the resulting PLI can deviate from the true value by as much as 25%, while the formal errors are an order of magnitude smaller.

Fig. A.1. Left: Survival functions for three N = 100 samples of simulated random variables distributed according to Eq. A.1 with αt = −2. Right: PLI αm obtained with the least-squares method for the subsamples consisting of 4 ≤ N ≤ 100 elements sorted in decreasing order. The shaded regions mark the standard least-squares error (darker shade) and the bootstrapping error (lighter shade).

The uncertainties in PLI measurements calculated using the bootstrapping method are indicative at best. The distribution of α_m on bootstrapped subsamples is generally skewed toward shallower PLIs. The shape of the SF tail can be heavily influenced by a few large random variables, which constitute a small fraction of the entire sample. Depending on the chance presence of these variables in the bootstrapped subsample, the measured PLI can vary substantially, resulting in multi-modal bootstrapped PLI distributions. This has a significant effect on the calculated standard deviation of bootstrapped PLIs, causing large scatter in the bootstrapped error on α_m from one simulation to another. At the same time, the error is not large enough to describe the true discrepancy between α_t and α_m (Fig. A.1).

To further test the quality of the PL fitting, we examined the distributions of test statistics a, a = (α_m − α_t)/α_err for several fitting methods. We generated a series of samples with N ranging from 5 to 500, values for α_t and −2. Performing smaller tests using different values of α_t did not reveal any differences in the behavior of examined distributions.

Fig. A.2 shows the accuracy of the fit for three methods: graphical with its standard errors, graphical with the bootstrap errors, and MLE (where we emphasize here that the three relevant ordinates in that figure are on different scales). Following James et al. (2019), we used unbiased α estimates:

Fig. A.2. Accuracy of three methods of α estimation vs sample size N. Left: Graphical method with its least-squares error; center: Graphical method with a bootstrap error; and right: Unbiased MLE method according to Eqs. A.2–A.3. We note that the y-axis scale decreases by over a factor of 10 from left to right. The pink line corresponds to the median of (αm − αt)/αerr, and the shaded regions mark percentiles corresponding to ±1, ±2, ±3 standard deviations of the normal distribution.

$$\begin{array}{c}\hfill \frac{1}{{\alpha }_{\mathrm{m}}}=\frac{1}{N-1}{\displaystyle \sum \ln }\left(\frac{X}{X_{\min }}\right).\end{array}$$(A.2)

For α_err we used the expression from the same work, despite it being valid only for large N:

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{err}}=\frac{{\alpha }_{\mathrm{m}}}{\sqrt{N-2}}.\end{array}$$(A.3)

In the ideal case a follows a normal distribution with a median of 0 and a standard deviation σ = 1. This means that α_m provides an unbiased estimate of α_t, and α_err has an intuitive Gaussian meaning. To illustrate how closely a obeys a normal distribution we plot, in Fig. A.2 (right sub-panel), the percentiles at levels corresponding to the median and its ±σ, ±2σ, and ±3σ counterparts. For the MLE method and relatively large sample sizes N ≳ 100, these levels are at a = 0, ±1, ±2, ±3, indicating that a very closely resembles the normal distribution. For smaller sample sizes a is skewed toward negative values, and accurately described by a gamma distribution (Eq. 2). We note that the median value of a stays close to 0 for all sample sizes.

In stark contrast to this, the graphical method leads to a median value of a on the order of 5 for the least-squares error and 1.5 for the bootstrap error, indicating that median α_m tends to be larger than α_t by a few α_err, regardless of sample size. The wide spread of the percentile curves in Fig. A.2 (left and middle sub-panel) demonstrates that α_err values are substantially underestimated when choosing to use the graphical method.

We performed similar accuracy estimates for PLIs determined with the Python package powerlaw¹⁹ (Alstott et al. 2014). In the original code, α and α_err are calculated as

$$\begin{array}{c}\hfill \frac{1}{{\alpha }_{\mathrm{m}}}=\frac{1}{N}{\displaystyle \sum \ln }\left(\frac{X}{X_{\min }}\right),\end{array}$$(A.4)

and

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{err}}=\frac{{\alpha }_{\mathrm{m}}}{\sqrt{N}}.\end{array}$$(A.5)

Fig. A.3 (left panel) shows that this estimate performs worse than the version modified according to Eqs. A.2-A.3 (center panel).

Fig. A.3. Similar to Fig. A.2, but for accuracy of α estimates using the powerlaw package with different modifications. Left: Original package version with biased estimates (Eqs. A.4–A.5) and fixed Xmin; center: Unbiased estimates with fixed Xmin; and right: Unbiased estimates with fitted Xmin. In all three cases a random variable was simulated with single-PL distribution according to Eq. A.1.

In real-life observations X_min is often not known. The powerlaw package offers fitting for the best X_min by minimizing the Kolmogorov-Smirnov distance between the data and the theoretical power-law fit. We have repeated the aforementioned simulations with fitting for X_min, restricting its range from min(X) to the third largest value. In this case, the median value of a is slightly more biased toward small negative values and the distribution is somewhat wider (Fig. A.3), but the MLE still provides a much better fit than the graphical method.

All Tables

All Figures

Fig. 1. Spectra and band-integrated sample intensity profiles for the 16 bursts from GMRT observations. Bursts G01–G06 were detected on 2021 June 01, the rest on 2022 February 06–07. All bursts are dedispersed to DM = 414.73 pc cm−3 determined from burst G05. For plotting, the spectra were normalized by the mean and standard deviation in the off-burst region in each sub-band and the colors were saturated at ±8σ for G05 and at ±4σ for the other bursts. All bursts from 2021 except for G05 are plotted with a time resolution of 5.24 ms and a frequency resolution of 22–25 MHz. Burst G05 is shown on a finer 2.62 ms/391 kHz spectro-temporal grid. For the 2022 bursts, the time resolution is 3.3–5.9 ms, and the frequency resolution is 24 MHz.

In the text

Fig. 2. Spectra (bottom subpanels) and intensity profiles, integrated over the frequency band (top subpanels), of the ten bursts discovered within the ALERT survey. For plotting, the spectra were normalized by the mean and the standard deviation in the off-burst region in each sub-band, and the colors were saturated at ±10σ for W01, W07, and W10, and at ±4σ for the other bursts. Bursts W01 and W10 are plotted with the original tres = 81.9 μs, burst W07 with 0.33 ms and the rest with 2.62 ms. The frequency resolution is 25 MHz for all bursts except W07 and W10 (6.25 MHz), and W01 (0.48 MHz). The time in ms is zeroed on the burst peak. Spectra at frequencies above 1450 MHz are not shown, since the signal there is mostly corrupted by RFI.

In the text

	Fig. 3. Compilation of DM measurements from previous studies: Xu et al. (2022, X21), Kumar et al. (2022, K22), Lanman et al. (2022, L22), and this work, plotted as a function of integrated S/N.
In the text

Fig. 4. Influence of noise on the DM estimate from DM_phase. Burst #377 from Fig. 6 in the extended data in Xu et al. (2022) was normalized in the off-burst window in each channel (such that the standard deviation of the noise std = 1), and Gaussian noise was added in each channel. Each simulation was then run through DM_phase with an automatic spectrum cut. Top row: Example of spectra before and after DM correction for three different values of noise std. Middle row: Corresponding DM_phase spectra. Bottom row:ΔDM vs. a set of noise realizations, six simulation runs per each noise std. For clarity, some jitter along the horizontal axis was added to each set of points.

In the text

	Fig. 5. Burst G05 dedispersed with structure-maximizing DM = 414.73 pc cm−3 (left) and with DM = 411 pc cm−3 (right), suggested by the brightest bursts observed by FAST.
In the text

	Fig. 6. Left column: Bursts W01, W07a, W07b, and W10 with potential quasiperiodic components. Right: Autocorrelation function from all four bursts.
In the text

	Fig. 7. Top: TOAs of individual sub-bursts of FRB W07 (red dots) identified as peaks of the signal smoothed with Gaussian kernel (thick line). The vertical lines show periodicity with best quasi-period P = 1.44 ms. Middle: Phase residuals with respect to model based on the best-fit quasi-period P. Bottom:Z2 statistic on the range of trial periods.
In the text

Fig. 8. Upper row: FRB fluence vs. MJD of arrival for three activity windows. The markers and colors correspond to the legend in the lower panel (see inset), which shows rate survival function vs. burst fluence. FRB information was collected from Lanman et al. (2022, Lan22), Xu et al. (2022, Xu21), Nimmo et al. (2022, Nim22), Hilmarsson et al. (2021, Hil21), Kumar et al. (2022, Kum22), Marthi et al. (2022, Mar22), Zhang et al. (2022, Zha22), and Kirsten et al. (2024, Kir24). Where possible, fluence distributions obtained from observations on different days and with different telescopes were computed separately.

In the text

Fig. 9. Left: Survival functions of the fluence distributions from the FAST and Onsala observations during the S21 activity window, together with MLE PL fits. The vertical lines mark the optimal minimum fluence as found by the powerlaw package, after modification according to Eqs. (A.2)–(A.3) in Appendix A. Top right: Posterior distributions of αt|αm. The gray lines indicate gamma distributions from Eq. (2). Bottom right: Probability of αt(Onsala)−αt(FAST) being smaller than the threshold value.

In the text

Fig. 10. Survival functions of the fluence distribution from the S21 FAST bursts (Xu et al. 2022) with fluences larger than 0.5 Jy ms, together with an example survival function from a random ten-pulse subsample. The vertical lines indicate the optimal minimum fluence values from the PL fit. The example survival function does not show apparent flattening. Inset: Distribution of αm for all 104 ten-pulse subsamples (green) and those that do not show flattening (orange).

In the text

Fig. 11. Un-normalized survival function of the fluence distributions (red line) for the bursts recorded at the end of F21 activity window (Zhang et al. 2022). The gray rectangles mark the regions in the fluence–burst number parameter space that contain bursts excluded from the PL fit (see Section 4.3.2 for details). The colored rectangles contain bursts included in the fit; the color represents the steepness of the best-fit power-law index (PLI) ranging from light yellow (αm = −0.5) to dark blue (αm = −2.5). Light blue bands of progressively lighter shades in the background correspond to regions containing 10, 100, and 1000 bursts.

In the text

	Fig. 12. Un-normalized survival function of the fluence distributions for the bursts recorded in the previous studies and current work (see caption of Fig. 11 for the explanation of the legend and caption of Fig. 8 for the literature abbreviations). Here the burst samples are ordered by the number of bursts in the best-fitting PL tail. MJD refers to the start of respective observing campaigns.
In the text

Fig. 13. Decorrelation bandwidth analysis. The top panel shows the ACF from 1220 MHz to 1453 MHz of the Apertif bursts W01 (orange), W07a (green), W07b (red), W10 (purple), and G05 (blue). The middle panel shows the ACF from 580 to 709 MHz of the uGMRT burst G05 in blue. The bottom panel shows how the decorrelation bandwidth evolves with frequency. The Apertif bandwidths were divided into five sub-bands, and the uGMRT bandwidth into four sub-bands. The gray dashed line with a shaded area shows the best-fitting power law, with an index α = 2.4 ± 0.2. As a reference, the black solid line shows an index α = 4 expected from a thin screen. The gray shaded region indicates decorrelation bandwidths that cannot be resolved by the uGMRT resolution.

In the text

Fig. 14. TOA and waiting time distribution for the most prolific FRB 20201124A observing session from Zhang et al. (2022) on MJD 59486 (upper row), and a series of log-uniform random variables (lower row) that simulate an FRB outburst. In our model only a small fraction of the outburst duration is covered by observations, long after the burst onset. The columns are: a) distribution of TOAs measured from the first burst arrival (brown) and from the hypothetical outburst start (violet), with the pale shade of violet showing the whole underlying mostly unobserved outburst; b) same distribution, but binned with logarithmic time bins; c) distribution of TOAs vs. waiting time between two bursts; d) distribution of waiting times. Bursts separated by less than 400 ms combined, producing a single TOA for the start of the combined burst.

In the text

Fig. 15. Upper panel: Overview of daily detection statistics at the end of S21 activity window from all telescopes combined. A filled rectangle indicates that the source was detected, a hatched rectangle for non-detected, and a white rectangle signifies there were no observations that day. Lower frequencies are shown at the top, following Sect. 5.1. Middle panel: FRB fluences, together with upper limits for non-detections vs. the integer part of the observation MJD. The time region around the end of the S21 activity window is shown. The red triangles mark GMRT observations at 400 or 650 MHz (bursts G01–G06 from Fig. 1). The darker and lighter red shades show detections (pentagons) or upper limits for bursts in similar frequency regions from the observations by CHIME/FRB (Lanman et al. 2022) and NC (Trudu et al. 2022). The results of higher-frequency L-band observations by FAST (Xu et al. 2022) and NSRT (Mao et al. 2022) are shown in darker and lighter blue, with dots marking FAST detections. Detections by Torun/Onsala/Stockert (T/O/S) dishes from Kirsten et al. (2024) are shown as blue triangles (see text for details on fluence upper limits). Lower panel: Observation time coverage.

In the text

Fig. 16. Maps of footpoint loops spanning all geometries. Shown are the coordinates of the footpoints of the magnetic field lines which are potentially visible for a LOS characterized by viewing angle θob, as it passes through a dipole field inclined to spin axis by an angle α. The line colors reflect the observing frequencies shown in Fig. 15. The blue loop corresponds to emission seen at WSRT frequencies νhi = 1300 MHz, coming from chosen emission altitude rhi = 5RNS; while the red loop contains the footpoints visible in the lower GMRT band νlow = 430 MHz, from rlow = 15RNS. The dots mark the quarters of the magnetar rotation. The coordinate system is centered on the magnetic dipole axis. The letters indicate the sets of geometry angles with emission regions plotted in Fig. 17.

In the text

Fig. 17. Field lines that are potentially visible, for a number of neutron-star geometries. The geometry determines the LOS, which is characterized by the viewing angle θob, and passes through a dipole field that is inclined to the spin axis by an angle α. These spin and magnetic axes are shown with thick black lines, and the NS coordinate lines are drawn in its rest frame. The frequency-color mapping follows the previous figures. Field lines with radio emission observable at 430 MHz are plotted in red, those observable at 1300 MHz in blue. Dots mark the emission points, at an altitude of 15RNS for 430 MHz and at 5RNS for 1300 MHz. The corresponding footpoint loops on the stellar surface are shown in Fig. 16. Left column: Side view with spin axis pointing vertically. Right column: Top view with spin axis pointing at the observer. Upper row: Almost aligned rotator, with our LOS passing close to the spin–magnetic axes. The emission points and field line footpoints form two almost concentric circles. Middle row: Inclined rotator with a LOS passing away from the magnetic axis, but close to the spin axis. The LOS samples a restricted range of longitudes and latitudes on the stellar surface, and the emission points and footpoints loops are completely separate for the two frequencies, meaning that the lower-frequency loop is not encircled by the higher-frequency one. Lower row: Inclined rotator with LOS passing away from the spin axis, but close to magnetic axis. The emission points and footpoints form concentric loops that are flattened on the side closer to the magnetic axis.

In the text

Fig. 18. Radius-to-frequency mapping of radio emission observed at two radio frequencies together with the footpoints of visible emission lines on the surface of the star for radio pulsars (left column) and FRB 20201124A (right column). Spin and magnetic axis are marked with vertical and inclined lines, respectively. The path of the LOS and corresponding line footpoints are shown with circle markers, gray for no emission, and red (resp. blue) for low (resp. high) frequency. For radio pulsars, only open field lines originating within the polar cap (green circle) can produce radio emission. At lower frequencies the LOS spends more time within the polar cap, resulting in a widening of observed radio profile. For FRB 20201124A, emission is not restricted to the open field lines and comes from an active region of unknown shape. At the end of the S21 activity window the LOS crossed the edge of an active region (right edge on the figure) in two places, with the edge point corresponding to the lower radio frequency extending farther along the LOS path.

In the text

Fig. 19. Top: Constraints on the magnetospheric geometry placed by constant median single-burst PA during the four-day observing epoch from the F21 activity window (Jiang et al. 2022). The shaded areas encompass inclination–viewing angle pairs related to three trial values of spin period P (see legend). For each of these areas, the maximum PA variation was less than 10° in at least one of the spin phase regions spanning 4 days/P. The stars mark three representative geometries whose PA curves are shown in the bottom plot. Geometries that approach aligned rotators and have small viewing angles predominantly show larger PA gradients than allowed by the observations, unless the spin period exceeds several months.

In the text

	Fig. A.1. Left: Survival functions for three N = 100 samples of simulated random variables distributed according to Eq. A.1 with αt = −2. Right: PLI αm obtained with the least-squares method for the subsamples consisting of 4 ≤ N ≤ 100 elements sorted in decreasing order. The shaded regions mark the standard least-squares error (darker shade) and the bootstrapping error (lighter shade).
In the text

Fig. A.2. Accuracy of three methods of α estimation vs sample size N. Left: Graphical method with its least-squares error; center: Graphical method with a bootstrap error; and right: Unbiased MLE method according to Eqs. A.2–A.3. We note that the y-axis scale decreases by over a factor of 10 from left to right. The pink line corresponds to the median of (αm − αt)/αerr, and the shaded regions mark percentiles corresponding to ±1, ±2, ±3 standard deviations of the normal distribution.

In the text

	Fig. A.3. Similar to Fig. A.2, but for accuracy of α estimates using the powerlaw package with different modifications. Left: Original package version with biased estimates (Eqs. A.4–A.5) and fixed Xmin; center: Unbiased estimates with fixed Xmin; and right: Unbiased estimates with fitted Xmin. In all three cases a random variable was simulated with single-PL distribution according to Eq. A.1.
In the text