FAST Observations of an Extremely Active Episode of FRB 20201124A. II. Energy Distribution

The data recorded by FAST are in fits format with 4096 frequency channels and a time resolution of 49.152 μs. We first perform de-dispersion processing on the time-frequency data according to different DM values to construct new time-DM data. We then use the trained object detection model to identify the arrival time and DM of the bursts in the data. According to the DM value and burst arrival time, the data on pulse candidates are cut from the original data, and a binary classification neural network is utilized to determine whether they are true bursts or false signals. We also searched for bursts in dedispersed time series using a matched filtering algorithm as implemented by the presto program (Ransom 2001). We combined the bursts detected by the two search methods as the final sample.

Since DM is the integral of electron density along the path, we assume that the DM does not change significantly within an hour. We first use a uniform DM value of 413 pc cm⁻³ to select the pulses detected in one day from the data. Under 8 times down-sampling of time, bursts with a signal-to-noise ratio (S/N) greater than 10 are selected. We then change the DM trials ranging from −5 to 5 pc cm⁻³ from the central value for de-dispersion to maximize the structure of all the burst profiles, and reapply this DM value to all bursts on this day.

For each burst, we estimate the peak flux density with the radiometer equation

$$\begin{eqnarray}&&{S}_{\nu }=\displaystyle \frac{\mathrm{SNR}\times {T}_{\mathrm{sys}}}{G\sqrt{{n}_{p}\times \mathrm{BW}\times {t}_{\mathrm{samp}}}},\end{eqnarray} \tag{ 1 }$$

where T_sys is the system temperature and G is the gain of the telescope, which is modeled as a function of the zenith angle and observation frequency in Jiang et al. (2020), n_p = 2 is the number of polarizations summed, t_samp is the sampling time and BW is the bandwidth used for calibration.

Repeating FRBs tend to have a narrow emission bandwidth compared with apparent non-repeating events (Pleunis et al. 2021). The bandwidth may be related to the telescope sensitivity. For these narrow band bursts, usually the energy is calculated using the flux density defined by the full bandwidth to be multiplied by the full bandwidth. Alternatively, the energy can be calculated using the flux defined with the narrow band in which the burst is detected and multiplied by this narrow band. We roughly estimate the emission bandwidth of each burst by dividing the whole bandpass into multiple 50 MHz sub-bands and identifying the sub-bands containing burst emission. If the calculations are performed correctly, the two methods should result in the same results. In order to check this, we calculate S_ν using two ways. The first is to assume that the burst energy is spread in the full bandwidth of the telescope BW = 500 MHz for all the bursts. The second is to use the true BW for each individual burst. When different bandwidths are adopted, the noise level and the S/N of a pulse change. Therefore, when utilizing the observed pulse bandwidth for flux calibration, we recalculated the S/N and the T_sys and G of the corresponding frequency bandwidth.

During the 17 day observing campaign, we detected 35, 72, 232 and 542 bursts from FRB 20201124A in the first four days respectively (i.e., from September 25 to 28), for a total of 881 bursts. The bursts are defined as significant signals separated by dips down to the noise level in this paper. ¹³ Properties of these bursts are listed in the Appendix Table A1. The observation and detection are depicted in Figure 1. The cumulative distribution is plotted in the linear-logarithmic scale. The linear line in the figure suggests that the accumulated number increases exponentially with time. Interestingly, the bursts suddenly disappeared the next day, which was unexpected. The peak burst rate on the third and fourth day reached 232 hr⁻¹ and 542 hr⁻¹, respectively, which exceeded all previously reported burst rates from this and other repeating FRBs, including 122 hr⁻¹ (Li et al. 2021) and 218 hr⁻¹ (Jahns et al. 2022) for FRB 20121102A, 4.5 hr⁻¹ for FRB 20190520B (Niu et al. 2022b) and 45.8 hr⁻¹ for FRB 20201124A itself during the previous active session (Xu et al. 2021).

Zoom In Zoom Out Reset image size

Waiting time is defined as the difference between the arrival time of two adjacent bursts. All the waiting times were calculated in the same observation session to avoid the observation gap of about 23 hr. The waiting time distribution of FRB 20201124A is plotted in Figure 2, which can be well-fitted by two log-normal functions peaking at 10.05 s and 51.22 ms, respectively. Note that the waiting time distribution of FRB 20201124A in the previous active period was also a bimodal distribution, but the characteristic timescale of the second (long-duration) peak was 135 s (Xu et al. 2021). Our observation suggests that this characteristic waiting time is not universal within one source, but rather depends on the activity level of the source. The bimodal form of waiting time has also been found in FRB 20121102A (Aggarwal et al. 2021; Li et al. 2021; Hewitt et al. 2022). The timescale of the second peak varies also significantly for telescopes with different sensitivities, with shorter waiting times for more sensitive telescopes (which detect more bursts).

Zoom In Zoom Out Reset image size

We calculate the isotropic equivalent burst energy following the equation

$$\begin{eqnarray}&&E={10}^{39}\mathrm{erg}\displaystyle \frac{4\pi }{1+z}{\left(\displaystyle \frac{{D}_{L}}{{10}^{28}\mathrm{cm}}\right)}^{2}\left(\displaystyle \frac{{F}_{\nu }}{\mathrm{Jy}\cdot \mathrm{ms}}\right)\left(\displaystyle \frac{{\rm{\Delta }}\nu }{\mathrm{GHz}}\right),\end{eqnarray} \tag{ 2 }$$

where D_L = 453.3 Mpc is the luminosity distance of FRB 20201124A corresponding to the redshift z = 0.09795 acquired in Xu et al. (2021), F_ν = S_ν × W_eq is the specific fluence, S_ν is the peak specific flux and Δν is the bandwidth of each pulse. The adoption of Δν is more relevant for narrow band FRBs, which are typically the case for the bursts in repeaters. As discussed in Section 2.3, we adopt two choices of BW to calculate S_ν and compare the derived E values.

Figure 3(A) features the specific fluence (F_ν ) and the energy distribution of bursts calculated with Δν = 500 MHz and the observed bandwidth for individual bursts. The specific fluence derived using the observed narrow bandwidth is greater than that derived assuming the full bandwidth. This is natural because the definition of F_ν is the average over frequency channels. Since the energy contained in the observed bandwidth is similar to that in the full bandwidth but the observed bandwidth is narrower than the full bandwidth, the average value for the full bandwidth is smaller. Figure 3(B) presents the energy distribution of the two methods after multiplying by the respective bandwidth. It can be seen that the derived energy distributions between the two methods appear to be similar to each other, as is expected. The p-value from the Kolmogorov-Smirnov (K-S) test of the energy distribution calculated with two types of bandwidth is 0.35, which cannot reject the hypothesis that the two energy distributions are consistent with each other. Therefore, we conclude that the energy calculation is not significantly affected by the adoption of the bandwidth, and the normal method using the full bandwidth to calculate the energy is reliable.

Zoom In Zoom Out Reset image size

Li et al. (2021) introduced the detection of 1652 pulses from FRB 20121102A with FAST and proposed a bimodal energy distribution. The calculation there follows Equation (9) from Zhang (2018), which differs from Equation (3) by a constant ν_c /Δν, if we calibrate the fluence both with the full bandwidth. ¹⁴ Aggarwal (2021) claimed that he recalculated the energy using the bandwidth in Table 1 of Li et al. (2021), and the resulting distribution did not show a significant bimodality. As we have demonstrated in Figure 3(A), the specific fluence varies greatly with the bandwidth selection. The calculation in Aggarwal (2021) simply replaced ν_c with the apparent observed bandwidth for individual bursts, but still used the fluence values derived from the full bandwidth to do the calculation. Therefore the calculations of Aggarwal (2021) resulted in incorrect results. The bimodal energy distribution for FRB 20121102A as claimed by Li et al. (2021) remains intact.

Table 1. Energy Budget of Three Repeating FRBs

Name	Ob Days	Ob Times	Total Observed Energy a	Averaged Energy b	Total Energy c
	(day)	Tobs (hr)	Eradio (erg)	${\bar{L}}_{\mathrm{radio}}$ (erg s−1)	Ebursts (erg)
FRB 20201124A	4	4	1.60 × 1041	1.11 × 1037	3.85 × 1046
FRB 20201124A-0928	1	1	1.02 × 1041	2.84 × 1037	2.46 × 1046
FRB 20121102A	47	59.5	3.41 × 1041	1.59 × 1036	6.47 × 1046
FRB 20190520B	11	18.5	1.10 × 1040	1.65 × 1035	1.56 × 1045

Notes.

^a Sum of the observed isotropic radio energies of all bursts. ^b Defined as E_radio/T_obs. ^c ${E}_{\mathrm{bursts}}\,=\,{E}_{\mathrm{radio}}\times \left[{(\zeta \,=\,\mathrm{ObTimes}/\mathrm{ObDays})}^{-1}\times {(\eta \,=\,{10}^{-5})}^{-1}\times ({F}_{b}\,=\,0.1),.\right]$

Download table as:  ASCIITypeset image

Energy/luminosity function can provide important hints on the emission mechanism of the source. As we demonstrated in the previous subsection, the use of limited bandwidth or full bandwidth does not change the energy distribution. Therefore, we uniformly use Equation (2) and the full bandwidth fluence for energy calculation in the following analysis.

As affirmed in Figure 4(C), the energy distribution of FRB 20201124A also presents a clear bimodal distribution, which can be well-fitted by two log-normal functions peaking at 2.27 × 10³⁷ erg and 2.28 × 10³⁸ erg, respectively. This may indicate the existence of more than one emission mechanism to generate FRBs. The distribution of FRB 20121102A could be also fitted by a similar double log-normal distribution (Extended Data Figure 5 in Li et al. 2021), but the peak energies are not the same for the two FRB sources.

Zoom In Zoom Out Reset image size

For a power law energy/luminosity distribution ${dN}/{dE}\,={{AE}}^{\alpha }$, the cumulative energy distribution can be calculated as

$$\begin{eqnarray}\begin{array}{rcl}N(\gt E) & = & {\displaystyle \int }_{E}^{{E}_{\max }}{dN}={\displaystyle \int }_{E}^{{E}_{\max }}{{AE}}^{\alpha }{dE}\\ & = & \displaystyle \frac{A}{-\alpha -1}\left[{E}^{\alpha +1}-{E}_{\max }^{\alpha +1}\right]=\displaystyle \frac{A}{-\hat{\alpha }}\left[{E}^{\hat{\alpha }}-{E}_{\max }^{\hat{\alpha }}\right],\end{array}\end{eqnarray} \tag{ 3 }$$

where ${E}_{\max }$ is the maximum energy for a certain sample, and $\hat{\alpha }=\alpha +1$. One can see that, in general, the cumulative energy distribution is also roughly a power law with index $\hat{\alpha }$ being $\hat{\alpha }\lt 0$. Different types of astrophysical sources have different $\hat{\alpha }$ values. For example, solar flare events have $-0.9\lt \hat{\alpha }\lt -0.5$ (Maehara et al. 2015) and giant pulses of the Crab pulsar have $-2.8\lt \hat{\alpha }\lt -1.1$ (Mickaliger et al. 2012; Lyu et al. 2021). Different repeating FRBs are observed to have different $\hat{\alpha }$ values: e.g., $-1.8\lt \hat{\alpha }\lt -0.7$ for FRB 20121102A (Law et al. 2017; Gourdji et al. 2019; Cruces et al. 2021), and $\hat{\alpha }\sim -1.3$ for FRB 20180916B (CHIME/FRB Collaboration et al. 2020). For FRB 20201124A, $\hat{\alpha }\sim -1.2$ was measured for the uGMRT bursts (Marthi et al. 2021) and $\hat{\alpha }\sim -3.6$ was measured for the CHIME bursts (Lanman et al. 2022). A broken power law was also used to model the cumulative energy distribution of FRB 20121102A, and the indexes $\hat{\alpha }=-1.4$ and $\hat{\beta }=-1.8$ around the turning point 2.3 × 10³⁷ erg (Aggarwal et al. 2021; Hewitt et al. 2022). The bursts in our FAST sample also require a two-segment power law, but the turning point is not obvious. Therefore, we apply an exponentially connected broken-power law function (which is called the Band function proposed by Band et al. (1993) to describe the gamma-ray spectra of gamma-ray bursts), which reads

$$\begin{eqnarray}N(\gt E)=\left\{\begin{array}{ll}{{AE}}^{\hat{\alpha }}\,{e}^{\left(-E/{E}_{0}\right)}\quad & E\leqslant (\hat{\alpha }-\hat{\beta }){E}_{0}\\ {{AE}}^{\hat{\beta }}{\left[\displaystyle \frac{(\hat{\alpha }-\hat{\beta }){E}_{0}}{e}\right]}^{\hat{\alpha }-\hat{\beta }}\quad & E\geqslant (\hat{\alpha }-\hat{\beta }){E}_{0},\end{array}\right.\end{eqnarray} \tag{ 4 }$$

where $\hat{\alpha }=\alpha +1$ and $\hat{\beta }=\beta +1$ are the power law indices of the lower and higher cumulative energy distributions, and E₀ is a characteristic energy. Note that strictly speaking the cumulative energy distribution of two log-normal functions should have a complicated shape. The simple power law function is insufficient. We therefore adopt the next simplest function, i.e., a smoothly joined two-segment broken power law function (the Band function) to fit it. As shown in Figure 4(A), this function fits the data reasonably well. In principle, the turnover in the low-energy end of the low-energy log-normal component would demand an additional break in the function to fit the cumulative distribution and indeed there is a deviation at the low-energy end between the data and fitting function. In any case, when considering possible selection effect near the sensitivity threshold (which would compensate the deficit of low-energy bursts from the data), we believe that the Band function presents a reasonable description of the cumulative distribution function.

We exclude the bursts that are below the 90% detection threshold (E_th = 3.6 × 10³⁶ erg), which make up 1% of the total sample. We then randomly select 90% of the bursts in the sample to generate a cumulative energy distribution. We repeat this process 1000 times and fit these 1000 cumulative energy distributions with a power law function and the Band function. The central limit theorem makes the fitting parameters converge to a Gaussian distribution. We use the mean and standard deviation of the Gaussian distribution as the standard value and error of the fitting parameters. Figure 4(A) plots the fitting results of the cumulative energy distribution. The power law index is $\hat{\alpha }=-0.08\pm 0.01$. The index values in the Band function are $\hat{\alpha }=-0.22\pm 0.01$ and $\hat{\beta }=-3.27\pm 0.34$, with the turning point at (1.1 ± 0.2) × 10³⁹ erg. Note that the power law indices and the break energy are all different from those derived from the CHIME bursts. Considering the uncertainties, the index of the higher energy end is consistent with the CHIME result.

We repeat the fitting process with a varying energy threshold. The index $\hat{\alpha }$ of a single-power law function does not converge to a single value in panel B of Figure 4. This is another indication that a single power law distribution is inadequate to describe the cumulative energy distribution of this FRB source.

FRB 20121102A, FRB 20190520B and FRB 20201124A are three repeating FRBs which have been extensively studied by FAST. FRB 20121102A is the first known repeating FRB (Spitler et al. 2014, 2016). FRB 20190520B is a repeating FRB detected by Commensal Radio Astronomy FAST Survey (CRAFTS) (Li et al. 2019; Niu et al. 2022b). Both FRB 20121102A and FRB 20190520B are found to be associated with a compact persistent radio source (PRS) (Marcote et al. 2017; Chatterjee et al. 2017; Niu et al. 2022b). FRB 20201124A is also associated with persistent radio emission, but it is not a point source and is rather extended in space (Piro et al. 2021; Ravi et al. 2022). The FAST-measured energy distributions of these three FRB sources can be directly compared against each other to reduce the possible selection biases caused by using different telescope systems. Figure 5 displays the energy distributions of the three FRBs in the form of both power density function (PDF) and cumulative density function (CDF). Both the PDF and CDF show that the energy distributions of the three FRB sources are not consistent.

Zoom In Zoom Out Reset image size

Table 1 lists the observed isotropic radio energy and the inferred burst energy of the three FRBs. The total observed energy E_radio is the sum of all the observed burst energies in the radio band, which is corrected to the total energy released during the bursts via Equation (5).

$$\begin{eqnarray}&&{E}_{\mathrm{bursts}}={E}_{\mathrm{radio}}\times {F}_{b}\times {\eta }^{-1}\times {\zeta }^{-1},\end{eqnarray} \tag{ 5 }$$

where E_bursts is the total energy emitted by the FRB source during bursts, F_b is the beaming factor, η is the radio efficiency which is normalized to ∼10⁻⁵ (similar order to FRB 20200428) and ζ is the observation duty cycle. For example, the observations of FRB 20201124A this time were carried out in four sessions over 4 days with 1 hr per session. The duty cycle ζ is therefore estimated as (4 × 1 hr)/(4 × 24 hr) = 1/24.

Compared with other repeaters, FRB 20201124A has the largest averaged energy, especially on September 28. Taking the nominal values of η = 10⁻⁵ and F_b = 0.1, the total burst energy released on September 28 reached $(2.46\times {10}^{46}\ \mathrm{erg}){\eta }_{-5}^{-1}{F}_{b.-1}$. Compared with the total dipolar magnetic energy of a magnetar ${E}_{\mathrm{mag}}=(1/6){B}_{p}^{2}{R}^{3}\sim (1.7\times {10}^{47}\ \mathrm{erg}){B}_{p,15}^{2}{R}_{6}^{3}$, the burst energy emitted on this day already exceeded $14.3 \% \ {\eta }_{-5}^{-1}{F}_{b,-1}$ of the available magnetar energy. This raises an even more severe energy budget problem compared with FRB 20121102A (Li et al. 2021). This requires that the radio emission efficiency should be quite high for the magnetar model or that the magnetar engine does not apply. For magnetar models, since the synchrotron maser models invoking relativistic shocks are quite inefficient, the data favor magnetospheric models for FRBs, which can have a higher efficiency as observed in radio pulsars (Kumar et al. 2017; Yang & Zhang 2018; Lu et al. 2020; Zhang 2022). The synchrotron maser model may be still relevant if the central engine is not a magnetar, but rather an accreting black hole system (e.g., Sridhar et al. 2021), whose total energy budget is not limited by the total magnetic energy of a magnetar.

To quantify the specific morphological characteristics of the bursts, we differentiate two kinds of pulse profile classification schemes. The first classification scheme is based on a judgment of pulse flux threshold, taking the pulse profile above the detection threshold as a separated burst. For conservative considerations, we take the detection threshold using a 3σ noise level of the baseline due to the complicated effect on the actual sensitivity, e.g., the bandwidth limited structure of the bursts and radio frequency interference (RFI) events. This definition is also the pulse definition used in previous work, and is the same as the definition used for FRB 20121102A (Li et al. 2021) and FRB 20190520B (Niu et al. 2022b). In the following text, we call this definition Combined Bursts Definition (CBD).

For another method, we try to differentiate the bursts by referring to the morphological clustering characteristics of the pulses in the dynamic spectrum, which means that pulses that can be distinguished in the time-frequency diagram are treated as different bursts, like Figure 6. We call this definition Separated Bursts Definition (SBD). As a result, the event in Figure 6 contains two bursts under SBD but one burst under CBD. All bursts under SBD are artificially isolated and flux calibration is performed independently.

Zoom In Zoom Out Reset image size

Figure 7(A) depicts the waiting time distribution of bursts with SBD, which is still a bimodal distribution and can be modeled with two log-normal functions peaking at 13.56 ms and 10.05 s, respectively. The second peak time around 10 s is consistent with the waiting time distribution of CBD. However, the first peak time is significantly reduced. The distributions of specific fluence and energy exhibit distinct distributions under the two burst definitions. We apply the K-S test to check whether the energies calculated from the two burst definition schemes belong to the same distribution. The p-value is 0.2%, meaning we could rule out the hypothesis that they were from the same distribution. We also use the Band function to fit the cumulative energy distribution of the bursts with SBD. The index values in the Band function are $\hat{\alpha }=-0.18\pm 0.01$ and $\hat{\beta }=-3.27\pm 0.16$. The index in the high energy regime is the same with CBD, while the index in the lower energy regime is slightly smaller than that of CBD. The turning point is (6.4 ± 0.1) × 10³⁸ erg, which is smaller than that for CBD. While we still do not know the true nature of FRBs, different definitions of bursts may lead to somewhat different analysis results. We caution that when studying burst energy distributions of FRBs, one needs to specify a burst definition scheme.

Zoom In Zoom Out Reset image size